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.
24 #include "scribbleimpl.h"
26 #include "hre_internal.h"
27 #include "li_recognizer_internal.h"
31 #define LI_MAGIC 0xACCBADDD
33 #define CHECK_LI_MAGIC(_a) \
34 ((_a) != nil && ((li_recognizer*)(_a))->li_magic == LI_MAGIC)
37 static void lialg_initialize(rClassifier *);
38 static int lialg_read_classifier_digest(rClassifier *);
39 static int lialg_canonicalize_examples(rClassifier *);
40 static char *lialg_recognize_stroke(rClassifier *, point_list *);
41 static void lialg_compute_lpf_parameters(void);
44 char* li_err_msg = nil;
46 #define bcopy(s1,s2,n) memmove(s2,s1,n)
48 /*Freeing classifier*/
51 free_rClassifier(rClassifier* rc);
58 add_example(point_list* l,int npts,pen_point* pts)
60 pen_point* lpts = mallocz(npts*sizeof(pen_point), 1);
61 point_list *p = malloc(sizeof(point_list));
65 p->next = l; /*Order doesn't matter, so we stick on end.*/
69 bcopy(pts, lpts, npts * sizeof(pen_point));
76 delete_examples(point_list* l)
80 for( ; l != nil; l = p ) {
92 * recognize_internal-Form Vector, use Classifier to classify, return char.
96 recognize_internal(rClassifier* rec, Stroke* str, int*)
101 stroke = add_example(nil, str->npts, str->pts);
102 if (stroke == nil) return(nil);
104 res = lialg_recognize_stroke(rec, stroke);
106 delete_examples(stroke);
111 * file_path-Construct pathname, check for proper extension.
115 file_path(char* dir,char* filename,char* pathname)
119 /*Check for proper extension on file name.*/
121 dot = strrchr(filename,'.');
127 /*Determine whether a gesture or character classifier.*/
129 if( strcmp(dot,LI_CLASSIFIER_EXTENSION) != 0 ) {
133 /*Concatenate directory and filename into pathname.*/
135 strcpy(pathname,dir);
136 strcat(pathname,"/");
137 strcat(pathname,filename);
142 /*read_classifier_points-Read points so classifier can be extended.*/
145 read_classifier_points(FILE* fd,int nclss,point_list** ex,char** cnames)
150 char* names[MAXSCLASSES];
151 point_list* examples[MAXSCLASSES];
157 for( i = 0; i < MAXSCLASSES; i++ ) {
164 for( k = 0; k < nclss; k++ ) {
166 /*Read class name and number of examples.*/
168 if( fscanf(fd,"%d %s",&nex,buf) != 2 )
173 names[k] = strdup(buf);
177 for( i = 0; i < nex; i++ ) {
179 /*Read number of points.*/
181 if( fscanf(fd,"%d",&npts) != 1 )
182 goto unallocate; /*Boy would I like exceptions!*/
184 /*Allocate array for points.*/
186 if( (pts = mallocz(npts*sizeof(pen_point), 1)) == nil )
191 for( j = 0; j < npts; j++ ) {
193 if( fscanf(fd,"%d %d",&x,&y) != 2 ) {
194 delete_pen_point_array(pts);
197 pts[j].Point = Pt(x, y);
202 if( (examples[k] = add_example(examples[k],npts,pts)) == nil ) {
203 delete_pen_point_array(pts);
207 delete_pen_point_array(pts);
212 /* ari -- end of list of classes */
213 /* fprint(2, "]\n"); */
215 /*Transfer to recognizer.*/
217 bcopy(examples,ex,sizeof(examples));
218 bcopy(names,cnames,sizeof(names));
222 /*Error. Deallocate memory and return.*/
225 for( ; k >= 0; k-- ) {
226 delete_examples(examples[k]);
233 /*read_classifier-Read a classifier file.*/
235 static int read_classifier(FILE* fd,rClassifier* rc)
240 /*Read in classifier file.*/
242 if(fscanf(fd, "%d", &rc->nclasses) != 1)
245 /*Read in the example points, so classifier can be extended.*/
247 if( read_classifier_points(fd,rc->nclasses,rc->ex,rc->cnames) != 0 )
254 * Extension Functions
257 /* getClasses and clearState are by Ari */
260 recognizer_getClasses (recognizer r, char ***list, int *nc)
266 rec = (li_recognizer*)r->recognizer_specific;
268 /*Check for LI recognizer.*/
270 if( !CHECK_LI_MAGIC(rec) ) {
271 li_err_msg = "Not a LI recognizer";
275 *nc = nclasses = rec->li_rc.nclasses;
276 ret = malloc(nclasses*sizeof(char*));
278 for (i = 0; i < nclasses; i++)
279 ret[i] = rec->li_rc.cnames[i]; /* only the 1st char of the cname */
285 recognizer_clearState (recognizer)
287 /*This operation isn't supported by the LI recognizer.*/
289 li_err_msg = "Clearing state is not supported by the LI recognizer";
294 static bool isa_li(recognizer r)
295 { return(CHECK_LI_MAGIC(r)); }
298 recognizer_train(recognizer, rc*, uint, Stroke*, rec_element*, bool)
300 /*This operation isn't supported by the LI recognizer.*/
302 li_err_msg = "Training is not supported by the LI recognizer";
308 li_recognizer_get_example (recognizer r,
315 li_recognizer *rec = (li_recognizer*)r->recognizer_specific;
316 int nclasses = rec->li_rc.nclasses;
319 if( !CHECK_LI_MAGIC(rec) ) {
320 li_err_msg = "Not a LI recognizer";
323 if (class > nclasses)
325 pl = rec->li_rc.canonex[class];
326 while (instance && pl)
333 *name = rec->li_rc.cnames[class];
336 return pl->npts; /* I hope [sjm] */
344 /*li_recognizer_load-Load a classifier file.*/
346 static int li_recognizer_load(recognizer r, char* dir, char* filename)
353 rec = (li_recognizer*)r->recognizer_specific;
355 /*Make sure recognizer's OK*/
357 if( !CHECK_LI_MAGIC(rec) ) {
358 li_err_msg = "Not a LI recognizer";
364 /*Check parameters.*/
366 if( filename == nil ) {
367 li_err_msg = "Invalid parameters";
371 /*We let the directory be null.*/
373 if( dir == nil || (int)strlen(dir) <= 0 ) {
377 if(0)fprint(2, "dir = %s, filename = %s\n",
380 /*Make full pathname and check filename*/
382 pathname = malloc((strlen(dir) + strlen(filename) + 2)*sizeof(char));
383 if( file_path(dir, filename, pathname) == -1 ) {
385 li_err_msg = "Not a LI recognizer classifier file";
389 /* Try to short-circuit the full classifier-file processing. */
390 rc->file_name = pathname;
391 if (lialg_read_classifier_digest(rc) == 0)
397 if( (fd = fopen(pathname,"r")) == nil ) {
398 li_err_msg = "Can't open classifier file";
399 if(0)fprint(2, "Can't open %s.\n", pathname);
404 /*If rClassifier is OK, then delete it first.*/
406 if( rc->file_name != nil ) {
407 free_rClassifier(rc);
412 if( read_classifier(fd,rc) < 0 ) {
421 /*Add classifier name.*/
423 rc->file_name = pathname;
425 /* Canonicalize examples. */
426 if (lialg_canonicalize_examples(rc) != 0) {
435 /*li_recognizer_save-Save a classifier file.*/
437 static int li_recognizer_save(recognizer, char*, char*)
439 /*This operation isn't supported by the LI recognizer.*/
441 li_err_msg = "Saving is not supported by the LI recognizer";
446 li_recognizer_load_dictionary(recognizer, char*, char*)
448 /*This operation isn't supported by the LI recognizer.*/
450 li_err_msg = "Dictionaries are not supported by the LI recognizer";
455 li_recognizer_save_dictionary(recognizer, char*, char*, wordset)
457 /*This operation isn't supported by the LI recognizer.*/
458 li_err_msg = "Dictionaries are not supported by the LI recognizer";
463 li_recognizer_free_dictionary(recognizer, wordset)
465 /*This operation isn't supported by the LI recognizer.*/
467 li_err_msg = "Dictionaries are not supported by the LI recognizer";
474 li_recognizer_add_to_dictionary(recognizer, letterset*, wordset)
476 /*This operation isn't supported by the LI recognizer.*/
477 li_err_msg = "Dictionaries are not supported by the LI recognizer";
482 li_recognizer_delete_from_dictionary(recognizer, letterset*, wordset)
484 /*This operation isn't supported by the LI recognizer.*/
485 li_err_msg = "Dictionaries are not supported by the LI recognizer";
490 li_recognizer_error(recognizer rec)
492 char* ret = li_err_msg;
494 /*Check for LI recognizer.*/
496 if( !CHECK_LI_MAGIC(rec->recognizer_specific) ) {
497 li_err_msg = "Not a LI recognizer";
505 li_recognizer_clear(recognizer r, bool)
509 rec = (li_recognizer*)r->recognizer_specific;
510 /*Check for LI recognizer.*/
511 if( !CHECK_LI_MAGIC(rec) ) {
512 li_err_msg = "Not a LI recognizer";
519 li_recognizer_set_context(recognizer, rc*)
521 /*This operation isn't supported by the LI recognizer.*/
522 li_err_msg = "Contexts are not supported by the LI recognizer";
527 li_recognizer_get_context(recognizer)
529 /*This operation isn't supported by the LI recognizer.*/
530 li_err_msg = "Contexts are not supported by the LI recognizer";
535 li_recognizer_get_buffer(recognizer, uint*, Stroke**)
537 /*This operation isn't supported by the LI recognizer.*/
538 li_err_msg = "Buffer get/set are not supported by the LI recognizer";
543 li_recognizer_set_buffer(recognizer, uint, Stroke*)
545 /*This operation isn't supported by the LI recognizer.*/
546 li_err_msg = "Buffer get/set are not supported by the LI recognizer";
551 li_recognizer_translate(recognizer r, uint ncs, Stroke* tps, bool, int* nret, rec_alternative** ret)
558 rec = (li_recognizer*)r->recognizer_specific;
563 /*Check for LI recognizer.*/
565 if( !CHECK_LI_MAGIC(rec) ) {
566 li_err_msg = "Not a LI recognizer";
572 /*Check for valid parameters.*/
574 li_err_msg = "Invalid parameters: ncs";
578 li_err_msg = "Invalid parameters: tps";
582 li_err_msg = "Invalid parameters: nret";
586 li_err_msg = "Invalid parameters: ret";
591 * Go through the stroke array and recognize. Since this is a single
592 * stroke recognizer, each stroke is treated as a separate
593 * character or gesture. We allow only characters or gestures
594 * to be recognized at one time, since otherwise, handling
595 * the display of segmentation would be difficult.
597 clss = recognize_internal(rc,tps,&conf);
610 li_recognizer_get_extension_functions(recognizer rec)
614 /*Check for LI recognizer.*/
616 if( !CHECK_LI_MAGIC(rec->recognizer_specific) ) {
617 li_err_msg = "Not a LI recognizer";
621 ret = make_rec_fn_array(LI_NUM_EX_FNS);
623 /* ari -- clearState & getClasses are mine */
624 ret[LI_GET_CLASSES] = (rec_fn)recognizer_getClasses;
625 ret[LI_CLEAR] = (rec_fn)recognizer_clearState;
626 ret[LI_ISA_LI] = (rec_fn)isa_li;
627 ret[LI_TRAIN] = (rec_fn)recognizer_train;
633 li_recognizer_get_gesture_names(recognizer)
635 /*This operation isn't supported by the LI recognizer.*/
636 li_err_msg = "Gestures are not supported by the LI recognizer";
641 li_recognizer_set_gesture_action(recognizer, char*, xgesture, void*)
643 /*This operation isn't supported by the LI recognizer.*/
644 li_err_msg = "Gestures are not supported by the LI recognizer";
652 /*RECOGNIZER_INITIALIZE-Initialize the recognizer.*/
654 /* note from ari: this expands via pre-processor to
656 * recognizer __recognizer_internal_initialize(rec_info* ri)
659 RECOGNIZER_INITIALIZE(ri)
665 /*Check that locale matches.*/
667 if( strcmp(ri->ri_locale,LI_SUPPORTED_LOCALE) != 0 ) {
668 li_err_msg = "Not a supported locale";
669 /* fprint(2, "Locale error.\n");*/
674 * Check that character sets match. Note that this is only approximate,
675 * since the classifier file will have more information.
678 if( ri->ri_subset != nil ) {
679 for(i = 0; ri->ri_subset[i] != nil; i++ ) {
681 if( strcmp(ri->ri_subset[i],UPPERCASE) != 0 &&
682 strcmp(ri->ri_subset[i],LOWERCASE) != 0 &&
683 strcmp(ri->ri_subset[i],DIGITS) != 0 &&
684 strcmp(ri->ri_subset[i],GESTURE) != 0 ) {
685 li_err_msg = "Not a supported character set";
686 /* fprint(2, "charset error.\n"); */
694 r = make_recognizer(ri);
695 /* fprint(2, "past make_recognizer.\n"); */
698 li_err_msg = "Can't allocate storage";
703 /*Make a LI recognizer structure.*/
706 /* rec = (li_recognizer*)safe_malloc(sizeof(li_recognizer))) == nil ); */
708 rec = malloc(sizeof(li_recognizer));
710 r->recognizer_specific = rec;
712 rec->li_rc.file_name = nil;
713 rec->li_rc.nclasses = 0;
715 /*Initialize the recognizer struct.*/
717 r->recognizer_load_state = li_recognizer_load;
718 r->recognizer_save_state = li_recognizer_save;
719 r->recognizer_load_dictionary = li_recognizer_load_dictionary;
720 r->recognizer_save_dictionary = li_recognizer_save_dictionary;
721 r->recognizer_free_dictionary = li_recognizer_free_dictionary;
722 r->recognizer_add_to_dictionary = li_recognizer_add_to_dictionary;
723 r->recognizer_delete_from_dictionary = li_recognizer_delete_from_dictionary;
724 r->recognizer_error = li_recognizer_error;
725 r->recognizer_translate = li_recognizer_translate;
726 r->recognizer_get_context = li_recognizer_get_context;
727 r->recognizer_set_context = li_recognizer_set_context;
728 r->recognizer_get_buffer = li_recognizer_get_buffer;
729 r->recognizer_set_buffer = li_recognizer_set_buffer;
730 r->recognizer_clear = li_recognizer_clear;
731 r->recognizer_get_extension_functions =
732 li_recognizer_get_extension_functions;
733 r->recognizer_get_gesture_names = li_recognizer_get_gesture_names;
734 r->recognizer_set_gesture_action =
735 li_recognizer_set_gesture_action;
737 /*Initialize LI Magic Number.*/
739 rec->li_magic = LI_MAGIC;
741 /*Initialize rClassifier.*/
743 rec->li_rc.file_name = nil;
745 for( i = 0; i < MAXSCLASSES; i++ ) {
746 rec->li_rc.ex[i] = nil;
747 rec->li_rc.cnames[i] = nil;
750 lialg_initialize(&rec->li_rc);
752 /*Get rid of error message. Not needed here.*/
758 /*free_rClassifier-Free the rClassifier.*/
761 free_rClassifier(rClassifier* rc)
765 if( rc->file_name != nil) {
769 for( i = 0; rc->ex[i] != nil; i++) {
770 delete_examples(rc->ex[i]);
776 /*RECOGNIZER_FINALIZE-Deallocate the recognizer, finalize.*/
778 RECOGNIZER_FINALIZE(r)
780 li_recognizer* rec = (li_recognizer*)r->recognizer_specific;
782 /*Make sure this is a li_recognizer first*/
783 if( !CHECK_LI_MAGIC(rec) ) {
784 li_err_msg = "Not a LI recognizer";
788 /*Deallocate rClassifier resources.*/
790 free_rClassifier(&(rec->li_rc));
792 /*Deallocate the li_recognizer struct.*/
795 /*Deallocate the recognizer*/
796 delete_recognizer(r);
802 /* **************************************************
804 Implementation of the Li/Yeung recognition algorithm
806 ************************************************** */
808 #define WORST_SCORE 0x7fffffff
810 /* Dynamic programming parameters */
813 #define MAX_DIST 0x7fffffff
814 #define SIM_THLD 60 /* x 100 */
815 #define DIST_THLD 3200 /* x 100 */
817 /* Low-pass filter parameters -- empirically derived */
818 #define LP_FILTER_WIDTH 6
819 #define LP_FILTER_ITERS 8
820 #define LP_FILTER_THLD 250 /* x 100 */
821 #define LP_FILTER_MIN 5
823 /* Pseudo-extrema parameters -- empirically derived */
824 #define PE_AL_THLD 1500 /* x 100 */
825 #define PE_ATCR_THLD 135 /* x 100 */
827 /* Contour-angle derivation parameters */
831 /* Pre-processing and canonicalization parameters */
832 #define CANONICAL_X 108
833 #define CANONICAL_Y 128
834 #define DIST_SQ_THRESHOLD (3*3) /* copied from fv.h */
835 #define NCANONICAL 50
837 /* Tap-handling parameters */
839 #define TAP_TIME_THLD 150 /* msec */
840 #define TAP_DIST_THLD 75 /* dx * dx + dy * dy */
841 #define TAP_PATHLEN 1000 /* x 100 */
846 #define RGN_CONCAVE 1
851 typedef struct RegionList {
855 struct RegionList *next;
859 /* direction-code table; indexed by dx, dy */
860 static int lialg_dctbl[3][3] = {{1, 0, 7}, {2, 0x7FFFFFFF, 6}, {3, 4, 5}};
862 /* low-pass filter weights */
863 static int lialg_lpfwts[2 * LP_FILTER_WIDTH + 1];
864 static int lialg_lpfconst = -1;
867 static int lialg_preprocess_stroke(point_list *);
868 static point_list *lialg_compute_dominant_points(point_list *);
869 static point_list *lialg_interpolate_points(point_list *);
870 static void lialg_bresline(pen_point *, pen_point *, point_list *, int *);
871 static void lialg_compute_chain_code(point_list *);
872 static void lialg_compute_unit_chain_code(point_list *);
873 static region_list *lialg_compute_regions(point_list *);
874 static point_list *lialg_compute_dompts(point_list *, region_list *);
875 static int *lialg_compute_contour_angle_set(point_list *, region_list *);
876 static void lialg_score_stroke(point_list *, point_list *, int *, int *);
877 static int lialg_compute_similarity(point_list *, point_list *);
878 static int lialg_compute_distance(point_list *, point_list *);
880 static int lialg_read_classifier_digest(rClassifier *);
882 static int lialg_canonicalize_examples(rClassifier *);
883 static int lialg_canonicalize_example_stroke(point_list *);
884 static int lialg_compute_equipoints(point_list *);
886 static int lialg_compute_pathlen(point_list *);
887 static int lialg_compute_pathlen_subset(point_list *, int, int);
888 static int lialg_filter_points(point_list *);
889 static int lialg_translate_points(point_list *, int, int, int, int);
890 static void lialg_get_bounding_box(point_list *, int *, int *, int *, int *);
891 static void lialg_compute_lpf_parameters();
892 static int isqrt(int);
893 static int likeatan(int, int);
894 static int quadr(int);
897 /*************************************************************
899 Core routines for the Li/Yeung recognition algorithm
901 *************************************************************/
903 static void lialg_initialize(rClassifier *rec) {
906 /* Initialize the dompts arrays. */
907 for (i = 0; i < MAXSCLASSES; i++) {
908 rec->dompts[i] = nil;
914 * Main recognition routine -- called by HRE API.
916 static char *lialg_recognize_stroke(rClassifier *rec, point_list *stroke) {
919 point_list *input_dompts = nil;
920 char *best_name = nil;
921 int best_score = WORST_SCORE;
923 point_list *curr_dompts;
925 /* (void)gettimeofday(&stv, nil);*/
927 if (stroke->npts < 1) goto done;
931 /* First thing is to filter out ``close points.'' */
932 if (lialg_filter_points(stroke) != 0) return(nil);
934 /* Unfortunately, we don't have the actual time that each point */
935 /* was recorded (i.e., dt is invalid). Hence, we have to use a */
936 /* heuristic based on total distance and the number of points. */
937 if (stroke->npts == 1 || lialg_compute_pathlen(stroke) < TAP_PATHLEN)
940 /* Pre-process input stroke. */
941 if (lialg_preprocess_stroke(stroke) != 0) goto done;
943 /* Compute its dominant points. */
944 input_dompts = lialg_compute_dominant_points(stroke);
945 if (input_dompts == nil) goto done;
947 /* Score input stroke against every class in classifier. */
948 for (i = 0, curr_name = rec->cnames[i], curr_dompts = rec->dompts[i];
949 i < MAXSCLASSES && curr_name != nil && curr_dompts != nil;
950 i++, curr_name = rec->cnames[i], curr_dompts = rec->dompts[i]) {
955 lialg_score_stroke(input_dompts, curr_dompts, &sim, &dist);
958 if (lidebug && curr_score < DIST_THLD)
959 fprint(2, "(%s, %d, %d) ", curr_name, sim, dist);
961 /* Is it the best so far? */
962 if (curr_score < best_score && curr_score <= DIST_THLD) {
963 best_score = curr_score;
964 best_name = curr_name;
975 delete_examples(input_dompts);
980 static int lialg_preprocess_stroke(point_list *points) {
981 int minx, miny, maxx, maxy, xrange, yrange, scale, xoff, yoff;
983 /* Filter out points that are too close. */
984 /* We did this earlier, when we checked for a tap. */
986 if (lialg_filter_points(points) != 0) return(-1);
989 /* assert(points->npts > 0);*/
991 /* Scale up to avoid conversion errors. */
992 lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
993 xrange = maxx - minx;
994 yrange = maxy - miny;
995 scale = ( ((100 * xrange + CANONICAL_X / 2) / CANONICAL_X) >
996 ((100 * yrange + CANONICAL_Y / 2) / CANONICAL_Y))
997 ? (100 * CANONICAL_X + xrange / 2) / xrange
998 : (100 * CANONICAL_Y + yrange / 2) / yrange;
999 if (lialg_translate_points(points, minx, miny, scale, scale) != 0)
1002 /* Center the stroke. */
1003 lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
1004 xrange = maxx - minx;
1005 yrange = maxy - miny;
1006 xoff = -((CANONICAL_X - xrange + 1) / 2);
1007 yoff = -((CANONICAL_Y - yrange + 1) / 2);
1008 if (lialg_translate_points(points, xoff, yoff, 100, 100) != 0) return(-1);
1010 /* Store the x and y ranges in the point list. */
1011 xrange = maxx - minx;
1012 yrange = maxy - miny;
1013 points->xrange = xrange;
1014 points->yrange = yrange;
1018 fprint(2, "After pre-processing: %d %d %d %d\n",
1019 minx, miny, maxx, maxy);
1020 for (i = 0; i < points->npts; i++)
1021 fprint(2, " (%P)\n", points->pts[i].Point);
1029 static point_list *lialg_compute_dominant_points(point_list *points) {
1031 region_list *regions;
1034 /* Interpolate points. */
1035 ipts = lialg_interpolate_points(points);
1036 if (ipts == nil) return(nil);
1039 fprint(2, "After interpolation: %d ipts\n", ipts->npts);
1040 for (j = 0; j < ipts->npts; j++) {
1041 fprint(2, " (%P), %lud\n", ipts->pts[j].Point, ipts->pts[j].chaincode);
1046 /* Compute regions. */
1047 regions = lialg_compute_regions(ipts);
1048 /* assert(regions != nil);*/
1050 /* Compute dominant points. */
1051 dpts = lialg_compute_dompts(ipts, regions);
1054 fprint(2, "Dominant points: ");
1055 for (j = 0; j < dpts->npts; j++) {
1056 fprint(2, "%P (%lud) ", dpts->pts[j].Point, dpts->pts[j].chaincode);
1062 /* Delete region data structure. */
1064 region_list *curr, *next;
1065 for (curr = regions; curr != nil; ) {
1071 delete_examples(ipts);
1075 /* Input points are assumed to be integer-valued! */
1076 static point_list *lialg_interpolate_points(point_list *points) {
1081 /* Compute an upper-bound on the number of interpolated points. */
1083 for (i = 0; i < (points->npts - 1); i++) {
1084 pen_point *pta = &(points->pts[i]);
1085 pen_point *ptb = &(points->pts[i+1]);
1086 maxpts += abs(pta->x - ptb->x) + abs(pta->y - ptb->y);
1089 /* Allocate an array of the requisite size. */
1090 maxpts += points->npts;
1091 /* newpts = (point_list *)safe_malloc(sizeof(point_list)); */
1092 newpts = malloc(sizeof(point_list));
1093 newpts->pts = mallocz(maxpts*sizeof(pen_point), 1);
1094 if (newpts->pts == nil) {
1098 newpts->npts = maxpts;
1101 /* Interpolate each of the segments. */
1103 for (i = 0; i < (points->npts - 1); i++) {
1104 pen_point *startpt = &(points->pts[i]);
1105 pen_point *endpt = &(points->pts[i+1]);
1107 lialg_bresline(startpt, endpt, newpts, &j);
1109 j--; /* end point gets recorded as start point of next segment! */
1112 /* Add-in last point. */
1113 newpts->pts[j++] = points->pts[points->npts - 1];
1116 /* Compute the chain code for P (the list of points). */
1117 lialg_compute_unit_chain_code(newpts);
1123 /* This implementation is due to Kenny Hoff. */
1124 static void lialg_bresline(pen_point *startpt, pen_point *endpt,
1125 point_list *newpts, int *j) {
1126 int Ax, Ay, Bx, By, dX, dY, Xincr, Yincr;
1133 /* INITIALIZE THE COMPONENTS OF THE ALGORITHM THAT ARE NOT AFFECTED */
1134 /* BY THE SLOPE OR DIRECTION OF THE LINE */
1135 dX = abs(Bx-Ax); /* store the change in X and Y of the line endpoints */
1138 /* DETERMINE "DIRECTIONS" TO INCREMENT X AND Y (REGARDLESS OF DECISION) */
1139 if (Ax > Bx) { Xincr=-1; } else { Xincr=1; } /* which direction in X? */
1140 if (Ay > By) { Yincr=-1; } else { Yincr=1; } /* which direction in Y? */
1142 /* DETERMINE INDEPENDENT VARIABLE (ONE THAT ALWAYS INCREMENTS BY 1 (OR -1) ) */
1143 /* AND INITIATE APPROPRIATE LINE DRAWING ROUTINE (BASED ON FIRST OCTANT */
1144 /* ALWAYS). THE X AND Y'S MAY BE FLIPPED IF Y IS THE INDEPENDENT VARIABLE. */
1145 if (dX >= dY) { /* if X is the independent variable */
1146 int dPr = dY<<1; /* amount to increment decision if right is chosen (always) */
1147 int dPru = dPr - (dX<<1); /* amount to increment decision if up is chosen */
1148 int P = dPr - dX; /* decision variable start value */
1150 /* process each point in the line one at a time (just use dX) */
1151 for (; dX>=0; dX--) {
1152 newpts->pts[*j].x = Ax;
1153 newpts->pts[*j].y = Ay;
1156 if (P > 0) { /* is the pixel going right AND up? */
1157 Ax+=Xincr; /* increment independent variable */
1158 Ay+=Yincr; /* increment dependent variable */
1159 P+=dPru; /* increment decision (for up) */
1160 } else { /* is the pixel just going right? */
1161 Ax+=Xincr; /* increment independent variable */
1162 P+=dPr; /* increment decision (for right) */
1165 } else { /* if Y is the independent variable */
1166 int dPr = dX<<1; /* amount to increment decision if right is chosen (always) */
1167 int dPru = dPr - (dY<<1); /* amount to increment decision if up is chosen */
1168 int P = dPr - dY; /* decision variable start value */
1170 /* process each point in the line one at a time (just use dY) */
1171 for (; dY>=0; dY--) {
1172 newpts->pts[*j].x = Ax;
1173 newpts->pts[*j].y = Ay;
1176 if (P > 0) { /* is the pixel going up AND right? */
1177 Ax+=Xincr; /* increment dependent variable */
1178 Ay+=Yincr; /* increment independent variable */
1179 P+=dPru; /* increment decision (for up) */
1180 } else { /* is the pixel just going up? */
1181 Ay+=Yincr; /* increment independent variable */
1182 P+=dPr; /* increment decision (for right) */
1188 static void lialg_compute_chain_code(point_list *pts) {
1191 for (i = 0; i < (pts->npts - 1); i++) {
1192 pen_point *startpt = &(pts->pts[i]);
1193 pen_point *endpt = &(pts->pts[i+1]);
1194 int dx = endpt->x - startpt->x;
1195 int dy = endpt->y - startpt->y;
1196 int tmp = quadr(likeatan(dy, dx));
1197 int dircode = (12 - tmp) % 8;
1199 startpt->chaincode = dircode;
1204 static void lialg_compute_unit_chain_code(point_list *pts) {
1207 for (i = 0; i < (pts->npts - 1); i++) {
1208 pen_point *startpt = &(pts->pts[i]);
1209 pen_point *endpt = &(pts->pts[i+1]);
1210 int dx = endpt->x - startpt->x;
1211 int dy = endpt->y - startpt->y;
1212 int dircode = lialg_dctbl[dx+1][dy+1];
1214 startpt->chaincode = dircode;
1219 static region_list *lialg_compute_regions(point_list *pts) {
1220 region_list *regions;
1221 region_list *curr_reg;
1222 int *R[2 + LP_FILTER_ITERS];
1227 /* Initialize low-pass filter parameters if necessary. */
1228 if (lialg_lpfconst == -1)
1229 lialg_compute_lpf_parameters();
1231 /* Allocate a 2 x pts->npts array for use in computing the (filtered) Angle set, A_n. */
1232 junk = malloc((2 + LP_FILTER_ITERS) * pts->npts*sizeof(int));
1233 for (i = 0; i < (2 + LP_FILTER_ITERS); i++)
1234 R[i] = junk + (i * pts->npts);
1237 /* Compute the Angle set, A, in the first element of array R. */
1238 /* Values in R are in degrees, x 100. */
1239 curr[0] = 18000; /* a_0 */
1240 for (i = 1; i < (pts->npts - 1); i++) {
1241 int d_i = pts->pts[i].chaincode;
1242 int d_iminusone = pts->pts[i-1].chaincode;
1245 if (d_iminusone < d_i)
1248 a_i = (d_iminusone - d_i) % 8;
1250 /* convert to degrees, x 100 */
1251 curr[i] = ((12 - a_i) % 8) * 45 * 100;
1253 curr[pts->npts - 1] = 18000; /* a_L-1 */
1255 /* Perform a number of filtering iterations. */
1257 for (j = 0; j < LP_FILTER_ITERS; j++, curr = R[j], next = R[j+1]) {
1258 for (i = 0; i < pts->npts; i++) {
1263 for (k = i - LP_FILTER_WIDTH; k <= i + LP_FILTER_WIDTH; k++) {
1264 int oldval = (k < 0 || k >= pts->npts) ? 18000 : curr[k];
1265 next[i] += oldval * lialg_lpfwts[k - (i - LP_FILTER_WIDTH)]; /* overflow? */
1268 next[i] /= lialg_lpfconst;
1272 /* Do final thresholding around PI. */
1273 /* curr and next are set-up correctly at end of previous loop! */
1274 for (i = 0; i < pts->npts; i++)
1275 next[i] = (abs(curr[i] - 18000) < LP_FILTER_THLD) ? 18000 : curr[i];
1280 for (i = 0; i < pts->npts; i++) {
1281 fprint(2, "%3d: (%P) %lud ",
1282 i, pts->pts[i].Point, pts->pts[i].chaincode);
1283 for (j = 0; j < 2 + LP_FILTER_ITERS; j++)
1284 fprint(2, "%d ", R[j][i]);
1289 /* Do the region segmentation. */
1294 #define RGN_TYPE(val) (((val)==18000)?RGN_PLAIN:((val)<18000?RGN_CONCAVE:RGN_CONVEX))
1297 currtype = RGN_TYPE(curr[0]);
1298 regions = malloc(sizeof(region_list));
1300 curr_reg->start = start;
1302 curr_reg->type = currtype;
1303 curr_reg->next = nil;
1304 for (i = 1; i < pts->npts; i++) {
1305 int nexttype = RGN_TYPE(curr[i]);
1307 if (nexttype != currtype) {
1308 region_list *next_reg;
1311 curr_reg->end = end;
1313 fprint(2, " (%d, %d), %d\n", start, end, currtype);
1316 currtype = nexttype;
1317 next_reg = malloc(sizeof(region_list));
1318 next_reg->start = start;
1320 next_reg->type = nexttype;
1321 next_reg->next = nil;
1323 curr_reg->next = next_reg;
1324 curr_reg = next_reg;
1328 curr_reg->end = end;
1330 fprint(2, " (%d, %d), %d\n", start, end, currtype);
1332 /* Filter out convex/concave regions that are too short. */
1333 for (curr_reg = regions; curr_reg; curr_reg = curr_reg->next)
1334 if (curr_reg->type == RGN_PLAIN) {
1335 region_list *next_reg;
1337 for (next_reg = curr_reg->next;
1339 (next_reg->end - next_reg->start) < LP_FILTER_MIN;
1340 next_reg = curr_reg->next) {
1341 /* next_reg must not be plain, and it must be followed by a plain */
1342 /* assert(next_reg->type != RGN_PLAIN); */
1343 /* assert(next_reg->next != nil && (next_reg->next)->type == RGN_PLAIN); */
1345 curr_reg->next = (next_reg->next)->next;
1346 curr_reg->end = (next_reg->next)->end;
1348 free(next_reg->next);
1353 /* Add-in pseudo-extremes. */
1355 region_list *tmp, *prev_reg;
1360 for (curr_reg = tmp; curr_reg; curr_reg = curr_reg->next) {
1361 if (curr_reg->type == RGN_PLAIN) {
1362 int arclen = lialg_compute_pathlen_subset(pts,
1365 int dx = pts->pts[curr_reg->end].x -
1366 pts->pts[curr_reg->start].x;
1367 int dy = pts->pts[curr_reg->end].y -
1368 pts->pts[curr_reg->start].y;
1369 int chordlen = isqrt(10000 * (dx * dx + dy * dy));
1370 int atcr = (chordlen == 0) ? 0 : (100 * arclen + chordlen / 2) / chordlen;
1373 fprint(2, "%d, %d, %d\n", arclen, chordlen, atcr);
1375 /* Split region if necessary. */
1376 if (arclen >= PE_AL_THLD && atcr >= PE_ATCR_THLD) {
1377 int mid = curr_reg->start + (curr_reg->end - curr_reg->start) / 2;
1378 int end = curr_reg->end;
1379 region_list *saved_next = curr_reg->next;
1381 curr_reg->end = mid - 1;
1382 if (prev_reg == nil)
1385 prev_reg->next = curr_reg;
1386 prev_reg = curr_reg;
1388 /* curr_reg = (region_list *)safe_malloc(sizeof(region_list));*/
1389 curr_reg = malloc(sizeof(region_list));
1390 curr_reg->start = mid;
1391 curr_reg->end = mid;
1392 curr_reg->type = RGN_PSEUDO;
1393 curr_reg->next = nil;
1394 prev_reg->next = curr_reg;
1395 prev_reg = curr_reg;
1397 /* curr_reg = (region_list *)malloc(sizeof(region_list)); */
1398 curr_reg = malloc(sizeof(region_list));
1399 curr_reg->start = mid + 1;
1400 curr_reg->end = end;
1401 curr_reg->type = RGN_PLAIN;
1402 curr_reg->next = nil;
1403 prev_reg->next = curr_reg;
1404 prev_reg = curr_reg;
1406 curr_reg->next = saved_next;
1411 if (prev_reg == nil)
1414 prev_reg->next = curr_reg;
1415 prev_reg = curr_reg;
1425 static point_list *lialg_compute_dompts(point_list *pts, region_list *regions) {
1436 /* Compute contour angle set. */
1437 cas = lialg_compute_contour_angle_set(pts, regions);
1439 /* Dominant points include: start_pt, end_pt, extrema_of_non_plain_regions, midpts of the preceding. */
1441 for (r = regions; r != nil; r = r->next)
1442 if (r->type != RGN_PLAIN)
1444 ndpts = 2 * (2 + nonplain) - 1;
1445 /* dpts = (point_list *)safe_malloc(sizeof(point_list)); */
1446 dpts = malloc(sizeof(point_list));
1447 dpts->pts = mallocz(ndpts*sizeof(pen_point), 1);
1448 if (dpts->pts == nil) {
1455 /* Pick out dominant points. */
1457 /* Record start point. */
1460 dpts->pts[dp++] = pts->pts[previx];
1462 for (curr = regions; curr != nil; curr = curr->next)
1463 if (curr->type != RGN_PLAIN) {
1465 int min_v = 0x7fffffff; /* maxint */
1470 for (i = curr->start; i <= curr->end; i++) {
1472 if (v > max_v) { max_v = v; max_ix = i; }
1473 if (v < min_v) { min_v = v; min_ix = i; }
1475 fprint(2, " %d\n", v);
1478 currix = (curr->type == RGN_CONVEX ? max_ix : min_ix);
1480 /* Record midpoint. */
1481 dpts->pts[dp++] = pts->pts[previx + (currix - previx) / 2];
1483 /* Record extreme point. */
1484 dpts->pts[dp++] = pts->pts[currix];
1489 /* Record last mid-point and end point. */
1490 currix = pts->npts - 1;
1491 dpts->pts[dp++] = pts->pts[previx + (currix - previx) / 2];
1492 dpts->pts[dp] = pts->pts[currix];
1494 /* Compute chain-code. */
1495 lialg_compute_chain_code(dpts);
1502 static int *lialg_compute_contour_angle_set(point_list *pts,
1503 region_list *regions) {
1505 region_list *curr_reg;
1508 V = malloc(pts->npts*sizeof(int));
1511 for (curr_reg = regions; curr_reg != nil; curr_reg = curr_reg->next) {
1512 for (i = curr_reg->start; i <= curr_reg->end; i++) {
1513 if (curr_reg->type == RGN_PLAIN) {
1516 /* For now, simply choose the mid-point. */
1517 int isMidPt = i == (curr_reg->start +
1518 (curr_reg->end - curr_reg->start) / 2);
1519 V[i] = (curr_reg->type == RGN_CONVEX)
1520 ? (isMidPt ? 18000 : 0)
1521 : (isMidPt ? 0 : 18000);
1525 V[pts->npts - 1] = 18000;
1532 * First compute the similarity between the two strings.
1533 * If it's above a threshold, compute the distance between
1534 * the two and return it as the ``score.''
1535 * Otherwise, return the constant WORST_SCORE.
1538 static void lialg_score_stroke(point_list *input_dompts, point_list *curr_dompts, int *sim, int *dist) {
1542 *sim = lialg_compute_similarity(input_dompts, curr_dompts);
1543 if (*sim < SIM_THLD) goto done;
1545 *dist = lialg_compute_distance(input_dompts, curr_dompts);
1549 fprint(2, "%d, %d\n", *sim, *dist);
1553 static int lialg_compute_similarity(point_list *input_dompts, point_list *curr_dompts) {
1561 /* A is the longer sequence, length N. */
1562 /* B is the shorter sequence, length M. */
1563 if (input_dompts->npts >= curr_dompts->npts) {
1565 N = input_dompts->npts;
1567 M = curr_dompts->npts;
1570 N = curr_dompts->npts;
1572 M = input_dompts->npts;
1575 /* Allocate and initialize the Gain matrix, G. */
1576 /* The size of G is M x (N + 1). */
1577 /* Note that row 0 is unused. */
1578 /* Similarities are x 10. */
1579 G = malloc(M*sizeof(int *));
1580 junk = malloc(M * (N + 1) * sizeof(int));
1581 for (i = 0; i < M; i++)
1582 G[i] = junk + (i * (N + 1));
1584 for (i = 1; i < M; i++) {
1585 int bval = B->pts[i-1].chaincode;
1587 /* Source column. */
1590 for (j = 1; j < N; j++) {
1591 int aval = A->pts[j-1].chaincode;
1592 int diff = abs(bval - aval);
1593 if (diff > 4) diff = 8 - diff;
1595 G[i][j] = (diff == 0)
1606 /* Do the DP algorithm. */
1607 /* Proceed in column order, from highest column to the lowest. */
1608 /* Within each column, proceed from the highest row to the lowest. */
1609 /* Skip the highest column. */
1610 for (j = N - 1; j >= 0; j--)
1611 for (i = M - 1; i > 0; i--) {
1612 int max = G[i][j + 1];
1615 int tmp = G[i + 1][j + 1];
1616 if (tmp > max) max = tmp;
1622 sim = (10 * G[1][0] + (N - 1) / 2) / (N - 1);
1632 static int lialg_compute_distance(point_list *input_dompts,
1633 point_list *curr_dompts) {
1643 /* A is the longer sequence, length N. */
1644 /* B is the shorter sequence, length M. */
1645 if (input_dompts->npts >= curr_dompts->npts) {
1647 N = input_dompts->npts;
1649 M = curr_dompts->npts;
1653 N = curr_dompts->npts;
1655 M = input_dompts->npts;
1658 /* Construct the helper vectors, BE and TE, which say for each column */
1659 /* what are the ``bottom'' and ``top'' rows of interest. */
1660 BE = malloc((N + 1)*sizeof(int));
1661 TE = malloc((N + 1)*sizeof(int));
1663 for (j = 1; j <= N; j++) {
1666 bot = j + (M - DP_BAND);
1667 if (bot > M) bot = M;
1670 top = j - (N - DP_BAND);
1671 if (top < 1) top = 1;
1675 /* Allocate and initialize the Cost matrix, C. */
1676 /* The size of C is (M + 1) x (N + 1). */
1677 /* Note that row and column 0 are unused. */
1678 /* Costs are x 100. */
1679 /* C = (int **)safe_malloc((M + 1) * sizeof(int *)); */
1680 C = malloc((M + 1)*sizeof( int *));
1681 junk = malloc((M + 1) * (N + 1)*sizeof(int));
1682 for (i = 0; i <= M; i++)
1683 C[i] = junk + (i * (N + 1));
1685 for (i = 1; i <= M; i++) {
1686 int bx = B->pts[i-1].x;
1687 int by = B->pts[i-1].y;
1689 for (j = 1; j <= N; j++) {
1690 int ax = A->pts[j-1].x;
1691 int ay = A->pts[j-1].y;
1694 int dist = isqrt(10000 * (dx * dx + dy * dy));
1700 /* Do the DP algorithm. */
1701 /* Proceed in column order, from highest column to the lowest. */
1702 /* Within each column, proceed from the highest row to the lowest. */
1703 for (j = N; j > 0; j--)
1704 for (i = M; i > 0; i--) {
1707 if (i > BE[j] || i < TE[j] || (j == N && i == M))
1712 int tmp = C[i][j+1];
1718 int tmp = C[i+1][j+1];
1725 int tmp = C[i+1][j];
1726 if (tmp < min) min = tmp;
1732 dist = (C[1][1] + N / 2) / N;
1734 if (C != nil) free(C);
1735 if (junk != nil) free(junk);
1736 if (BE != nil) free(BE);
1737 if (TE != nil) free(TE);
1742 /*************************************************************
1744 Digest-processing routines
1746 *************************************************************/
1748 static int lialg_read_classifier_digest(rClassifier *rec) {
1752 /* Try to open the corresponding digest file. */
1757 /* Get a copy of the filename, with some room on the end. */
1758 /* clx_path = safe_malloc(strlen(rec->file_name) + 5); */
1759 clx_path = malloc((strlen(rec->file_name) + 5) *sizeof(char));
1760 strcpy(clx_path, rec->file_name);
1762 /* Truncate the path after the last dot. */
1763 dot = strrchr(clx_path, '.');
1764 if (dot == nil) { free(clx_path); return(-1); }
1767 /* Append the classifier-digest extension. */
1768 strcat(clx_path, "clx");
1770 fp = fopen(clx_path, "r");
1779 /* Read-in the name and dominant points for each class. */
1780 for (nclasses = 0; !feof(fp); nclasses++) {
1781 point_list *dpts = nil;
1786 if (fscanf(fp, "%s %d", class, &npts) != 2) {
1787 if (feof(fp)) break;
1791 rec->cnames[nclasses] = strdup(class);
1793 /* Allocate a dominant-points list. */
1794 /* dpts = (point_list *)safe_malloc(sizeof(point_list)); */
1795 dpts = malloc(sizeof(point_list));
1796 dpts->pts = mallocz(npts*sizeof(pen_point), 1);
1797 if (dpts->pts == nil) goto failed;
1801 /* Read in each point. */
1802 for (j = 0; j < npts; j++) {
1805 if (fscanf(fp, "%d %d", &x, &y) != 2) goto failed;
1810 /* Compute the chain-code. */
1811 lialg_compute_chain_code(dpts);
1813 /* Store the list in the rec data structure. */
1814 rec->dompts[nclasses] = dpts;
1819 fprint(2, "read_classifier_digest failed...\n");
1820 for (; nclasses >= 0; nclasses--) {
1821 if (rec->cnames[nclasses] != nil) {
1822 free(rec->cnames[nclasses]);
1823 rec->cnames[nclasses] = nil;
1825 if (rec->dompts[nclasses] != nil) {
1826 delete_examples(rec->dompts[nclasses]);
1827 rec->dompts[nclasses] = nil;
1831 delete_examples(dpts);
1841 /*************************************************************
1843 Canonicalization routines
1845 *************************************************************/
1847 static int lialg_canonicalize_examples(rClassifier *rec) {
1852 fprint(2, "lialg_canonicalize_examples working on %s\n",
1855 /* Initialize canonical-example arrays. */
1856 for (i = 0; i < MAXSCLASSES; i++) {
1857 rec->canonex[i] = nil;
1860 /* Figure out number of classes. */
1862 nclasses < MAXSCLASSES && rec->cnames[nclasses] != nil;
1866 /* Canonicalize the examples for each class. */
1867 for (i = 0; i < nclasses; i++) {
1870 point_list *pts, *tmp, *avg;
1871 int maxxrange, maxyrange;
1872 int minx, miny, maxx, maxy;
1873 int avgxrange, avgyrange, avgxoff, avgyoff, avgscale;
1877 fprint(2, "lialg_canonicalize_examples working on class %s\n",
1880 /* Make a copy of the examples. */
1883 for (nex = 0; tmp != nil; nex++, tmp = tmp->next) {
1884 if ((pts = add_example(pts, tmp->npts, tmp->pts)) == nil) {
1885 delete_examples(pts);
1890 /* Canonicalize each example, and derive the max x and y ranges. */
1893 for (j = 0, tmp = pts; j < nex; j++, tmp = tmp->next) {
1894 if (lialg_canonicalize_example_stroke(tmp) != 0) {
1896 fprint(2, "lialg_canonicalize_example_stroke returned error\n");
1901 if (tmp->xrange > maxxrange) maxxrange = tmp->xrange;
1902 if (tmp->yrange > maxyrange) maxyrange = tmp->yrange;
1905 /* Normalize max ranges. */
1906 if (((100 * maxxrange + CANONICAL_X / 2) / CANONICAL_X) >
1907 ((100 * maxyrange + CANONICAL_Y / 2) / CANONICAL_Y)) {
1908 maxyrange = (maxyrange * CANONICAL_X + maxxrange / 2) / maxxrange;
1909 maxxrange = CANONICAL_X;
1912 maxxrange = (maxxrange * CANONICAL_Y + maxyrange / 2) / maxyrange;
1913 maxyrange = CANONICAL_Y;
1916 /* Re-scale each example to max ranges. */
1917 for (j = 0, tmp = pts; j < nex; j++, tmp = tmp->next) {
1918 int scalex = (tmp->xrange == 0) ? 100 : (100 * maxxrange + tmp->xrange / 2) / tmp->xrange;
1919 int scaley = (tmp->yrange == 0) ? 100 : (100 * maxyrange + tmp->yrange / 2) / tmp->yrange;
1920 if (lialg_translate_points(tmp, 0, 0, scalex, scaley) != 0) {
1921 delete_examples(pts);
1926 /* Average the examples; leave average in first example. */
1927 avg = pts; /* careful aliasing!! */
1928 for (k = 0; k < NCANONICAL; k++) {
1932 for (j = 0, tmp = pts; j < nex; j++, tmp = tmp->next) {
1933 xsum += tmp->pts[k].x;
1934 ysum += tmp->pts[k].y;
1937 avg->pts[k].x = (xsum + j / 2) / j;
1938 avg->pts[k].y = (ysum + j / 2) / j;
1941 /* Compute BB of averaged stroke and re-scale. */
1942 lialg_get_bounding_box(avg, &minx, &miny, &maxx, &maxy);
1943 avgxrange = maxx - minx;
1944 avgyrange = maxy - miny;
1945 avgscale = (((100 * avgxrange + CANONICAL_X / 2) / CANONICAL_X) >
1946 ((100 * avgyrange + CANONICAL_Y / 2) / CANONICAL_Y))
1947 ? (100 * CANONICAL_X + avgxrange / 2) / avgxrange
1948 : (100 * CANONICAL_Y + avgyrange / 2) / avgyrange;
1949 if (lialg_translate_points(avg, minx, miny, avgscale, avgscale) != 0) {
1950 delete_examples(pts);
1954 /* Re-compute the x and y ranges and center the stroke. */
1955 lialg_get_bounding_box(avg, &minx, &miny, &maxx, &maxy);
1956 avgxrange = maxx - minx;
1957 avgyrange = maxy - miny;
1958 avgxoff = -((CANONICAL_X - avgxrange + 1) / 2);
1959 avgyoff = -((CANONICAL_Y - avgyrange + 1) / 2);
1960 if (lialg_translate_points(avg, avgxoff, avgyoff, 100, 100) != 0) {
1961 delete_examples(pts);
1965 /* Create a point list to serve as the ``canonical representation. */
1966 if ((rec->canonex[i] = add_example(nil, avg->npts, avg->pts)) == nil) {
1967 delete_examples(pts);
1970 (rec->canonex[i])->xrange = maxx - minx;
1971 (rec->canonex[i])->yrange = maxy - miny;
1974 fprint(2, "%s, avgpts = %d\n", rec->cnames[i], avg->npts);
1975 for (j = 0; j < avg->npts; j++) {
1976 fprint(2, " (%P)\n", avg->pts[j].Point);
1980 /* Compute dominant points of canonical representation. */
1981 rec->dompts[i] = lialg_compute_dominant_points(avg);
1984 delete_examples(pts);
1988 for (i = 0; i < nclasses; i++) {
1989 char *best_name = lialg_recognize_stroke(rec, rec->canonex[i]);
1991 if (best_name != rec->cnames[i])
1992 fprint(2, "%s, best = %s\n", rec->cnames[i], best_name);
1999 static int lialg_canonicalize_example_stroke(point_list *points) {
2000 int minx, miny, maxx, maxy, xrange, yrange, scale;
2002 /* Filter out points that are too close. */
2003 if (lialg_filter_points(points) != 0) return(-1);
2005 /* Must be at least two points! */
2006 if (points->npts < 2) {
2008 fprint(2, "lialg_canonicalize_example_stroke: npts=%d\n",
2014 /* Scale up to avoid conversion errors. */
2015 lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
2016 xrange = maxx - minx;
2017 yrange = maxy - miny;
2018 scale = (((100 * xrange + CANONICAL_X / 2) / CANONICAL_X) >
2019 ((100 * yrange + CANONICAL_Y / 2) / CANONICAL_Y))
2020 ? (100 * CANONICAL_X + xrange / 2) / xrange
2021 : (100 * CANONICAL_Y + yrange / 2) / yrange;
2022 if (lialg_translate_points(points, minx, miny, scale, scale) != 0) {
2024 fprint(2, "lialg_translate_points (minx=%d,miny=%d,scale=%d) returned error\n", minx, miny, scale);
2029 /* Compute an equivalent stroke with equi-distant points. */
2030 if (lialg_compute_equipoints(points) != 0) return(-1);
2032 /* Re-translate the points to the origin. */
2033 lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
2034 if (lialg_translate_points(points, minx, miny, 100, 100) != 0) {
2036 fprint(2, "lialg_translate_points (minx=%d,miny=%d) returned error\n", minx, miny);
2041 /* Store the x and y ranges in the point list. */
2042 xrange = maxx - minx;
2043 yrange = maxy - miny;
2044 points->xrange = xrange;
2045 points->yrange = yrange;
2049 fprint(2, "Canonicalized: %d, %d, %d, %d\n", minx, miny, maxx, maxy);
2050 for (i = 0; i < points->npts; i++)
2051 fprint(2, " (%P)\n", points->pts[i].Point);
2059 static int lialg_compute_equipoints(point_list *points) {
2060 pen_point *equipoints = mallocz(NCANONICAL*sizeof(pen_point), 1);
2061 int nequipoints = 0;
2062 int pathlen = lialg_compute_pathlen(points);
2063 int equidist = (pathlen + (NCANONICAL - 1) / 2) / (NCANONICAL - 1);
2065 int dist_since_last_eqpt;
2066 int remaining_seglen;
2067 int dist_to_next_eqpt;
2069 if (equipoints == nil) {
2070 fprint(2, "can't allocate memory in lialg_compute_equipoints");
2075 fprint(2, "compute_equipoints: npts = %d, pathlen = %d, equidist = %d\n",
2076 points->npts, pathlen, equidist);
2080 /* First original point is an equipoint. */
2081 equipoints[0] = points->pts[0];
2083 dist_since_last_eqpt = 0;
2085 for (i = 1; i < points->npts; i++) {
2086 int dx1 = points->pts[i].x - points->pts[i-1].x;
2087 int dy1 = points->pts[i].y - points->pts[i-1].y;
2088 int endx = 100 * points->pts[i-1].x;
2089 int endy = 100 * points->pts[i-1].y;
2090 remaining_seglen = isqrt(10000 * (dx1 * dx1 + dy1 * dy1));
2091 dist_to_next_eqpt = equidist - dist_since_last_eqpt;
2093 while (remaining_seglen >= dist_to_next_eqpt) {
2095 /* x-coordinate stays the same */
2097 endy += dist_to_next_eqpt;
2099 endy -= dist_to_next_eqpt;
2102 int slope = (100 * dy1 + dx1 / 2) / dx1;
2103 int tmp = isqrt(10000 + slope * slope);
2104 int dx = (100 * dist_to_next_eqpt + tmp / 2) / tmp;
2105 int dy = (slope * dx + 50) / 100;
2107 if (dy < 0) dy = -dy;
2118 equipoints[nequipoints].x = (endx + 50) / 100;
2119 equipoints[nequipoints].y = (endy + 50) / 100;
2121 /* assert(nequipoints <= NCANONICAL);*/
2122 dist_since_last_eqpt = 0;
2123 remaining_seglen -= dist_to_next_eqpt;
2124 dist_to_next_eqpt = equidist;
2127 dist_since_last_eqpt += remaining_seglen;
2130 /* Take care of last equipoint. */
2131 if (nequipoints == NCANONICAL) {
2133 } else if (nequipoints == (NCANONICAL - 1)) {
2134 /* Make last original point the last equipoint. */
2135 equipoints[nequipoints] = points->pts[points->npts - 1];
2138 fprint(2,"lialg_compute_equipoints: nequipoints = %d\n",
2145 points->npts = NCANONICAL;
2146 delete_pen_point_array(points->pts);
2147 points->pts = equipoints;
2152 /*************************************************************
2156 *************************************************************/
2158 /* Result is x 100. */
2159 static int lialg_compute_pathlen(point_list *points) {
2160 return(lialg_compute_pathlen_subset(points, 0, points->npts - 1));
2164 /* Result is x 100. */
2165 static int lialg_compute_pathlen_subset(point_list *points,
2166 int start, int end) {
2171 for (i = start + 1; i <= end; i++) {
2172 int dx = points->pts[i].x - points->pts[i-1].x;
2173 int dy = points->pts[i].y - points->pts[i-1].y;
2174 int dist = isqrt(10000 * (dx * dx + dy * dy));
2182 /* Note that this does NOT update points->xrange and points->yrange! */
2183 static int lialg_filter_points(point_list *points) {
2185 pen_point *filtered_pts = mallocz(points->npts*sizeof(pen_point), 1);
2188 if (filtered_pts == nil) {
2189 fprint(2, "can't allocate memory in lialg_filter_points");
2193 filtered_pts[0] = points->pts[0];
2195 for (i = 1; i < points->npts; i++) {
2196 int j = filtered_npts - 1;
2197 int dx = points->pts[i].x - filtered_pts[j].x;
2198 int dy = points->pts[i].y - filtered_pts[j].y;
2199 int magsq = dx * dx + dy * dy;
2201 if (magsq >= DIST_SQ_THRESHOLD) {
2202 filtered_pts[filtered_npts] = points->pts[i];
2207 points->npts = filtered_npts;
2208 delete_pen_point_array(points->pts);
2209 points->pts = filtered_pts;
2214 /* scalex and scaley are x 100. */
2215 /* Note that this does NOT update points->xrange and points->yrange! */
2216 static int lialg_translate_points(point_list *points,
2218 int scalex, int scaley) {
2221 for (i = 0; i < points->npts; i++) {
2222 points->pts[i].x = ((points->pts[i].x - minx) * scalex + 50) / 100;
2223 points->pts[i].y = ((points->pts[i].y - miny) * scaley + 50) / 100;
2230 static void lialg_get_bounding_box(point_list *points,
2231 int *pminx, int *pminy,
2232 int *pmaxx, int *pmaxy) {
2233 int minx, miny, maxx, maxy;
2236 minx = maxx = points->pts[0].x;
2237 miny = maxy = points->pts[0].y;
2238 for (i = 1; i < points->npts; i++) {
2239 pen_point *pt = &(points->pts[i]);
2240 if (pt->x < minx) minx = pt->x;
2241 else if (pt->x > maxx) maxx = pt->x;
2242 if (pt->y < miny) miny = pt->y;
2243 else if (pt->y > maxy) maxy = pt->y;
2253 int wtvals[] = {100, 104, 117, 143, 189, 271, 422};
2255 static void lialg_compute_lpf_parameters(void) {
2258 for (i = LP_FILTER_WIDTH; i >= 0; i--) {
2259 // double x = 0.04 * (i * i);
2260 // double tmp = 100.0 * exp(x);
2261 // int wt = floor((double)tmp);
2263 lialg_lpfwts[LP_FILTER_WIDTH - i] = wt;
2264 lialg_lpfwts[LP_FILTER_WIDTH + i] = wt;
2267 for (i = 0; i < (2 * LP_FILTER_WIDTH + 1); i++) {
2268 lialg_lpfconst += lialg_lpfwts[i];
2273 /* Code from Joseph Hall (jnhall@sat.mot.com). */
2274 static int isqrt(int n) {
2276 register long k0, k1, nn;
2278 for (nn = i = n, k0 = 2; i > 0; i >>= 2, k0 <<= 1)
2282 k1 = (nn / k0 + k0) >> 1;
2283 if (((k0 ^ k1) & ~1) == 0)
2287 return (int) ((k1 + 1) >> 1);
2291 /* Helper routines from Mark Hayter. */
2292 static int likeatan(int tantop, int tanbot) {
2294 /* Use tan(theta)=top/bot --> order for t */
2295 /* t in range 0..0x40000 */
2297 if ((tantop == 0) && (tanbot == 0))
2301 t = (tantop << 16) / (abs(tantop) + abs(tanbot));
2305 if (tantop < 0) t = 0x40000 + t;
2311 static int quadr(int t) {
2312 return (8 - (((t + 0x4000) >> 15) & 7)) & 7;