4 * Copyright 2000 Compaq Computer Corporation.
5 * Copying or modifying this code for any purpose is permitted,
6 * provided that this copyright notice is preserved in its entirety
7 * in all copies or modifications.
8 * COMPAQ COMPUTER CORPORATION MAKES NO WARRANTIES, EXPRESSED OR
9 * IMPLIED, AS TO THE USEFULNESS OR CORRECTNESS OF THIS CODE OR
12 * Adapted from cmu_recognizer.c by Jay Kistler.
14 * Where is the CMU copyright???? Gotta track it down - Jim Gettys
16 * Credit to Dean Rubine, Jim Kempf, and Ari Rapkin.
20 #include <sys/types.h>
28 #include <hre_internal.h>
33 #include "li_recognizer.h"
34 #include "li_recognizer_internal.h"
41 #define LI_MAGIC 0xACCBADDD
43 #define CHECK_LI_MAGIC(_a) \
44 ((_a) != NULL && ((li_recognizer*)(_a))->li_magic == LI_MAGIC)
47 static void lialg_initialize(rClassifier *);
48 static int lialg_read_classifier_digest(rClassifier *);
49 static int lialg_canonicalize_examples(rClassifier *);
50 static char *lialg_recognize_stroke(rClassifier *, point_list *);
53 char* li_err_msg = NULL;
54 char _zdebug_flag[128];
57 // This is standard - defined in <stdlib.h>
58 #define bcopy(s1,s2,n) memcpy(s2,s1,n)
61 #if 0 /* was #ifdef mips*/
62 char *strdup(char* from) {
64 int len = strlen(from) + 1;
66 /* to = (char *) safe_malloc( len * sizeof(char) );*/
67 to = allocate(len, char);
68 memcpy(to, from, len);
74 /*Freeing classifier*/
77 free_rClassifier(rClassifier* rc);
84 add_example(point_list* l,int npts,pen_point* pts)
86 pen_point* lpts = make_pen_point_array(npts);
87 /* point_list* p = (point_list*)safe_malloc(sizeof(point_list));*/
88 point_list *p = allocate(1, point_list);
92 p->next = l; /*Order doesn't matter, so we stick on end.*/
96 bcopy(pts,lpts,npts * sizeof(pen_point));
103 delete_examples(point_list* l)
107 for( ; l != NULL; l = p ) {
119 * recognize_internal-Form Vector, use Classifier to classify, return char.
123 recognize_internal(rClassifier* rec,pen_stroke* str,int* rconf)
126 point_list *stroke = NULL;
128 stroke = add_example(NULL, str->ps_npts, str->ps_pts);
129 if (stroke == NULL) return(NULL);
131 res = lialg_recognize_stroke(rec, stroke);
133 delete_examples(stroke);
138 * file_path-Construct pathname, check for proper extension.
142 file_path(char* dir,char* filename,char* pathname)
146 /*Check for proper extension on file name.*/
148 dot = strrchr(filename,'.');
154 /*Determine whether a gesture or character classifier.*/
156 if( strcmp(dot,LI_CLASSIFIER_EXTENSION) != 0 ) {
160 /*Concatenate directory and filename into pathname.*/
162 strcpy(pathname,dir);
163 strcat(pathname,"/");
164 strcat(pathname,filename);
169 /*read_classifier_points-Read points so classifier can be extended.*/
172 read_classifier_points(FILE* fd,int nclss,point_list** ex,char** cnames)
177 char* names[MAXSCLASSES];
178 point_list* examples[MAXSCLASSES];
184 for( i = 0; i < MAXSCLASSES; i++ ) {
192 /* fprintf(stderr, "Classes: [ "); */
194 for( k = 0; k < nclss; k++ ) {
196 /*Read class name and number of examples.*/
198 if( fscanf(fd,"%d %s",&nex,buf) != 2 ) {
199 printf("%s *FAILED* - line: %d\n", __FUNCTION__, __LINE__);
205 names[k] = strdup(buf);
207 /* fprintf(stderr, "%s ", buf); */
211 for( i = 0; i < nex; i++ ) {
213 /*Read number of points.*/
215 if( fscanf(fd,"%d",&npts) != 1 ) {
216 printf("%s *FAILED* - line: %d\n", __FUNCTION__, __LINE__);
217 goto unallocate; /*Boy would I like exceptions!*/
220 /*Allocate array for points.*/
222 if( (pts = make_pen_point_array(npts)) == NULL ) {
223 printf("%s *FAILED* - line: %d\n", __FUNCTION__, __LINE__);
229 for( j = 0; j < npts; j++ ) {
232 if( fscanf(fd,"%d %d",&x,&y) != 2 ) {
233 delete_pen_point_array(pts);
234 printf("%s *FAILED* - line: %d\n", __FUNCTION__, __LINE__);
235 printf("class = %d/%d/%s, ex = %d/%d, pt: %d/%d\n",
236 k, nclss, names[k], i, nex, j, npts);
237 for (jj = 0; jj < j; jj++) {
238 printf("pts[%d] = %d/%d\n", jj, pts[jj].x, pts[jj].y);
248 if( (examples[k] = add_example(examples[k],npts,pts)) == NULL ) {
249 delete_pen_point_array(pts);
250 printf("%s *FAILED* - line: %d\n", __FUNCTION__, __LINE__);
254 delete_pen_point_array(pts);
259 /* ari -- end of list of classes */
260 /* fprintf(stderr, "]\n"); */
262 /*Transfer to recognizer.*/
264 bcopy(examples,ex,sizeof(examples));
265 bcopy(names,cnames,sizeof(names));
269 /*Error. Deallocate memory and return.*/
273 for( ; k >= 0; k-- ) {
274 delete_examples(examples[k]);
278 error("Error in reading example points from classifier file");
282 /*read_classifier-Read a classifier file.*/
284 static int read_classifier(FILE* fd,rClassifier* rc)
290 /*Read in classifier file.*/
292 if( (sc = sRead(fd)) == NULL ) {
296 /*Read in the example points, so classifier can be extended.*/
298 if( read_classifier_points(fd,sc->nclasses,rc->ex,rc->cnames) != 0 ) {
303 /*Transfer sClassifier to the rClassifier*/
311 * Extension Functions
314 /* getClasses and clearState are by Ari */
317 recognizer_getClasses (recognizer r, char ***list, int *nc)
324 rec = (li_recognizer*)r->recognizer_specific;
326 /*Check for LI recognizer.*/
328 if( !CHECK_LI_MAGIC(rec) ) {
329 li_err_msg = "Not a LI recognizer";
334 *nc = nclasses = sc->nclasses;
335 /* ret = (char **) safe_malloc (nclasses * sizeof(char*));*/
336 ret = allocate(nclasses, char*);
338 for (i = 0; i < nclasses; i++) {
339 ret[i] = rec->li_rc.cnames[i]; /* only the 1st char of the cname */
346 recognizer_clearState (recognizer r)
348 /*This operation isn't supported by the LI recognizer.*/
350 li_err_msg = "Clearing state is not supported by the LI recognizer";
355 static bool isa_li(recognizer r)
356 { return(CHECK_LI_MAGIC(r)); }
359 recognizer_train(recognizer r,rc* rec_xt,u_int nstrokes,
360 pen_stroke* strokes,rec_element* re,
363 /*This operation isn't supported by the LI recognizer.*/
365 li_err_msg = "Training is not supported by the LI recognizer";
371 li_recognizer_get_example (recognizer r,
378 li_recognizer *rec = (li_recognizer*)r->recognizer_specific;
379 sClassifier sc = rec->li_rc.sc;
382 if( !CHECK_LI_MAGIC(rec) ) {
383 li_err_msg = "Not a LI recognizer";
386 if (class > sc->nclasses)
388 pl = rec->li_rc.canonex[class];
389 while (instance && pl)
396 *name = rec->li_rc.cnames[class];
407 /*li_recognizer_load-Load a classifier file.*/
409 static int li_recognizer_load(recognizer r,char* dir,char* filename)
416 rec = (li_recognizer*)r->recognizer_specific;
418 /*Make sure recognizer's OK*/
420 if( !CHECK_LI_MAGIC(rec) ) {
421 li_err_msg = "Not a LI recognizer";
427 /*Check parameters.*/
429 if( filename == NULL ) {
430 li_err_msg = "Invalid parameters";
434 /*We let the directory be null.*/
436 if( dir == NULL || (int)strlen(dir) <= 0 ) {
440 /*Make full pathname and check filename*/
442 /* pathname = (char*)safe_malloc(strlen(dir) + strlen(filename) + 2)); */
444 pathname = allocate(strlen(dir) + strlen(filename) + 2, char);
445 if( file_path(dir,filename,pathname) == -1 ) {
447 li_err_msg = "Not a LI recognizer classifier file";
451 /* Try to short-circuit the full classifier-file processing. */
452 rc->file_name = pathname;
453 if (lialg_read_classifier_digest(rc) == 0)
455 rc->file_name = NULL;
459 if( (fd = fopen(pathname,"r")) == NULL ) {
461 li_err_msg = "Can't open classifier file";
463 /* fprintf(stderr, "Trying to open %s.\n", pathname); */
468 /*If rClassifier is OK, then delete it first.*/
470 if( rc->file_name != NULL ) {
471 free_rClassifier(rc);
476 if( read_classifier(fd,rc) < 0 ) {
485 /*Add classifier name.*/
487 rc->file_name = pathname;
489 /* Canonicalize examples. */
490 if (lialg_canonicalize_examples(rc) != 0) {
492 rc->file_name = NULL;
499 /*li_recognizer_save-Save a classifier file.*/
501 static int li_recognizer_save(recognizer r,char* dir,char* filename)
503 /*This operation isn't supported by the LI recognizer.*/
505 li_err_msg = "Saving is not supported by the LI recognizer";
511 li_recognizer_load_dictionary(recognizer rec,char* directory,char* name)
513 /*This operation isn't supported by the LI recognizer.*/
515 li_err_msg = "Dictionaries are not supported by the LI recognizer";
522 li_recognizer_save_dictionary(recognizer rec,
527 /*This operation isn't supported by the LI recognizer.*/
529 li_err_msg = "Dictionaries are not supported by the LI recognizer";
536 li_recognizer_free_dictionary(recognizer rec,wordset dict)
538 /*This operation isn't supported by the LI recognizer.*/
540 li_err_msg = "Dictionaries are not supported by the LI recognizer";
547 li_recognizer_add_to_dictionary(recognizer rec,letterset* word,wordset dict)
549 /*This operation isn't supported by the LI recognizer.*/
551 li_err_msg = "Dictionaries are not supported by the LI recognizer";
558 li_recognizer_delete_from_dictionary(recognizer rec,
562 /*This operation isn't supported by the LI recognizer.*/
564 li_err_msg = "Dictionaries are not supported by the LI recognizer";
571 li_recognizer_error(recognizer rec)
573 char* ret = li_err_msg;
575 /*Check for LI recognizer.*/
577 if( !CHECK_LI_MAGIC(rec->recognizer_specific) ) {
578 li_err_msg = "Not a LI recognizer";
588 li_recognizer_clear(recognizer r,bool delete_points_p)
592 rec = (li_recognizer*)r->recognizer_specific;
594 /*Check for LI recognizer.*/
596 if( !CHECK_LI_MAGIC(rec) ) {
597 li_err_msg = "Not a LI recognizer";
605 li_recognizer_set_context(recognizer r,rc* rec_xt)
607 /*This operation isn't supported by the LI recognizer.*/
609 li_err_msg = "Contexts are not supported by the LI recognizer";
615 li_recognizer_get_context(recognizer r)
617 /*This operation isn't supported by the LI recognizer.*/
619 li_err_msg = "Contexts are not supported by the LI recognizer";
625 li_recognizer_get_buffer(recognizer r, u_int* nstrokes,pen_stroke** strokes)
627 /*This operation isn't supported by the LI recognizer.*/
629 li_err_msg = "Buffer get/set are not supported by the LI recognizer";
635 li_recognizer_set_buffer(recognizer r,u_int nstrokes,pen_stroke* strokes)
637 /*This operation isn't supported by the LI recognizer.*/
639 li_err_msg = "Buffer get/set are not supported by the LI recognizer";
645 li_recognizer_translate(recognizer r,
650 rec_alternative** ret)
657 rec = (li_recognizer*)r->recognizer_specific;
662 /*Check for LI recognizer.*/
664 if( !CHECK_LI_MAGIC(rec) ) {
665 li_err_msg = "Not a LI recognizer";
671 /*Check for valid parameters.*/
673 li_err_msg = "Invalid parameters: ncs";
677 li_err_msg = "Invalid parameters: tps";
681 li_err_msg = "Invalid parameters: nret";
685 li_err_msg = "Invalid parameters: ret";
689 /* if( ncs < 1 || tps == NULL || nret == NULL || ret == NULL) {
690 li_err_msg = "Invalid parameters";
695 /*Check for null classifier. It must have at least one.*/
697 if( rec->li_rc.sc == NULL ) {
698 li_err_msg = "No classifier";
704 * Go through the stroke array and recognize. Since this is a single
705 * stroke recognizer, each stroke is treated as a separate
706 * character or gesture. We allow only characters or gestures
707 * to be recognized at one time, since otherwise, handling
708 * the display of segmentation would be difficult.
710 clss = recognize_internal(rc,tps,&conf);
713 li_err_msg = "unrecognized character";
727 li_recognizer_get_extension_functions(recognizer rec)
731 /*Check for LI recognizer.*/
733 if( !CHECK_LI_MAGIC(rec->recognizer_specific) ) {
734 li_err_msg = "Not a LI recognizer";
738 ret = make_rec_fn_array(LI_NUM_EX_FNS);
740 /* ari -- clearState & getClasses are mine */
741 ret[LI_GET_CLASSES] = (rec_fn)recognizer_getClasses;
742 ret[LI_CLEAR] = (rec_fn)recognizer_clearState;
743 ret[LI_ISA_LI] = (rec_fn)isa_li;
744 ret[LI_TRAIN] = (rec_fn)recognizer_train;
750 li_recognizer_get_gesture_names(recognizer r)
752 /*This operation isn't supported by the LI recognizer.*/
754 li_err_msg = "Gestures are not supported by the LI recognizer";
760 li_recognizer_set_gesture_action(recognizer r,
765 /*This operation isn't supported by the LI recognizer.*/
767 li_err_msg = "Gestures are not supported by the LI recognizer";
776 /*RECOGNIZER_INITIALIZE-Initialize the recognizer.*/
778 /* note from ari: this expands via pre-processor to
780 * recognizer __recognizer_internal_initialize(rec_info* ri)
783 RECOGNIZER_INITIALIZE(ri)
789 /*Check that locale matches.*/
791 if( strcmp(ri->ri_locale,LI_SUPPORTED_LOCALE) != 0 ) {
792 li_err_msg = "Not a supported locale";
793 fprintf(stderr, "Locale error.\n");
800 * Check that character sets match. Note that this is only approximate,
801 * since the classifier file will have more information.
804 if( ri->ri_subset != NULL ) {
805 for(i = 0; ri->ri_subset[i] != NULL; i++ ) {
807 if( strcmp(ri->ri_subset[i],UPPERCASE) != 0 &&
808 strcmp(ri->ri_subset[i],LOWERCASE) != 0 &&
809 strcmp(ri->ri_subset[i],DIGITS) != 0 &&
810 strcmp(ri->ri_subset[i],GESTURE) != 0 ) {
811 li_err_msg = "Not a supported character set";
812 fprintf(stderr, "charset error.\n");
820 r = make_recognizer(ri);
821 /*fprintf(stderr, "past make_recognizer.\n");*/
824 li_err_msg = "Can't allocate storage";
829 /*Make a LI recognizer structure.*/
832 /* rec = (li_recognizer*)safe_malloc(sizeof(li_recognizer))) == NULL ); */
834 rec = allocate(1, li_recognizer);
836 r->recognizer_specific = rec;
838 rec->li_rc.file_name = NULL;
839 rec->li_rc.sc = NULL;
841 /*Initialize the recognizer struct.*/
843 r->recognizer_load_state = li_recognizer_load;
844 r->recognizer_save_state = li_recognizer_save;
845 r->recognizer_load_dictionary = li_recognizer_load_dictionary;
846 r->recognizer_save_dictionary = li_recognizer_save_dictionary;
847 r->recognizer_free_dictionary = li_recognizer_free_dictionary;
848 r->recognizer_add_to_dictionary = li_recognizer_add_to_dictionary;
849 r->recognizer_delete_from_dictionary = li_recognizer_delete_from_dictionary;
850 r->recognizer_error = li_recognizer_error;
851 r->recognizer_translate = li_recognizer_translate;
852 r->recognizer_get_context = li_recognizer_get_context;
853 r->recognizer_set_context = li_recognizer_set_context;
854 r->recognizer_get_buffer = li_recognizer_get_buffer;
855 r->recognizer_set_buffer = li_recognizer_set_buffer;
856 r->recognizer_clear = li_recognizer_clear;
857 r->recognizer_get_extension_functions =
858 li_recognizer_get_extension_functions;
859 r->recognizer_get_gesture_names = li_recognizer_get_gesture_names;
860 r->recognizer_set_gesture_action =
861 li_recognizer_set_gesture_action;
863 /*Initialize LI Magic Number.*/
865 rec->li_magic = LI_MAGIC;
867 /*Initialize rClassifier.*/
869 rec->li_rc.file_name = NULL;
871 for( i = 0; i < MAXSCLASSES; i++ ) {
872 rec->li_rc.ex[i] = NULL;
873 rec->li_rc.cnames[i] = NULL;
876 lialg_initialize(&rec->li_rc);
878 /*Get rid of error message. Not needed here.*/
884 /*free_rClassifier-Free the rClassifier.*/
887 free_rClassifier(rClassifier* rc)
891 if( rc->file_name != NULL) {
895 for( i = 0; rc->ex[i] != NULL; i++) {
896 delete_examples(rc->ex[i]);
900 if(rc->sc != NULL ) {
901 sFreeClassifier(rc->sc);
905 /*RECOGNIZER_FINALIZE-Deallocate the recognizer, finalize.*/
907 RECOGNIZER_FINALIZE(r)
909 li_recognizer* rec = (li_recognizer*)r->recognizer_specific;
911 /*Make sure this is a li_recognizer first*/
913 if( !CHECK_LI_MAGIC(rec) ) {
914 li_err_msg = "Not a LI recognizer";
918 /*Deallocate rClassifier resources.*/
920 free_rClassifier(&(rec->li_rc));
922 /*Deallocate the li_recognizer struct.*/
926 /*Deallocate the recognizer*/
928 delete_recognizer(r);
934 /* **************************************************
936 Implementation of the Li/Yeung recognition algorithm
938 ************************************************** */
940 /*#include <assert.h>*/
942 #define MAXINT 0x7FFFFFFF
946 #include <sys/time.h>
949 /* Ultrix doesn't have these declarations in math.h! */
950 extern double rint(double);
951 extern float expf(float);
955 extern double rint (double);
956 extern float expf (float); /* N.B. exp() appears to be broken on ELX! */
959 #define WORST_SCORE MAXINT
961 /* Dynamic programming parameters */
964 #define MAX_DIST MAXINT
965 #define SIM_THLD 60 /* x 100 */
966 #define DIST_THLD 3200 /* x 100 */
968 /* Low-pass filter parameters -- empirically derived */
969 #define LP_FILTER_WIDTH 6
970 #define LP_FILTER_ITERS 8
971 #define LP_FILTER_THLD 250 /* x 100 */
972 #define LP_FILTER_MIN 5
974 /* Pseudo-extrema parameters -- empirically derived */
975 #define PE_AL_THLD 1500 /* x 100 */
976 #define PE_ATCR_THLD 135 /* x 100 */
978 /* Contour-angle derivation parameters */
982 /* Pre-processing and canonicalization parameters */
983 #define CANONICAL_X 108
984 #define CANONICAL_Y 128
985 #define DIST_SQ_THRESHOLD (3*3) /* copied from fv.h */
986 #define NCANONICAL 50
988 /* Tap-handling parameters */
990 #define TAP_TIME_THLD 150 /* msec */
991 #define TAP_DIST_THLD 75 /* dx * dx + dy * dy */
992 #define TAP_PATHLEN 1000 /* x 100 */
995 /* Overload the time field of the pen_point struct with the chain-code. */
996 #define chaincode time
1000 #define RGN_CONCAVE 1
1002 #define RGN_PSEUDO 3
1005 typedef struct RegionList {
1009 struct RegionList *next;
1013 /* direction-code table; indexed by dx, dy */
1014 static int lialg_dctbl[3][3] = {{1, 0, 7}, {2, 0x7FFFFFFF, 6}, {3, 4, 5}};
1016 /* low-pass filter weights */
1017 static int lialg_lpfwts[2 * LP_FILTER_WIDTH + 1];
1018 static int lialg_lpfconst = -1;
1021 static int lialg_preprocess_stroke(point_list *);
1022 static point_list *lialg_compute_dominant_points(point_list *);
1023 static point_list *lialg_interpolate_points(point_list *);
1024 static void lialg_bresline(pen_point *, pen_point *, point_list *, int *);
1025 static void lialg_compute_chain_code(point_list *);
1026 static void lialg_compute_unit_chain_code(point_list *);
1027 static region_list *lialg_compute_regions(point_list *);
1028 static point_list *lialg_compute_dompts(point_list *, region_list *);
1029 static int *lialg_compute_contour_angle_set(point_list *, region_list *);
1030 static void lialg_score_stroke(point_list *, point_list *, int *, int *);
1031 static int lialg_compute_similarity(point_list *, point_list *);
1032 static int lialg_compute_distance(point_list *, point_list *);
1034 static int lialg_read_classifier_digest(rClassifier *);
1036 static int lialg_canonicalize_examples(rClassifier *);
1037 static int lialg_canonicalize_example_stroke(point_list *);
1038 static int lialg_compute_equipoints(point_list *);
1040 static int lialg_compute_pathlen(point_list *);
1041 static int lialg_compute_pathlen_subset(point_list *, int, int);
1042 static int lialg_filter_points(point_list *);
1043 static int lialg_translate_points(point_list *, int, int, int, int);
1044 static void lialg_get_bounding_box(point_list *, int *, int *, int *, int *);
1045 static void lialg_compute_lpf_parameters();
1046 static int isqrt(int);
1047 static int likeatan(int, int);
1048 static int quadr(int);
1051 /*************************************************************
1053 Core routines for the Li/Yeung recognition algorithm
1055 *************************************************************/
1057 static void lialg_initialize(rClassifier *rec) {
1060 /* Initialize the dompts arrays. */
1061 for (i = 0; i < MAXSCLASSES; i++) {
1062 rec->dompts[i] = NULL;
1068 * Main recognition routine -- called by HRE API.
1070 static char *lialg_recognize_stroke(rClassifier *rec, point_list *stroke) {
1073 point_list *input_dompts = NULL;
1074 char *best_name = NULL;
1075 int best_score = WORST_SCORE;
1077 point_list *curr_dompts = NULL;
1078 /*struct timeval stv, etv;
1081 /* (void)gettimeofday(&stv, NULL);*/
1083 if (stroke->npts < 1) goto done;
1085 /* Check for tap. */
1088 pen_point *startpt = &stroke->pts[0];
1089 pen_point *endpt = &stroke->pts[stroke->npts - 1];
1090 int dt = endpt->time - startpt->time;
1091 int dx = endpt->x - startpt->x;
1092 int dy = endpt->y - startpt->y;
1093 int magsq = dx * dx + dy * dy;
1096 /* First thing is to filter out ``close points.'' */
1097 if (lialg_filter_points(stroke) != 0) return(NULL);
1099 /* Unfortunately, we don't have the actual time that each point */
1100 /* was recorded (i.e., dt is invalid). Hence, we have to use a */
1101 /* heuristic based on total distance and the number of points. */
1102 if (stroke->npts == 1 || lialg_compute_pathlen(stroke) < TAP_PATHLEN)
1106 /* Pre-process input stroke. */
1107 if (lialg_preprocess_stroke(stroke) != 0) goto done;
1109 /* Compute its dominant points. */
1110 input_dompts = lialg_compute_dominant_points(stroke);
1111 if (input_dompts == NULL) goto done;
1113 /* Score input stroke against every class in classifier. */
1114 for (i = 0, curr_name = rec->cnames[i], curr_dompts = rec->dompts[i];
1115 i < MAXSCLASSES && curr_name != NULL && curr_dompts != NULL;
1116 i++, curr_name = rec->cnames[i], curr_dompts = rec->dompts[i]) {
1121 lialg_score_stroke(input_dompts, curr_dompts, &sim, &dist);
1124 if (lidebug && curr_score < DIST_THLD)
1125 fprintf(stderr, "(%s, %d, %d) ", curr_name, sim, dist);
1127 /* Is it the best so far? */
1128 if (curr_score < best_score && curr_score <= DIST_THLD) {
1129 best_score = curr_score;
1130 best_name = curr_name;
1135 fprintf(stderr, "\n");
1141 delete_examples(input_dompts);
1142 /* (void)gettimeofday(&etv, NULL);
1143 elapsed = (1000 * (etv.tv_sec - stv.tv_sec)) + ((etv.tv_usec - stv.tv_usec + 500) / 1000);
1144 fprintf(stderr, "elapsed = %d\n", elapsed);
1150 static int lialg_preprocess_stroke(point_list *points) {
1151 int minx, miny, maxx, maxy, xrange, yrange, scale, xoff, yoff;
1153 /* Filter out points that are too close. */
1154 /* We did this earlier, when we checked for a tap. */
1156 if (lialg_filter_points(points) != 0) return(-1);
1159 /* assert(points->npts > 0);*/
1161 /* Scale up to avoid conversion errors. */
1162 lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
1163 xrange = maxx - minx;
1164 yrange = maxy - miny;
1165 scale = ( ((100 * xrange + CANONICAL_X / 2) / CANONICAL_X) >
1166 ((100 * yrange + CANONICAL_Y / 2) / CANONICAL_Y))
1167 ? (100 * CANONICAL_X + xrange / 2) / xrange
1168 : (100 * CANONICAL_Y + yrange / 2) / yrange;
1169 if (lialg_translate_points(points, minx, miny, scale, scale) != 0)
1172 /* Center the stroke. */
1173 lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
1174 xrange = maxx - minx;
1175 yrange = maxy - miny;
1176 xoff = -((CANONICAL_X - xrange + 1) / 2);
1177 yoff = -((CANONICAL_Y - yrange + 1) / 2);
1178 if (lialg_translate_points(points, xoff, yoff, 100, 100) != 0) return(-1);
1180 /* Store the x and y ranges in the point list. */
1181 xrange = maxx - minx;
1182 yrange = maxy - miny;
1183 points->xrange = xrange;
1184 points->yrange = yrange;
1188 fprintf(stderr, "After pre-processing: %d %d %d %d\n",
1189 minx, miny, maxx, maxy);
1190 for (i = 0; i < points->npts; i++)
1191 fprintf(stderr, " (%d %d)\n",
1192 points->pts[i].x, points->pts[i].y);
1200 static point_list *lialg_compute_dominant_points(point_list *points) {
1201 point_list *ipts = NULL;
1202 region_list *regions = NULL;
1203 point_list *dpts = NULL;
1205 /* Interpolate points. */
1206 ipts = lialg_interpolate_points(points);
1207 if (ipts == NULL) return(NULL);
1210 fprintf(stderr, "After interpolation: %d ipts\n", ipts->npts);
1211 for (j = 0; j < ipts->npts; j++) {
1212 fprintf(stderr, " (%d, %d), %ld\n",
1213 ipts->pts[j].x, ipts->pts[j].y, ipts->pts[j].chaincode);
1218 /* Compute regions. */
1219 regions = lialg_compute_regions(ipts);
1220 /* assert(regions != NULL);*/
1222 /* Compute dominant points. */
1223 dpts = lialg_compute_dompts(ipts, regions);
1226 fprintf(stderr, "Dominant points: ");
1227 for (j = 0; j < dpts->npts; j++) {
1228 fprintf(stderr, "%d %d (%ld) ",
1229 dpts->pts[j].x, dpts->pts[j].y, dpts->pts[j].chaincode);
1231 fprintf(stderr, "\n");
1235 /* Delete region data structure. */
1237 region_list *curr, *next;
1238 for (curr = regions; curr != NULL; ) {
1244 delete_examples(ipts);
1249 /* Input points are assumed to be integer-valued! */
1250 static point_list *lialg_interpolate_points(point_list *points) {
1255 /* Compute an upper-bound on the number of interpolated points. */
1257 for (i = 0; i < (points->npts - 1); i++) {
1258 pen_point *pta = &(points->pts[i]);
1259 pen_point *ptb = &(points->pts[i+1]);
1260 maxpts += abs(pta->x - ptb->x) + abs(pta->y - ptb->y);
1263 /* Allocate an array of the requisite size. */
1264 maxpts += points->npts;
1265 /* newpts = (point_list *)safe_malloc(sizeof(point_list)); */
1266 newpts = allocate(1, point_list);
1267 newpts->pts = make_pen_point_array(maxpts);
1268 if (newpts->pts == NULL) {
1272 newpts->npts = maxpts;
1273 newpts->next = NULL;
1275 /* Interpolate each of the segments. */
1277 for (i = 0; i < (points->npts - 1); i++) {
1278 pen_point *startpt = &(points->pts[i]);
1279 pen_point *endpt = &(points->pts[i+1]);
1281 lialg_bresline(startpt, endpt, newpts, &j);
1283 j--; /* end point gets recorded as start point of next segment! */
1286 /* Add-in last point. */
1287 newpts->pts[j++] = points->pts[points->npts - 1];
1290 /* Compute the chain code for P (the list of points). */
1291 lialg_compute_unit_chain_code(newpts);
1297 /* This implementation is due to Kenny Hoff. */
1298 static void lialg_bresline(pen_point *startpt, pen_point *endpt,
1299 point_list *newpts, int *j) {
1300 int Ax, Ay, Bx, By, dX, dY, Xincr, Yincr;
1307 /* INITIALIZE THE COMPONENTS OF THE ALGORITHM THAT ARE NOT AFFECTED */
1308 /* BY THE SLOPE OR DIRECTION OF THE LINE */
1309 dX = abs(Bx-Ax); /* store the change in X and Y of the line endpoints */
1312 /* DETERMINE "DIRECTIONS" TO INCREMENT X AND Y (REGARDLESS OF DECISION) */
1313 if (Ax > Bx) { Xincr=-1; } else { Xincr=1; } /* which direction in X? */
1314 if (Ay > By) { Yincr=-1; } else { Yincr=1; } /* which direction in Y? */
1316 /* DETERMINE INDEPENDENT VARIABLE (ONE THAT ALWAYS INCREMENTS BY 1 (OR -1) ) */
1317 /* AND INITIATE APPROPRIATE LINE DRAWING ROUTINE (BASED ON FIRST OCTANT */
1318 /* ALWAYS). THE X AND Y'S MAY BE FLIPPED IF Y IS THE INDEPENDENT VARIABLE. */
1319 if (dX >= dY) { /* if X is the independent variable */
1320 int dPr = dY<<1; /* amount to increment decision if right is chosen (always) */
1321 int dPru = dPr - (dX<<1); /* amount to increment decision if up is chosen */
1322 int P = dPr - dX; /* decision variable start value */
1324 /* process each point in the line one at a time (just use dX) */
1325 for (; dX>=0; dX--) {
1326 newpts->pts[*j].x = Ax;
1327 newpts->pts[*j].y = Ay;
1330 if (P > 0) { /* is the pixel going right AND up? */
1331 Ax+=Xincr; /* increment independent variable */
1332 Ay+=Yincr; /* increment dependent variable */
1333 P+=dPru; /* increment decision (for up) */
1335 else { /* is the pixel just going right? */
1336 Ax+=Xincr; /* increment independent variable */
1337 P+=dPr; /* increment decision (for right) */
1341 else { /* if Y is the independent variable */
1342 int dPr = dX<<1; /* amount to increment decision if right is chosen (always) */
1343 int dPru = dPr - (dY<<1); /* amount to increment decision if up is chosen */
1344 int P = dPr - dY; /* decision variable start value */
1346 /* process each point in the line one at a time (just use dY) */
1347 for (; dY>=0; dY--) {
1348 newpts->pts[*j].x = Ax;
1349 newpts->pts[*j].y = Ay;
1352 if (P > 0) { /* is the pixel going up AND right? */
1353 Ax+=Xincr; /* increment dependent variable */
1354 Ay+=Yincr; /* increment independent variable */
1355 P+=dPru; /* increment decision (for up) */
1357 else { /* is the pixel just going up? */
1358 Ay+=Yincr; /* increment independent variable */
1359 P+=dPr; /* increment decision (for right) */
1366 static void lialg_compute_chain_code(point_list *pts) {
1369 for (i = 0; i < (pts->npts - 1); i++) {
1370 pen_point *startpt = &(pts->pts[i]);
1371 pen_point *endpt = &(pts->pts[i+1]);
1372 int dx = endpt->x - startpt->x;
1373 int dy = endpt->y - startpt->y;
1375 int tmp = rint(4.0 * atan2((double)dx, (double)dy) / M_PI);
1376 int dircode = (10 + tmp) % 8;
1378 int tmp = quadr(likeatan(dy, dx));
1379 int dircode = (12 - tmp) % 8;
1381 startpt->chaincode = dircode;
1386 static void lialg_compute_unit_chain_code(point_list *pts) {
1389 for (i = 0; i < (pts->npts - 1); i++) {
1390 pen_point *startpt = &(pts->pts[i]);
1391 pen_point *endpt = &(pts->pts[i+1]);
1392 int dx = endpt->x - startpt->x;
1393 int dy = endpt->y - startpt->y;
1394 int dircode = lialg_dctbl[dx+1][dy+1];
1396 /* assert(dircode < 8);*/
1397 startpt->chaincode = dircode;
1402 static region_list *lialg_compute_regions(point_list *pts) {
1403 region_list *regions = NULL;
1404 region_list *curr_reg = NULL;
1405 int *R[2 + LP_FILTER_ITERS];
1410 /* Initialize low-pass filter parameters if necessary. */
1411 if (lialg_lpfconst == -1)
1412 lialg_compute_lpf_parameters();
1414 /* Allocate a 2 x pts->npts array for use in computing the (filtered) Angle set, A_n. */
1415 /* junk = (int *)safe_malloc((2 + LP_FILTER_ITERS) * pts->npts * sizeof(int)); */
1416 junk = allocate((2 + LP_FILTER_ITERS) * pts->npts, int);
1417 for (i = 0; i < (2 + LP_FILTER_ITERS); i++)
1418 R[i] = junk + (i * pts->npts);
1421 /* Compute the Angle set, A, in the first element of array R. */
1422 /* Values in R are in degrees, x 100. */
1423 curr[0] = 18000; /* a_0 */
1424 for (i = 1; i < (pts->npts - 1); i++) {
1425 int d_i = pts->pts[i].chaincode;
1426 int d_iminusone = pts->pts[i-1].chaincode;
1429 if (d_iminusone < d_i)
1432 a_i = (d_iminusone - d_i) % 8;
1434 /* convert to degrees, x 100 */
1435 curr[i] = ((12 - a_i) % 8) * 45 * 100;
1437 curr[pts->npts - 1] = 18000; /* a_L-1 */
1439 /* Perform a number of filtering iterations. */
1441 for (j = 0; j < LP_FILTER_ITERS; j++, curr = R[j], next = R[j+1]) {
1442 for (i = 0; i < pts->npts; i++) {
1447 for (k = i - LP_FILTER_WIDTH; k <= i + LP_FILTER_WIDTH; k++) {
1448 int oldval = (k < 0 || k >= pts->npts) ? 18000 : curr[k];
1449 next[i] += oldval * lialg_lpfwts[k - (i - LP_FILTER_WIDTH)]; /* overflow? */
1452 next[i] /= lialg_lpfconst;
1456 /* Do final thresholding around PI. */
1457 /* curr and next are set-up correctly at end of previous loop! */
1458 for (i = 0; i < pts->npts; i++) {
1459 next[i] = (abs(curr[i] - 18000) < LP_FILTER_THLD)
1467 for (i = 0; i < pts->npts; i++) {
1468 fprintf(stderr, "%3d: (%3d, %3d) %ld ",
1469 i, pts->pts[i].x, pts->pts[i].y, pts->pts[i].chaincode);
1470 for (j = 0; j < 2 + LP_FILTER_ITERS; j++)
1471 fprintf(stderr, "%d ", R[j][i]);
1472 fprintf(stderr, "\n");
1476 /* Do the region segmentation. */
1481 #define RGN_TYPE(val)\
1484 : ((val) < 18000 ? RGN_CONCAVE : RGN_CONVEX))
1487 currtype = RGN_TYPE(curr[0]);
1488 /* regions = (region_list *)safe_malloc(sizeof(region_list));*/
1489 regions = allocate(1, region_list);
1491 curr_reg->start = start;
1493 curr_reg->type = currtype;
1494 curr_reg->next = NULL;
1495 for (i = 1; i < pts->npts; i++) {
1496 int nexttype = RGN_TYPE(curr[i]);
1498 if (nexttype != currtype) {
1499 region_list *next_reg = NULL;
1502 curr_reg->end = end;
1504 fprintf(stderr, " (%d, %d), %d\n", start, end, currtype);
1507 currtype = nexttype;
1508 /* next_reg = (region_list *)safe_malloc(sizeof(region_list));*/
1509 next_reg = allocate(1, region_list);
1510 next_reg->start = start;
1512 next_reg->type = nexttype;
1513 next_reg->next = NULL;
1515 curr_reg->next = next_reg;
1516 curr_reg = next_reg;
1520 curr_reg->end = end;
1522 fprintf(stderr, " (%d, %d), %d\n", start, end, currtype);
1524 /* Filter out convex/concave regions that are too short. */
1525 for (curr_reg = regions; curr_reg; curr_reg = curr_reg->next)
1526 if (curr_reg->type == RGN_PLAIN) {
1527 region_list *next_reg;
1529 for (next_reg = curr_reg->next;
1531 (next_reg->end - next_reg->start) < LP_FILTER_MIN;
1532 next_reg = curr_reg->next) {
1533 /* next_reg must not be plain, and it must be followed by a plain */
1534 /* assert(next_reg->type != RGN_PLAIN); */
1535 /* assert(next_reg->next != NULL && (next_reg->next)->type == RGN_PLAIN); */
1537 curr_reg->next = (next_reg->next)->next;
1538 curr_reg->end = (next_reg->next)->end;
1540 free(next_reg->next);
1545 /* Add-in pseudo-extremes. */
1547 region_list *tmp, *prev_reg;
1552 for (curr_reg = tmp; curr_reg; curr_reg = curr_reg->next) {
1553 if (curr_reg->type == RGN_PLAIN) {
1554 int arclen = lialg_compute_pathlen_subset(pts,
1557 int dx = pts->pts[curr_reg->end].x -
1558 pts->pts[curr_reg->start].x;
1559 int dy = pts->pts[curr_reg->end].y -
1560 pts->pts[curr_reg->start].y;
1561 int chordlen = isqrt(10000 * (dx * dx + dy * dy));
1562 int atcr = (chordlen == 0) ? 0 : (100 * arclen + chordlen / 2) / chordlen;
1565 fprintf(stderr, "%d, %d, %d\n", arclen, chordlen, atcr);
1567 /* Split region if necessary. */
1568 if (arclen >= PE_AL_THLD && atcr >= PE_ATCR_THLD) {
1569 int mid = curr_reg->start + (curr_reg->end - curr_reg->start) / 2;
1570 int end = curr_reg->end;
1571 region_list *saved_next = curr_reg->next;
1573 curr_reg->end = mid - 1;
1574 if (prev_reg == NULL)
1577 prev_reg->next = curr_reg;
1578 prev_reg = curr_reg;
1580 /* curr_reg = (region_list *)safe_malloc(sizeof(region_list));*/
1581 curr_reg = allocate(1, region_list);
1582 curr_reg->start = mid;
1583 curr_reg->end = mid;
1584 curr_reg->type = RGN_PSEUDO;
1585 curr_reg->next = NULL;
1586 prev_reg->next = curr_reg;
1587 prev_reg = curr_reg;
1589 /* curr_reg = (region_list *)safe_malloc(sizeof(region_list)); */
1590 curr_reg = allocate(1, region_list);
1591 curr_reg->start = mid + 1;
1592 curr_reg->end = end;
1593 curr_reg->type = RGN_PLAIN;
1594 curr_reg->next = NULL;
1595 prev_reg->next = curr_reg;
1596 prev_reg = curr_reg;
1598 curr_reg->next = saved_next;
1603 if (prev_reg == NULL)
1606 prev_reg->next = curr_reg;
1607 prev_reg = curr_reg;
1617 static point_list *lialg_compute_dompts(point_list *pts, region_list *regions) {
1618 point_list *dpts = NULL;
1624 /* Compute contour angle set. */
1625 cas = lialg_compute_contour_angle_set(pts, regions);
1626 /* assert(cas != NULL);*/
1628 /* Dominant points include: start_pt, end_pt, extrema_of_non_plain_regions, midpts of the preceding. */
1630 for (r = regions; r != NULL; r = r->next)
1631 if (r->type != RGN_PLAIN) nonplain++;
1632 ndpts = 2 * (2 + nonplain) - 1;
1633 /* dpts = (point_list *)safe_malloc(sizeof(point_list)); */
1634 dpts = allocate(1, point_list);
1635 dpts->pts = make_pen_point_array(ndpts);
1636 if (dpts->pts == NULL) {
1643 /* Pick out dominant points. */
1650 /* Record start point. */
1653 dpts->pts[dp++] = pts->pts[previx];
1655 for (curr = regions; curr != NULL; curr = curr->next)
1656 if (curr->type != RGN_PLAIN) {
1663 for (i = curr->start; i <= curr->end; i++) {
1665 if (v > max_v) { max_v = v; max_ix = i; }
1666 if (v < min_v) { min_v = v; min_ix = i; }
1668 fprintf(stderr, " %d\n", v);
1671 currix = (curr->type == RGN_CONVEX ? max_ix : min_ix);
1673 /* Record midpoint. */
1674 dpts->pts[dp++] = pts->pts[previx + (currix - previx) / 2];
1676 /* Record extreme point. */
1677 dpts->pts[dp++] = pts->pts[currix];
1682 /* Record last mid-point and end point. */
1683 currix = pts->npts - 1;
1684 dpts->pts[dp++] = pts->pts[previx + (currix - previx) / 2];
1685 dpts->pts[dp++] = pts->pts[currix];
1688 /* Compute chain-code. */
1689 lialg_compute_chain_code(dpts);
1696 static int *lialg_compute_contour_angle_set(point_list *pts,
1697 region_list *regions) {
1699 region_list *curr_reg, *prev_reg;
1702 /* V = (int *)safe_malloc(pts->npts * sizeof(int));*/
1703 V = allocate(pts->npts, int);
1706 for (curr_reg = regions; curr_reg != NULL;
1707 prev_reg = curr_reg, curr_reg = curr_reg->next) {
1708 for (i = curr_reg->start; i <= curr_reg->end; i++) {
1709 if (curr_reg->type == RGN_PLAIN) {
1714 /* XXX - eliminate floating point */
1715 region_list *next_reg = curr_reg->next;
1716 int b = curr_reg->start;
1717 int h = prev_reg->start;
1718 int t = next_reg->end;
1719 int pts_before = i - h;
1720 int pts_after = t - i;
1721 int min_pts = (pts_before < pts_after)
1724 int k = (min_pts < T_ONE)
1731 for (j = 1; j <= k; j++) {
1733 int ptB = i + j - 1;
1734 int d_A = pts->pts[ptA].chaincode;
1735 int d_B = pts->pts[ptB].chaincode;
1741 a_i = (d_A - d_B) % 8;
1743 /* convert to radians */
1744 if (a_i == 4 && curr_reg->type == RGN_CONVEX)
1747 sum += (float)((12 - a_i) % 8) / 4.0 * M_PI;
1749 V[i] = sum / (float)k;
1751 /* For now, simply choose the mid-point. */
1752 int isMidPt = (i == (curr_reg->start +
1753 (curr_reg->end - curr_reg->start) / 2));
1754 V[i] = (curr_reg->type == RGN_CONVEX)
1755 ? (isMidPt ? 18000 : 0)
1756 : (isMidPt ? 0 : 18000);
1761 V[pts->npts - 1] = 18000;
1768 * First compute the similarity between the two strings.
1769 * If it's above a threshold, compute the distance between
1770 * the two and return it as the ``score.''
1771 * Otherwise, return the constant WORST_SCORE.
1774 static void lialg_score_stroke(point_list *input_dompts, point_list *curr_dompts, int *sim, int *dist) {
1778 *sim = lialg_compute_similarity(input_dompts, curr_dompts);
1779 if (*sim < SIM_THLD) goto done;
1781 *dist = lialg_compute_distance(input_dompts, curr_dompts);
1785 fprintf(stderr, "%d, %d\n", *sim, *dist);
1789 static int lialg_compute_similarity(point_list *input_dompts,
1790 point_list *curr_dompts) {
1798 /* A is the longer sequence, length N. */
1799 /* B is the shorter sequence, length M. */
1800 if (input_dompts->npts >= curr_dompts->npts) {
1802 N = input_dompts->npts;
1804 M = curr_dompts->npts;
1808 N = curr_dompts->npts;
1810 M = input_dompts->npts;
1813 /* Allocate and initialize the Gain matrix, G. */
1814 /* The size of G is M x (N + 1). */
1815 /* Note that row 0 is unused. */
1816 /* Similarities are x 10. */
1818 /* G = (int **)safe_malloc(M * sizeof(int *));*/
1819 G = allocate(M, int *);
1820 /* junk = (int *)safe_malloc(M * (N + 1) * sizeof(int)); */
1821 junk = allocate(M * (N + 1), int);
1822 for (i = 0; i < M; i++)
1823 G[i] = junk + (i * (N + 1));
1825 for (i = 1; i < M; i++) {
1826 int bval = B->pts[i-1].chaincode;
1828 /* Source column. */
1831 for (j = 1; j < N; j++) {
1832 int aval = A->pts[j-1].chaincode;
1833 int diff = abs(bval - aval);
1834 if (diff > 4) diff = 8 - diff;
1836 G[i][j] = (diff == 0)
1848 /* Do the DP algorithm. */
1849 /* Proceed in column order, from highest column to the lowest. */
1850 /* Within each column, proceed from the highest row to the lowest. */
1851 /* Skip the highest column. */
1853 for (j = N - 1; j >= 0; j--)
1854 for (i = M - 1; i > 0; i--) {
1855 int max = G[i][j + 1];
1858 int tmp = G[i + 1][j + 1];
1859 if (tmp > max) max = tmp;
1865 sim = (10 * G[1][0] + (N - 1) / 2) / (N - 1);
1868 if (G != NULL) free(G);
1869 if (junk != NULL) free(junk);
1874 static int lialg_compute_distance(point_list *input_dompts,
1875 point_list *curr_dompts) {
1876 int dist = MAX_DIST;
1885 /* A is the longer sequence, length N. */
1886 /* B is the shorter sequence, length M. */
1887 if (input_dompts->npts >= curr_dompts->npts) {
1889 N = input_dompts->npts;
1891 M = curr_dompts->npts;
1895 N = curr_dompts->npts;
1897 M = input_dompts->npts;
1900 /* Construct the helper vectors, BE and TE, which say for each column */
1901 /* what are the ``bottom'' and ``top'' rows of interest. */
1903 /* BE = (int *)safe_malloc((N + 1) * sizeof(int));*/
1904 BE = allocate((N + 1), int);
1905 /* TE = (int *)safe_malloc((N + 1) * sizeof(int)); */
1906 TE = allocate((N + 1), int);
1908 for (j = 1; j <= N; j++) {
1911 bot = j + (M - DP_BAND);
1912 if (bot > M) bot = M;
1915 top = j - (N - DP_BAND);
1916 if (top < 1) top = 1;
1921 /* Allocate and initialize the Cost matrix, C. */
1922 /* The size of C is (M + 1) x (N + 1). */
1923 /* Note that row and column 0 are unused. */
1924 /* Costs are x 100. */
1926 /* C = (int **)safe_malloc((M + 1) * sizeof(int *)); */
1927 C = allocate((M + 1), int *);
1928 /* junk = (int *)safe_malloc((M + 1) * (N + 1) * sizeof(int)); */
1929 junk = allocate((M + 1) * (N + 1), int);
1930 for (i = 0; i <= M; i++)
1931 C[i] = junk + (i * (N + 1));
1933 for (i = 1; i <= M; i++) {
1934 int bx = B->pts[i-1].x;
1935 int by = B->pts[i-1].y;
1937 for (j = 1; j <= N; j++) {
1938 int ax = A->pts[j-1].x;
1939 int ay = A->pts[j-1].y;
1942 int dist = isqrt(10000 * (dx * dx + dy * dy));
1949 /* Do the DP algorithm. */
1950 /* Proceed in column order, from highest column to the lowest. */
1951 /* Within each column, proceed from the highest row to the lowest. */
1953 for (j = N; j > 0; j--)
1954 for (i = M; i > 0; i--) {
1957 if (i > BE[j] || i < TE[j] || (j == N && i == M))
1962 int tmp = C[i][j+1];
1963 if (tmp < min) min = tmp;
1967 int tmp = C[i+1][j+1];
1968 if (tmp < min) min = tmp;
1973 int tmp = C[i+1][j];
1974 if (tmp < min) min = tmp;
1980 dist = (C[1][1] + N / 2) / N;
1983 if (C != NULL) free(C);
1984 if (junk != NULL) free(junk);
1985 if (BE != NULL) free(BE);
1986 if (TE != NULL) free(TE);
1991 /*************************************************************
1993 Digest-processing routines
1995 *************************************************************/
1997 static int lialg_read_classifier_digest(rClassifier *rec) {
2001 /* Try to open the corresponding digest file. */
2006 /* Get a copy of the filename, with some room on the end. */
2007 /* clx_path = safe_malloc(strlen(rec->file_name) + 5); */
2008 clx_path = allocate(strlen(rec->file_name) + 5, char);
2009 strcpy(clx_path, rec->file_name);
2011 /* Truncate the path after the last dot. */
2012 dot = strrchr(clx_path, '.');
2013 if (dot == NULL) { free(clx_path); return(-1); }
2016 /* Append the classifier-digest extension. */
2017 strcat(clx_path, "clx");
2019 fp = fopen(clx_path, "r");
2020 if (fp == NULL) { free(clx_path); return(-1); }
2025 /* Read-in the name and dominant points for each class. */
2026 for (nclasses = 0; !feof(fp); nclasses++) {
2027 point_list *dpts = NULL;
2032 if (fscanf(fp, "%s %d", class, &npts) != 2) {
2033 if (feof(fp)) break;
2037 rec->cnames[nclasses] = strdup(class);
2039 /* Allocate a dominant-points list. */
2040 /* dpts = (point_list *)safe_malloc(sizeof(point_list)); */
2041 dpts = allocate(1, point_list);
2042 dpts->pts = make_pen_point_array(npts);
2043 if (dpts->pts == NULL) goto failed;
2047 /* Read in each point. */
2048 for (j = 0; j < npts; j++) {
2051 if (fscanf(fp, "%d %d", &x, &y) != 2) goto failed;
2056 /* Compute the chain-code. */
2057 lialg_compute_chain_code(dpts);
2059 /* Store the list in the rec data structure. */
2060 rec->dompts[nclasses] = dpts;
2065 fprintf(stderr, "read_classifier_digest failed...\n");
2066 for (; nclasses >= 0; nclasses--) {
2067 if (rec->cnames[nclasses] != NULL) {
2068 free(rec->cnames[nclasses]);
2069 rec->cnames[nclasses] = NULL;
2071 if (rec->dompts[nclasses] != NULL) {
2072 delete_examples(rec->dompts[nclasses]);
2073 rec->dompts[nclasses] = NULL;
2077 delete_examples(dpts);
2087 /*************************************************************
2089 Canonicalization routines
2091 *************************************************************/
2093 static int lialg_canonicalize_examples(rClassifier *rec) {
2098 fprintf(stderr, "lialg_canonicalize_examples working on %s\n",
2101 /* Initialize canonical-example arrays. */
2102 for (i = 0; i < MAXSCLASSES; i++) {
2103 rec->canonex[i] = NULL;
2106 /* Figure out number of classes. */
2108 nclasses < MAXSCLASSES && rec->cnames[nclasses] != NULL;
2112 /* Canonicalize the examples for each class. */
2113 for (i = 0; i < nclasses; i++) {
2116 point_list *pts, *tmp, *avg;
2117 int maxxrange, maxyrange;
2118 int minx, miny, maxx, maxy;
2119 int avgxrange, avgyrange, avgxoff, avgyoff, avgscale;
2123 fprintf(stderr, "lialg_canonicalize_examples working on class %s\n",
2126 /* Make a copy of the examples. */
2129 for (nex = 0; tmp != NULL; nex++, tmp = tmp->next) {
2130 if ((pts = add_example(pts, tmp->npts, tmp->pts)) == NULL) {
2131 delete_examples(pts);
2136 /* Canonicalize each example, and derive the max x and y ranges. */
2139 for (j = 0, tmp = pts; j < nex; j++, tmp = tmp->next) {
2140 if (lialg_canonicalize_example_stroke(tmp) != 0) {
2142 fprintf(stderr, "lialg_canonicalize_example_stroke returned error\n");
2147 if (tmp->xrange > maxxrange) maxxrange = tmp->xrange;
2148 if (tmp->yrange > maxyrange) maxyrange = tmp->yrange;
2151 /* Normalize max ranges. */
2152 if (((100 * maxxrange + CANONICAL_X / 2) / CANONICAL_X) >
2153 ((100 * maxyrange + CANONICAL_Y / 2) / CANONICAL_Y)) {
2154 maxyrange = (maxyrange * CANONICAL_X + maxxrange / 2) / maxxrange;
2155 maxxrange = CANONICAL_X;
2158 maxxrange = (maxxrange * CANONICAL_Y + maxyrange / 2) / maxyrange;
2159 maxyrange = CANONICAL_Y;
2162 /* Re-scale each example to max ranges. */
2163 for (j = 0, tmp = pts; j < nex; j++, tmp = tmp->next) {
2164 int scalex = (tmp->xrange == 0) ? 100 : (100 * maxxrange + tmp->xrange / 2) / tmp->xrange;
2165 int scaley = (tmp->yrange == 0) ? 100 : (100 * maxyrange + tmp->yrange / 2) / tmp->yrange;
2166 if (lialg_translate_points(tmp, 0, 0, scalex, scaley) != 0) {
2167 delete_examples(pts);
2172 /* Average the examples; leave average in first example. */
2173 avg = pts; /* careful aliasing!! */
2174 for (k = 0; k < NCANONICAL; k++) {
2178 for (j = 0, tmp = pts; j < nex; j++, tmp = tmp->next) {
2179 xsum += tmp->pts[k].x;
2180 ysum += tmp->pts[k].y;
2183 avg->pts[k].x = (xsum + j / 2) / j;
2184 avg->pts[k].y = (ysum + j / 2) / j;
2187 /* Compute BB of averaged stroke and re-scale. */
2188 lialg_get_bounding_box(avg, &minx, &miny, &maxx, &maxy);
2189 avgxrange = maxx - minx;
2190 avgyrange = maxy - miny;
2191 avgscale = (((100 * avgxrange + CANONICAL_X / 2) / CANONICAL_X) >
2192 ((100 * avgyrange + CANONICAL_Y / 2) / CANONICAL_Y))
2193 ? (100 * CANONICAL_X + avgxrange / 2) / avgxrange
2194 : (100 * CANONICAL_Y + avgyrange / 2) / avgyrange;
2195 if (lialg_translate_points(avg, minx, miny, avgscale, avgscale) != 0) {
2196 delete_examples(pts);
2200 /* Re-compute the x and y ranges and center the stroke. */
2201 lialg_get_bounding_box(avg, &minx, &miny, &maxx, &maxy);
2202 avgxrange = maxx - minx;
2203 avgyrange = maxy - miny;
2204 avgxoff = -((CANONICAL_X - avgxrange + 1) / 2);
2205 avgyoff = -((CANONICAL_Y - avgyrange + 1) / 2);
2206 if (lialg_translate_points(avg, avgxoff, avgyoff, 100, 100) != 0) {
2207 delete_examples(pts);
2211 /* Create a point list to serve as the ``canonical representation. */
2212 if ((rec->canonex[i] = add_example(NULL, avg->npts, avg->pts)) == NULL) {
2213 delete_examples(pts);
2216 (rec->canonex[i])->xrange = maxx - minx;
2217 (rec->canonex[i])->yrange = maxy - miny;
2220 fprintf(stderr, "%s, avgpts = %d\n", rec->cnames[i], avg->npts);
2221 for (j = 0; j < avg->npts; j++) {
2222 fprintf(stderr, " (%d, %d)\n",
2223 avg->pts[j].x, avg->pts[j].y);
2227 /* Compute dominant points of canonical representation. */
2228 rec->dompts[i] = lialg_compute_dominant_points(avg);
2231 delete_examples(pts);
2235 for (i = 0; i < nclasses; i++) {
2236 char *best_name = lialg_recognize_stroke(rec, rec->canonex[i]);
2238 if (best_name != rec->cnames[i])
2239 fprintf(stderr, "%s, best = %s\n", rec->cnames[i], best_name);
2246 static int lialg_canonicalize_example_stroke(point_list *points) {
2247 int minx, miny, maxx, maxy, xrange, yrange, scale;
2249 /* Filter out points that are too close. */
2250 if (lialg_filter_points(points) != 0) return(-1);
2252 /* Must be at least two points! */
2253 if (points->npts < 2) {
2255 fprintf(stderr, "lialg_canonicalize_example_stroke: npts=%d\n",
2261 /* Scale up to avoid conversion errors. */
2262 lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
2263 xrange = maxx - minx;
2264 yrange = maxy - miny;
2265 scale = (((100 * xrange + CANONICAL_X / 2) / CANONICAL_X) >
2266 ((100 * yrange + CANONICAL_Y / 2) / CANONICAL_Y))
2267 ? (100 * CANONICAL_X + xrange / 2) / xrange
2268 : (100 * CANONICAL_Y + yrange / 2) / yrange;
2269 if (lialg_translate_points(points, minx, miny, scale, scale) != 0) {
2271 fprintf(stderr, "lialg_translate_points (minx=%d,miny=%d,scale=%d) returned error\n", minx, miny, scale);
2276 /* Compute an equivalent stroke with equi-distant points. */
2277 if (lialg_compute_equipoints(points) != 0) return(-1);
2279 /* Re-translate the points to the origin. */
2280 lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
2281 if (lialg_translate_points(points, minx, miny, 100, 100) != 0) {
2283 fprintf(stderr, "lialg_translate_points (minx=%d,miny=%d) returned error\n", minx, miny);
2288 /* Store the x and y ranges in the point list. */
2289 xrange = maxx - minx;
2290 yrange = maxy - miny;
2291 points->xrange = xrange;
2292 points->yrange = yrange;
2296 fprintf(stderr, "Canonicalized: %d, %d, %d, %d\n", minx, miny, maxx, maxy);
2297 for (i = 0; i < points->npts; i++)
2298 fprintf(stderr, " (%d %d)\n",
2299 points->pts[i].x, points->pts[i].y);
2307 static int lialg_compute_equipoints(point_list *points) {
2308 pen_point *equipoints = make_pen_point_array(NCANONICAL);
2309 int nequipoints = 0;
2310 int pathlen = lialg_compute_pathlen(points);
2311 int equidist = (pathlen + (NCANONICAL - 1) / 2) / (NCANONICAL - 1);
2313 int dist_since_last_eqpt;
2314 int remaining_seglen;
2315 int dist_to_next_eqpt;
2317 if (equipoints == NULL) {
2318 error("can't allocate memory in lialg_compute_equipoints");
2323 fprintf(stderr, "compute_equipoints: npts = %d, pathlen = %d, equidist = %d\n",
2324 points->npts, pathlen, equidist);
2328 /* First original point is an equipoint. */
2329 equipoints[0] = points->pts[0];
2331 dist_since_last_eqpt = 0;
2333 for (i = 1; i < points->npts; i++) {
2334 int dx1 = points->pts[i].x - points->pts[i-1].x;
2335 int dy1 = points->pts[i].y - points->pts[i-1].y;
2336 int endx = 100 * points->pts[i-1].x;
2337 int endy = 100 * points->pts[i-1].y;
2338 remaining_seglen = isqrt(10000 * (dx1 * dx1 + dy1 * dy1));
2339 dist_to_next_eqpt = equidist - dist_since_last_eqpt;
2341 while (remaining_seglen >= dist_to_next_eqpt) {
2343 /* x-coordinate stays the same */
2345 endy += dist_to_next_eqpt;
2347 endy -= dist_to_next_eqpt;
2350 int slope = (100 * dy1 + dx1 / 2) / dx1;
2351 int tmp = isqrt(10000 + slope * slope);
2352 int dx = (100 * dist_to_next_eqpt + tmp / 2) / tmp;
2353 int dy = (slope * dx + 50) / 100;
2355 if (dy < 0) dy = -dy;
2366 equipoints[nequipoints].x = (endx + 50) / 100;
2367 equipoints[nequipoints].y = (endy + 50) / 100;
2369 /* assert(nequipoints <= NCANONICAL);*/
2370 dist_since_last_eqpt = 0;
2371 remaining_seglen -= dist_to_next_eqpt;
2372 dist_to_next_eqpt = equidist;
2375 dist_since_last_eqpt += remaining_seglen;
2378 /* Take care of last equipoint. */
2379 if (nequipoints == NCANONICAL) {
2381 } else if (nequipoints == (NCANONICAL - 1)) {
2382 /* Make last original point the last equipoint. */
2383 equipoints[nequipoints] = points->pts[points->npts - 1];
2386 fprintf(stderr,"lialg_compute_equipoints: nequipoints = %d\n",
2393 points->npts = NCANONICAL;
2394 delete_pen_point_array(points->pts);
2395 points->pts = equipoints;
2400 /*************************************************************
2404 *************************************************************/
2406 /* Result is x 100. */
2407 static int lialg_compute_pathlen(point_list *points) {
2408 return(lialg_compute_pathlen_subset(points, 0, points->npts - 1));
2412 /* Result is x 100. */
2413 static int lialg_compute_pathlen_subset(point_list *points,
2414 int start, int end) {
2419 for (i = start + 1; i <= end; i++) {
2420 int dx = points->pts[i].x - points->pts[i-1].x;
2421 int dy = points->pts[i].y - points->pts[i-1].y;
2422 int dist = isqrt(10000 * (dx * dx + dy * dy));
2430 /* Note that this does NOT update points->xrange and points->yrange! */
2431 static int lialg_filter_points(point_list *points) {
2433 pen_point *filtered_pts = make_pen_point_array(points->npts);
2436 if (filtered_pts == NULL) {
2437 error("can't allocate memory in lialg_filter_points");
2441 filtered_pts[0] = points->pts[0];
2443 for (i = 1; i < points->npts; i++) {
2444 int j = filtered_npts - 1;
2445 int dx = points->pts[i].x - filtered_pts[j].x;
2446 int dy = points->pts[i].y - filtered_pts[j].y;
2447 int magsq = dx * dx + dy * dy;
2449 if (magsq >= DIST_SQ_THRESHOLD) {
2450 filtered_pts[filtered_npts] = points->pts[i];
2455 points->npts = filtered_npts;
2456 delete_pen_point_array(points->pts);
2457 points->pts = filtered_pts;
2462 /* scalex and scaley are x 100. */
2463 /* Note that this does NOT update points->xrange and points->yrange! */
2464 static int lialg_translate_points(point_list *points,
2466 int scalex, int scaley) {
2469 for (i = 0; i < points->npts; i++) {
2470 points->pts[i].x = ((points->pts[i].x - minx) * scalex + 50) / 100;
2471 points->pts[i].y = ((points->pts[i].y - miny) * scaley + 50) / 100;
2478 static void lialg_get_bounding_box(point_list *points,
2479 int *pminx, int *pminy,
2480 int *pmaxx, int *pmaxy) {
2481 int minx, miny, maxx, maxy;
2484 minx = maxx = points->pts[0].x;
2485 miny = maxy = points->pts[0].y;
2486 for (i = 1; i < points->npts; i++) {
2487 pen_point *pt = &(points->pts[i]);
2488 if (pt->x < minx) minx = pt->x;
2489 if (pt->x > maxx) maxx = pt->x;
2490 if (pt->y < miny) miny = pt->y;
2491 if (pt->y > maxy) maxy = pt->y;
2504 return exp((double)x);
2509 static void lialg_compute_lpf_parameters() {
2512 for (i = LP_FILTER_WIDTH; i >= 0; i--) {
2513 float x = 0.04 * (i * i);
2514 #if defined(ARM_LINUX) || !defined(__GLIBC__)
2515 double tmp = 100.0 * exp((double)x);
2517 float tmp = 100.0 * expf(x);
2519 int wt = rint((double)tmp);
2521 lialg_lpfwts[LP_FILTER_WIDTH - i] = wt;
2522 lialg_lpfwts[LP_FILTER_WIDTH + i] = wt;
2525 for (i = 0; i < (2 * LP_FILTER_WIDTH + 1); i++) {
2526 lialg_lpfconst += lialg_lpfwts[i];
2531 /* Code from Joseph Hall (jnhall@sat.mot.com). */
2532 static int isqrt(int n) {
2534 register long k0, k1, nn;
2536 for (nn = i = n, k0 = 2; i > 0; i >>= 2, k0 <<= 1)
2540 k1 = (nn / k0 + k0) >> 1;
2541 if (((k0 ^ k1) & ~1) == 0)
2545 return (int) ((k1 + 1) >> 1);
2549 /* Helper routines from Mark Hayter. */
2550 static int likeatan(int tantop, int tanbot) {
2552 /* Use tan(theta)=top/bot --> order for t */
2553 /* t in range 0..0x40000 */
2555 if ((tantop == 0) && (tanbot == 0))
2559 t = (tantop << 16) / (abs(tantop) + abs(tanbot));
2563 if (tantop < 0) t = 0x40000 + t;
2569 static int quadr(int t) {
2570 return (8 - (((t + 0x4000) >> 15) & 7)) & 7;