sys/src/libscribble/li_recognizer.c

   1 /*
   2  *  li_recognizer.c
   3  *
   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
  10  *
  11  *
  12  *  Adapted from cmu_recognizer.c by Jay Kistler.
  13  *
  14  *  Where is the CMU copyright???? Gotta track it down - Jim Gettys
  15  *
  16  *  Credit to Dean Rubine, Jim Kempf, and Ari Rapkin.
  17  */
  18
  19 #include <u.h>
  20 #include <libc.h>
  21 #include <stdio.h>
  22 #include <draw.h>
  23 #include <scribble.h>
  24 #include "scribbleimpl.h"
  25
  26 #include "hre_internal.h"
  27 #include "li_recognizer_internal.h"
  28
  29 int lidebug = 0;
  30
  31 #define LI_MAGIC 0xACCBADDD
  32
  33 #define CHECK_LI_MAGIC(_a) \
  34   ((_a) != nil && ((li_recognizer*)(_a))->li_magic == LI_MAGIC)
  35
  36
  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);
  42
  43
  44 char* li_err_msg = nil;
  45
  46 #define bcopy(s1,s2,n) memmove(s2,s1,n)
  47
  48 /*Freeing classifier*/
  49
  50 static void
  51 free_rClassifier(rClassifier* rc);
  52
  53 /*
  54  * Point List Support
  55 */
  56
  57 static point_list*
  58 add_example(point_list* l,int npts,pen_point* pts)
  59 {
  60         pen_point* lpts = mallocz(npts*sizeof(pen_point), 1);
  61         point_list *p = malloc(sizeof(point_list));
  62
  63         p->npts = npts;
  64         p->pts = lpts;
  65         p->next = l;       /*Order doesn't matter, so we stick on end.*/
  66
  67         /*Copy points.*/
  68
  69         bcopy(pts, lpts, npts * sizeof(pen_point));
  70
  71         return(p);
  72 }
  73
  74
  75 static void
  76 delete_examples(point_list* l)
  77 {
  78         point_list* p;
  79
  80         for( ; l != nil; l = p ) {
  81                 p = l->next;
  82                 free(l->pts);
  83                 free(l);
  84         }
  85 }
  86
  87 /*
  88  * Local functions
  89  */
  90
  91 /*
  92  * recognize_internal-Form Vector, use Classifier to classify, return char.
  93  */
  94
  95 static char*
  96 recognize_internal(rClassifier* rec, Stroke* str, int*)
  97 {
  98         char *res;
  99         point_list *stroke;
 100
 101         stroke = add_example(nil, str->npts, str->pts);
 102         if (stroke == nil) return(nil);
 103
 104         res = lialg_recognize_stroke(rec, stroke);
 105
 106         delete_examples(stroke);
 107         return(res);
 108 }
 109
 110 /*
 111  * file_path-Construct pathname, check for proper extension.
 112  */
 113
 114 static int
 115   file_path(char* dir,char* filename,char* pathname)
 116 {
 117         char* dot;
 118
 119         /*Check for proper extension on file name.*/
 120
 121         dot = strrchr(filename,'.');
 122
 123         if( dot == nil ) {
 124                 return(-1);
 125         }
 126
 127         /*Determine whether a gesture or character classifier.*/
 128
 129         if( strcmp(dot,LI_CLASSIFIER_EXTENSION) != 0 ) {
 130                 return(-1);
 131         }
 132
 133         /*Concatenate directory and filename into pathname.*/
 134
 135         strcpy(pathname,dir);
 136         strcat(pathname,"/");
 137         strcat(pathname,filename);
 138
 139         return(0);
 140 }
 141
 142 /*read_classifier_points-Read points so classifier can be extended.*/
 143
 144 static int
 145 read_classifier_points(FILE* fd,int nclss,point_list** ex,char** cnames)
 146 {
 147         int i,j,k;
 148         char buf[BUFSIZ];
 149         int nex = 0;
 150         char* names[MAXSCLASSES];
 151         point_list* examples[MAXSCLASSES];
 152         pen_point* pts;
 153         int npts;
 154
 155         /*Initialize*/
 156
 157                 for( i = 0; i < MAXSCLASSES; i++ ) {
 158                                 names[i] = nil;
 159                                 examples[i] = nil;
 160                 }
 161
 162                 /*Go thru classes.*/
 163
 164                 for( k = 0; k < nclss; k++ ) {
 165
 166                                 /*Read class name and number of examples.*/
 167
 168                                 if( fscanf(fd,"%d %s",&nex,buf) != 2 )
 169                                         goto unallocate;
 170
 171                                 /*Save class name*/
 172
 173                                 names[k] = strdup(buf);
 174
 175                                 /*Read examples.*/
 176
 177                                 for( i = 0; i < nex; i++ ) {
 178
 179                                         /*Read number of points.*/
 180
 181                                         if( fscanf(fd,"%d",&npts) != 1 )
 182                                                                 goto unallocate; /*Boy would I like exceptions!*/
 183
 184                                         /*Allocate array for points.*/
 185
 186                                         if( (pts = mallocz(npts*sizeof(pen_point), 1)) == nil )
 187                                                                 goto unallocate;
 188
 189                                         /*Read in points.*/
 190
 191                                         for( j = 0; j < npts; j++ ) {
 192                                                                 int x,y;
 193                                                                 if( fscanf(fd,"%d %d",&x,&y) != 2 ) {
 194                                                                         delete_pen_point_array(pts);
 195                                                                         goto unallocate;
 196                                                                 }
 197                                                                 pts[j].Point = Pt(x, y);
 198                                         }
 199
 200                                         /*Add example*/
 201
 202                                         if( (examples[k] = add_example(examples[k],npts,pts)) == nil ) {
 203                                                                 delete_pen_point_array(pts);
 204                                                                 goto unallocate;
 205                                         }
 206
 207                                         delete_pen_point_array(pts);
 208
 209                                 }
 210                 }
 211
 212 /* ari -- end of list of classes */
 213 /* fprint(2, "]\n"); */
 214
 215                 /*Transfer to recognizer.*/
 216
 217                 bcopy(examples,ex,sizeof(examples));
 218                 bcopy(names,cnames,sizeof(names));
 219
 220                 return(0);
 221
 222                 /*Error. Deallocate memory and return.*/
 223
 224 unallocate:
 225                 for( ; k >= 0; k-- ) {
 226                                 delete_examples(examples[k]);
 227                                 free(names[k]);
 228                 }
 229
 230                 return(-1);
 231 }
 232
 233 /*read_classifier-Read a classifier file.*/
 234
 235 static int read_classifier(FILE* fd,rClassifier* rc)
 236 {
 237
 238         li_err_msg = nil;
 239
 240         /*Read in classifier file.*/
 241
 242         if(fscanf(fd, "%d", &rc->nclasses) != 1)
 243                 return -1;
 244
 245         /*Read in the example points, so classifier can be extended.*/
 246
 247         if( read_classifier_points(fd,rc->nclasses,rc->ex,rc->cnames) != 0 )
 248                 return -1;
 249
 250         return(0);
 251 }
 252
 253 /*
 254  * Extension Functions
 255 */
 256
 257 /* getClasses and clearState are by Ari */
 258
 259 static int
 260 recognizer_getClasses (recognizer r, char ***list, int *nc)
 261 {
 262         int i, nclasses;
 263         li_recognizer* rec;
 264         char **ret;
 265
 266                 rec = (li_recognizer*)r->recognizer_specific;
 267
 268                 /*Check for LI recognizer.*/
 269
 270                 if( !CHECK_LI_MAGIC(rec) ) {
 271                                 li_err_msg = "Not a LI recognizer";
 272                                 return(-1);
 273         }
 274
 275         *nc = nclasses = rec->li_rc.nclasses;
 276         ret = malloc(nclasses*sizeof(char*));
 277
 278         for (i = 0; i < nclasses; i++)
 279                                 ret[i] = rec->li_rc.cnames[i];   /* only the 1st char of the cname */
 280         *list = ret;
 281                 return 0;
 282 }
 283
 284 static int
 285 recognizer_clearState (recognizer)
 286 {
 287   /*This operation isn't supported by the LI recognizer.*/
 288
 289   li_err_msg = "Clearing state is not supported by the LI recognizer";
 290
 291   return(-1);
 292 }
 293
 294 static bool isa_li(recognizer r)
 295 { return(CHECK_LI_MAGIC(r)); }
 296
 297 static int
 298 recognizer_train(recognizer, rc*, uint, Stroke*, rec_element*, bool)
 299 {
 300   /*This operation isn't supported by the LI recognizer.*/
 301
 302   li_err_msg = "Training is not supported by the LI recognizer";
 303
 304   return(-1);
 305 }
 306
 307 int
 308 li_recognizer_get_example (recognizer r,
 309                                                    int          class,
 310                                                    int          instance,
 311                                                    char         **name,
 312                                                    pen_point    **points,
 313                                                    int          *npts)
 314 {
 315         li_recognizer   *rec = (li_recognizer*)r->recognizer_specific;
 316         int nclasses = rec->li_rc.nclasses;
 317         point_list          *pl;
 318
 319                 if( !CHECK_LI_MAGIC(rec) ) {
 320                                 li_err_msg = "Not a LI recognizer";
 321                                 return(-1);
 322                 }
 323                 if (class > nclasses)
 324                                 return -1;
 325                 pl = rec->li_rc.canonex[class];
 326                 while (instance && pl)
 327                 {
 328                                 pl = pl->next;
 329                                 instance--;
 330                 }
 331                 if (!pl)
 332                                 return -1;
 333                 *name = rec->li_rc.cnames[class];
 334                 *points = pl->pts;
 335                 *npts = pl->npts;
 336                 return pl->npts; /* I hope [sjm] */
 337 }
 338
 339 /*
 340  * API Functions
 341  */
 342
 343
 344 /*li_recognizer_load-Load a classifier file.*/
 345
 346 static int li_recognizer_load(recognizer r, char* dir, char* filename)
 347 {
 348                 FILE *fd;
 349                 char* pathname;
 350                 li_recognizer* rec;
 351                 rClassifier* rc;
 352
 353                 rec = (li_recognizer*)r->recognizer_specific;
 354
 355                 /*Make sure recognizer's OK*/
 356
 357                 if( !CHECK_LI_MAGIC(rec) ) {
 358                                 li_err_msg = "Not a LI recognizer";
 359                                 return(-1);
 360                 }
 361
 362                 rc = &(rec->li_rc);
 363
 364                 /*Check parameters.*/
 365
 366                 if( filename == nil ) {
 367                                 li_err_msg = "Invalid parameters";
 368                                 return(-1);
 369                 }
 370
 371                 /*We let the directory be null.*/
 372
 373                 if( dir == nil || (int)strlen(dir) <= 0 ) {
 374                                 dir = ".";
 375                 }
 376
 377                 if(0)fprint(2, "dir = %s, filename = %s\n",
 378                                 dir, filename);
 379
 380                 /*Make full pathname and check filename*/
 381
 382                 pathname = malloc((strlen(dir) + strlen(filename) + 2)*sizeof(char));
 383                 if( file_path(dir, filename, pathname) == -1 ) {
 384                                 free(pathname);
 385                                 li_err_msg = "Not a LI recognizer classifier file";
 386                                 return -1;
 387                 }
 388
 389                 /* Try to short-circuit the full classifier-file processing. */
 390                 rc->file_name = pathname;
 391                 if (lialg_read_classifier_digest(rc) == 0)
 392                                 return(0);
 393                 rc->file_name = nil;
 394
 395                 /*Open the file*/
 396
 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);
 400                                 free(pathname);
 401                                 return(-1);
 402                 }
 403
 404                 /*If rClassifier is OK, then delete it first.*/
 405
 406                 if( rc->file_name != nil ) {
 407                                 free_rClassifier(rc);
 408                 }
 409
 410                 /*Read classifier.*/
 411
 412                 if( read_classifier(fd,rc) < 0 ) {
 413                                 free(pathname);
 414                                 return(-1);
 415                 }
 416
 417                 /*Close file.*/
 418
 419                 fclose(fd);
 420
 421                 /*Add classifier name.*/
 422
 423                 rc->file_name = pathname;
 424
 425                 /* Canonicalize examples. */
 426                 if (lialg_canonicalize_examples(rc) != 0) {
 427                                 free(pathname);
 428                                 rc->file_name = nil;
 429                                 return -1;
 430                 }
 431
 432                 return(0);
 433 }
 434
 435 /*li_recognizer_save-Save a classifier file.*/
 436
 437 static int li_recognizer_save(recognizer, char*, char*)
 438 {
 439                 /*This operation isn't supported by the LI recognizer.*/
 440
 441                 li_err_msg = "Saving is not supported by the LI recognizer";
 442                 return -1;
 443 }
 444
 445 static wordset
 446 li_recognizer_load_dictionary(recognizer, char*, char*)
 447 {
 448                 /*This operation isn't supported by the LI recognizer.*/
 449
 450                 li_err_msg = "Dictionaries are not supported by the LI recognizer";
 451                 return nil;
 452 }
 453
 454 static int
 455 li_recognizer_save_dictionary(recognizer, char*, char*, wordset)
 456 {
 457                 /*This operation isn't supported by the LI recognizer.*/
 458                 li_err_msg = "Dictionaries are not supported by the LI recognizer";
 459                 return -1;
 460 }
 461
 462 static int
 463 li_recognizer_free_dictionary(recognizer, wordset)
 464 {
 465   /*This operation isn't supported by the LI recognizer.*/
 466
 467   li_err_msg = "Dictionaries are not supported by the LI recognizer";
 468
 469   return -1;
 470
 471 }
 472
 473 static int
 474 li_recognizer_add_to_dictionary(recognizer, letterset*, wordset)
 475 {
 476                 /*This operation isn't supported by the LI recognizer.*/
 477                 li_err_msg = "Dictionaries are not supported by the LI recognizer";
 478                 return -1;
 479 }
 480
 481 static int
 482 li_recognizer_delete_from_dictionary(recognizer, letterset*, wordset)
 483 {
 484                 /*This operation isn't supported by the LI recognizer.*/
 485                 li_err_msg = "Dictionaries are not supported by the LI recognizer";
 486                 return -1;
 487 }
 488
 489 static char*
 490 li_recognizer_error(recognizer rec)
 491 {
 492         char* ret = li_err_msg;
 493
 494         /*Check for LI recognizer.*/
 495
 496         if( !CHECK_LI_MAGIC(rec->recognizer_specific) ) {
 497                                 li_err_msg = "Not a LI recognizer";
 498                                 return nil;
 499         }
 500         li_err_msg = nil;
 501         return ret;
 502 }
 503
 504 static int
 505 li_recognizer_clear(recognizer r, bool)
 506 {
 507                 li_recognizer* rec;
 508
 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";
 513                                 return 0;
 514                 }
 515                 return 0;
 516 }
 517
 518 static int
 519 li_recognizer_set_context(recognizer, rc*)
 520 {
 521                 /*This operation isn't supported by the LI recognizer.*/
 522                 li_err_msg = "Contexts are not supported by the LI recognizer";
 523                 return -1;
 524 }
 525
 526 static rc*
 527 li_recognizer_get_context(recognizer)
 528 {
 529                 /*This operation isn't supported by the LI recognizer.*/
 530                 li_err_msg = "Contexts are not supported by the LI recognizer";
 531                 return nil;
 532 }
 533
 534 static int
 535 li_recognizer_get_buffer(recognizer, uint*, Stroke**)
 536 {
 537                 /*This operation isn't supported by the LI recognizer.*/
 538                 li_err_msg = "Buffer get/set are not supported by the LI recognizer";
 539                 return -1;
 540 }
 541
 542 static int
 543 li_recognizer_set_buffer(recognizer, uint, Stroke*)
 544 {
 545                 /*This operation isn't supported by the LI recognizer.*/
 546                 li_err_msg = "Buffer get/set are not supported by the LI recognizer";
 547                 return -1;
 548 }
 549
 550 static int
 551 li_recognizer_translate(recognizer r, uint ncs, Stroke* tps, bool, int* nret, rec_alternative** ret)
 552 {
 553         char* clss;
 554         li_recognizer* rec;
 555         int conf;
 556         rClassifier* rc;
 557
 558         rec = (li_recognizer*)r->recognizer_specific;
 559
 560         *nret = 0;
 561         *ret = nil;
 562
 563         /*Check for LI recognizer.*/
 564
 565         if( !CHECK_LI_MAGIC(rec) ) {
 566                                 li_err_msg = "Not a LI recognizer";
 567                                 return(-1);
 568         }
 569
 570                 rc = &(rec->li_rc);
 571
 572         /*Check for valid parameters.*/
 573         if (ncs < 1) {
 574                                 li_err_msg = "Invalid parameters: ncs";
 575                                 return(-1);
 576         }
 577         if( tps == nil) {
 578                                 li_err_msg = "Invalid parameters: tps";
 579                                 return(-1);
 580         }
 581         if( nret == nil) {
 582                                 li_err_msg = "Invalid parameters: nret";
 583                                 return(-1);
 584         }
 585         if( ret == nil) {
 586                                 li_err_msg = "Invalid parameters: ret";
 587                                 return(-1);
 588         }
 589
 590         /*
 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.
 596         */
 597         clss = recognize_internal(rc,tps,&conf);
 598         if (clss == nil) {
 599                                 *nret = 1;
 600                                 return(0);
 601         }
 602
 603         /*Return values.*/
 604         *nret = 1;
 605         return(*clss);
 606 }
 607
 608
 609 static rec_fn*
 610 li_recognizer_get_extension_functions(recognizer rec)
 611 {
 612         rec_fn* ret;
 613
 614         /*Check for LI recognizer.*/
 615
 616         if( !CHECK_LI_MAGIC(rec->recognizer_specific) ) {
 617                                 li_err_msg = "Not a LI recognizer";
 618                                 return(nil);
 619         }
 620
 621         ret = make_rec_fn_array(LI_NUM_EX_FNS);
 622
 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;
 628
 629         return(ret);
 630 }
 631
 632 static char**
 633 li_recognizer_get_gesture_names(recognizer)
 634 {
 635                 /*This operation isn't supported by the LI recognizer.*/
 636                 li_err_msg = "Gestures are not supported by the LI recognizer";
 637                 return nil;
 638 }
 639
 640 static xgesture
 641 li_recognizer_set_gesture_action(recognizer, char*, xgesture, void*)
 642 {
 643                 /*This operation isn't supported by the LI recognizer.*/
 644                 li_err_msg = "Gestures are not supported by the LI recognizer";
 645                 return nil;
 646 }
 647
 648 /*
 649  * Exported Functions
 650 */
 651
 652 /*RECOGNIZER_INITIALIZE-Initialize the recognizer.*/
 653
 654 /* note from ari:  this expands via pre-processor to
 655  *
 656  * recognizer __recognizer_internal_initialize(rec_info* ri)
 657  */
 658
 659 RECOGNIZER_INITIALIZE(ri)
 660 {
 661         recognizer r;
 662         li_recognizer* rec;
 663         int i;
 664
 665         /*Check that locale matches.*/
 666
 667         if( strcmp(ri->ri_locale,LI_SUPPORTED_LOCALE) != 0 ) {
 668                 li_err_msg = "Not a supported locale";
 669 /* fprint(2, "Locale error.\n");*/
 670                 return(nil);
 671         }
 672
 673         /*
 674          * Check that character sets match. Note that this is only approximate,
 675          * since the classifier file will have more information.
 676         */
 677
 678         if( ri->ri_subset != nil ) {
 679           for(i = 0; ri->ri_subset[i] != nil; i++ ) {
 680
 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"); */
 687
 688                   return(nil);
 689                 }
 690           }
 691         }
 692
 693 /* ari */
 694         r = make_recognizer(ri);
 695 /* fprint(2, "past make_recognizer.\n"); */
 696
 697         if( r == nil ) {
 698                 li_err_msg = "Can't allocate storage";
 699
 700                 return(nil);
 701         }
 702
 703         /*Make a LI recognizer structure.*/
 704
 705
 706         /* rec = (li_recognizer*)safe_malloc(sizeof(li_recognizer))) == nil ); */
 707
 708         rec = malloc(sizeof(li_recognizer));
 709
 710         r->recognizer_specific = rec;
 711
 712         rec->li_rc.file_name = nil;
 713         rec->li_rc.nclasses = 0;
 714
 715         /*Initialize the recognizer struct.*/
 716
 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;
 736
 737         /*Initialize LI Magic Number.*/
 738
 739         rec->li_magic = LI_MAGIC;
 740
 741         /*Initialize rClassifier.*/
 742
 743         rec->li_rc.file_name = nil;
 744
 745         for( i = 0; i < MAXSCLASSES; i++ ) {
 746                 rec->li_rc.ex[i] = nil;
 747                 rec->li_rc.cnames[i] = nil;
 748         }
 749
 750         lialg_initialize(&rec->li_rc);
 751
 752         /*Get rid of error message. Not needed here.*/
 753         li_err_msg = nil;
 754
 755         return(r);
 756 }
 757
 758 /*free_rClassifier-Free the rClassifier.*/
 759
 760 static void
 761 free_rClassifier(rClassifier* rc)
 762 {
 763         int i;
 764
 765         if( rc->file_name != nil) {
 766                 free(rc->file_name);
 767         }
 768
 769         for( i = 0; rc->ex[i] != nil; i++) {
 770                 delete_examples(rc->ex[i]);
 771                 free(rc->cnames[i]);
 772         }
 773
 774 }
 775
 776 /*RECOGNIZER_FINALIZE-Deallocate the recognizer, finalize.*/
 777
 778 RECOGNIZER_FINALIZE(r)
 779 {
 780                 li_recognizer* rec = (li_recognizer*)r->recognizer_specific;
 781
 782                 /*Make sure this is a li_recognizer first*/
 783                 if( !CHECK_LI_MAGIC(rec) ) {
 784                                 li_err_msg = "Not a LI recognizer";
 785                                 return(-1);
 786         }
 787
 788                 /*Deallocate rClassifier resources.*/
 789
 790                 free_rClassifier(&(rec->li_rc));
 791
 792                 /*Deallocate the li_recognizer struct.*/
 793                 free(rec);
 794
 795                 /*Deallocate the recognizer*/
 796                 delete_recognizer(r);
 797
 798                 return(0);
 799 }
 800
 801
 802 /* **************************************************
 803
 804   Implementation of the Li/Yeung recognition algorithm
 805
 806 ************************************************** */
 807
 808 #define WORST_SCORE     0x7fffffff
 809
 810 /* Dynamic programming parameters */
 811 #define DP_BAND         3
 812 #define MIN_SIM         0
 813 #define MAX_DIST        0x7fffffff
 814 #define SIM_THLD        60      /* x 100 */
 815 #define DIST_THLD       3200    /* x 100 */
 816
 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
 822
 823 /* Pseudo-extrema parameters -- empirically derived */
 824 #define PE_AL_THLD      1500    /* x 100 */
 825 #define PE_ATCR_THLD    135     /* x 100 */
 826
 827 /* Contour-angle derivation parameters */
 828 #define T_ONE           1
 829 #define T_TWO           20
 830
 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
 836
 837 /* Tap-handling parameters */
 838 #define TAP_CHAR        "."
 839 #define TAP_TIME_THLD   150         /* msec */
 840 #define TAP_DIST_THLD   75          /* dx * dx + dy * dy */
 841 #define TAP_PATHLEN     1000        /* x 100 */
 842
 843
 844 /* region types */
 845 #define RGN_CONVEX  0
 846 #define RGN_CONCAVE 1
 847 #define RGN_PLAIN   2
 848 #define RGN_PSEUDO  3
 849
 850
 851 typedef struct RegionList {
 852         int start;
 853         int end;
 854         int type;
 855         struct RegionList *next;
 856 } region_list;
 857
 858
 859 /* direction-code table; indexed by dx, dy */
 860 static int lialg_dctbl[3][3] = {{1, 0, 7}, {2, 0x7FFFFFFF, 6}, {3, 4, 5}};
 861
 862 /* low-pass filter weights */
 863 static int lialg_lpfwts[2 * LP_FILTER_WIDTH + 1];
 864 static int lialg_lpfconst = -1;
 865
 866
 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 *);
 879
 880 static int lialg_read_classifier_digest(rClassifier *);
 881
 882 static int lialg_canonicalize_examples(rClassifier *);
 883 static int lialg_canonicalize_example_stroke(point_list *);
 884 static int lialg_compute_equipoints(point_list *);
 885
 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);
 895
 896
 897 /*************************************************************
 898
 899   Core routines for the Li/Yeung recognition algorithm
 900
 901  *************************************************************/
 902
 903 static void lialg_initialize(rClassifier *rec) {
 904         int i;
 905
 906         /* Initialize the dompts arrays. */
 907         for (i = 0; i < MAXSCLASSES; i++) {
 908                 rec->dompts[i] = nil;
 909         }
 910 }
 911
 912
 913 /*
 914  *  Main recognition routine -- called by HRE API.
 915  */
 916 static char *lialg_recognize_stroke(rClassifier *rec, point_list *stroke) {
 917                 int i;
 918                 char *name = nil;
 919                 point_list *input_dompts = nil;
 920                 char *best_name = nil;
 921                 int best_score = WORST_SCORE;
 922                 char *curr_name;
 923                 point_list *curr_dompts;
 924
 925                 /*    (void)gettimeofday(&stv, nil);*/
 926
 927                 if (stroke->npts < 1) goto done;
 928
 929                 /* Check for tap. */
 930
 931                 /* First thing is to filter out ``close points.'' */
 932                 if (lialg_filter_points(stroke) != 0) return(nil);
 933
 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)
 938                         return(TAP_CHAR);
 939
 940                 /* Pre-process input stroke. */
 941                 if (lialg_preprocess_stroke(stroke) != 0) goto done;
 942
 943                 /* Compute its dominant points. */
 944                 input_dompts = lialg_compute_dominant_points(stroke);
 945                 if (input_dompts == nil) goto done;
 946
 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]) {
 951                                 int sim;
 952                                 int dist;
 953                                 int curr_score;
 954
 955                                 lialg_score_stroke(input_dompts, curr_dompts, &sim, &dist);
 956                                 curr_score = dist;
 957
 958                                 if (lidebug && curr_score < DIST_THLD)
 959                                 fprint(2, "(%s, %d, %d) ", curr_name, sim, dist);
 960
 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;
 965                                 }
 966                 }
 967
 968                 if (lidebug)
 969                                 fprint(2, "\n");
 970
 971                 /* No errors. */
 972                 name = best_name;
 973
 974 done:
 975                 delete_examples(input_dompts);
 976                 return(name);
 977 }
 978
 979
 980 static int lialg_preprocess_stroke(point_list *points) {
 981         int minx, miny, maxx, maxy, xrange, yrange, scale, xoff, yoff;
 982
 983         /* Filter out points that are too close. */
 984         /* We did this earlier, when we checked for a tap. */
 985 /*
 986         if (lialg_filter_points(points) != 0) return(-1);
 987 */
 988
 989 /*    assert(points->npts > 0);*/
 990
 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)
1000                 return(-1);
1001
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);
1009
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;
1015
1016         if (lidebug) {
1017                 int i;
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);
1022                 fflush(stderr);
1023         }
1024
1025         return(0);
1026 }
1027
1028
1029 static point_list *lialg_compute_dominant_points(point_list *points) {
1030         point_list *ipts;
1031         region_list *regions;
1032         point_list *dpts;
1033
1034         /* Interpolate points. */
1035         ipts = lialg_interpolate_points(points);
1036         if (ipts == nil) return(nil);
1037         if (lidebug) {
1038                                 int j;
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);
1042                                 }
1043                                 fflush(stderr);
1044         }
1045
1046         /* Compute regions. */
1047         regions = lialg_compute_regions(ipts);
1048 /*    assert(regions != nil);*/
1049
1050         /* Compute dominant points. */
1051         dpts = lialg_compute_dompts(ipts, regions);
1052         if (lidebug) {
1053                                 int j;
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);
1057                                 }
1058                                 fprint(2, "\n");
1059                                 fflush(stderr);
1060         }
1061
1062         /* Delete region data structure. */
1063         {
1064                                 region_list *curr, *next;
1065                                 for (curr = regions; curr != nil; ) {
1066                                         next = curr->next;
1067                                         free(curr);
1068                                         curr = next;
1069                                 }
1070         }
1071         delete_examples(ipts);
1072         return(dpts);
1073 }
1074
1075 /* Input points are assumed to be integer-valued! */
1076 static point_list *lialg_interpolate_points(point_list *points) {
1077         int i, j;
1078         int maxpts;
1079         point_list *newpts;
1080
1081         /* Compute an upper-bound on the number of interpolated points. */
1082         maxpts = 0;
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);
1087         }
1088
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) {
1095                                 free(newpts);
1096                                 return(nil);
1097         }
1098         newpts->npts = maxpts;
1099         newpts->next = nil;
1100
1101         /* Interpolate each of the segments. */
1102         j = 0;
1103         for (i = 0; i < (points->npts - 1); i++) {
1104                                 pen_point *startpt = &(points->pts[i]);
1105                                 pen_point *endpt = &(points->pts[i+1]);
1106
1107                                 lialg_bresline(startpt, endpt, newpts, &j);
1108
1109                                 j--;    /* end point gets recorded as start point of next segment! */
1110         }
1111
1112         /* Add-in last point. */
1113         newpts->pts[j++] = points->pts[points->npts - 1];
1114         newpts->npts = j;
1115
1116         /* Compute the chain code for P (the list of points). */
1117         lialg_compute_unit_chain_code(newpts);
1118
1119         return(newpts);
1120 }
1121
1122
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;
1127
1128         Ax = startpt->x;
1129         Ay = startpt->y;
1130         Bx = endpt->x;
1131         By = endpt->y;
1132
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 */
1136         dY = abs(By-Ay);
1137
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? */
1141
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 */
1149
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;
1154                                         (*j)++;
1155
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) */
1163                                         }
1164                                 }
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 */
1169
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;
1174                                         (*j)++;
1175
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) */
1183                                         }
1184                                 }
1185         }
1186 }
1187
1188 static void lialg_compute_chain_code(point_list *pts) {
1189         int i;
1190
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;
1198
1199                                 startpt->chaincode = dircode;
1200         }
1201 }
1202
1203
1204 static void lialg_compute_unit_chain_code(point_list *pts) {
1205         int i;
1206
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];
1213
1214                                 startpt->chaincode = dircode;
1215         }
1216 }
1217
1218
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];
1223                 int *junk;
1224                 int *curr, *next;
1225                 int i, j;
1226
1227                 /* Initialize low-pass filter parameters if necessary. */
1228                 if (lialg_lpfconst == -1)
1229                                 lialg_compute_lpf_parameters();
1230
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);
1235         curr = R[0];
1236
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;
1243                                 int a_i;
1244
1245                                 if (d_iminusone < d_i)
1246                                         d_iminusone += 8;
1247
1248                                 a_i = (d_iminusone - d_i) % 8;
1249
1250                                 /* convert to degrees, x 100 */
1251                                 curr[i] = ((12 - a_i) % 8) * 45 * 100;
1252         }
1253         curr[pts->npts - 1]     = 18000;                /* a_L-1 */
1254
1255         /* Perform a number of filtering iterations. */
1256         next = R[1];
1257         for (j = 0; j < LP_FILTER_ITERS; j++, curr = R[j], next = R[j+1]) {
1258                                 for (i = 0; i < pts->npts; i++) {
1259                                         int k;
1260
1261                                         next[i] = 0;
1262
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? */
1266                                         }
1267
1268                                         next[i] /= lialg_lpfconst;
1269                                 }
1270         }
1271
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];
1276         curr = next;
1277
1278         /* Debugging. */
1279         if (lidebug > 1) {
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]);
1285                                         fprint(2, "\n");
1286                                 }
1287         }
1288
1289         /* Do the region segmentation. */
1290         {
1291                                 int start, end;
1292                                 int currtype;
1293
1294 #define RGN_TYPE(val) (((val)==18000)?RGN_PLAIN:((val)<18000?RGN_CONCAVE:RGN_CONVEX))
1295
1296                                 start = 0;
1297                                 currtype = RGN_TYPE(curr[0]);
1298                                 regions = malloc(sizeof(region_list));
1299                                 curr_reg = regions;
1300                                 curr_reg->start = start;
1301                                 curr_reg->end = 0;
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]);
1306
1307                                         if (nexttype != currtype) {
1308                                                                 region_list *next_reg;
1309
1310                                                                 end = i - 1;
1311                                                                 curr_reg->end = end;
1312                                                                 if (lidebug > 1)
1313                                                                         fprint(2, "  (%d, %d), %d\n", start, end, currtype);
1314
1315                                                                 start = i;
1316                                                                 currtype = nexttype;
1317                                                                 next_reg = malloc(sizeof(region_list));
1318                                                                 next_reg->start = start;
1319                                                                 next_reg->end = 0;
1320                                                                 next_reg->type = nexttype;
1321                                                                 next_reg->next = nil;
1322
1323                                                                 curr_reg->next = next_reg;
1324                                                                 curr_reg = next_reg;
1325                                         }
1326                                 }
1327                                 end = i - 1;
1328                                 curr_reg->end = end;
1329                                 if (lidebug > 1)
1330                                         fprint(2, "  (%d, %d), %d\n", start, end, currtype);
1331
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;
1336
1337                                                                 for (next_reg = curr_reg->next;
1338                                                                          next_reg != nil &&
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); */
1344
1345                                                                         curr_reg->next = (next_reg->next)->next;
1346                                                                         curr_reg->end = (next_reg->next)->end;
1347
1348                                                                         free(next_reg->next);
1349                                                                         free(next_reg);
1350                                                                 }
1351                                         }
1352
1353                                 /* Add-in pseudo-extremes. */
1354                                 {
1355                                         region_list *tmp, *prev_reg;
1356
1357                                         tmp = regions;
1358                                         regions = nil;
1359                                         prev_reg = nil;
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,
1363                                                                                                                                                 curr_reg->start,
1364                                                                                                                                                 curr_reg->end);
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;
1371
1372                                                         if (lidebug)
1373                                                                 fprint(2, "%d, %d, %d\n", arclen, chordlen, atcr);
1374
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;
1380
1381                                                                 curr_reg->end = mid - 1;
1382                                                                 if (prev_reg == nil)
1383                                                                         regions = curr_reg;
1384                                                                 else
1385                                                                         prev_reg->next = curr_reg;
1386                                                                 prev_reg = curr_reg;
1387
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;
1396
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;
1405
1406                                                                 curr_reg->next = saved_next;
1407                                                                 continue;
1408                                                         }
1409                                                 }
1410
1411                                                 if (prev_reg == nil)
1412                                                         regions = curr_reg;
1413                                                 else
1414                                                         prev_reg->next = curr_reg;
1415                                                 prev_reg = curr_reg;
1416                                         }
1417                                 }
1418                 }
1419
1420         free(junk);
1421         return(regions);
1422 }
1423
1424
1425 static point_list *lialg_compute_dompts(point_list *pts, region_list *regions) {
1426         point_list *dpts;
1427         int ndpts;
1428         int *cas;
1429         int nonplain;
1430         region_list *r;
1431                 region_list *curr;
1432                 int dp;
1433                 int previx;
1434                 int currix;
1435
1436         /* Compute contour angle set. */
1437         cas = lialg_compute_contour_angle_set(pts, regions);
1438
1439         /* Dominant points include:  start_pt, end_pt, extrema_of_non_plain_regions, midpts of the preceding. */
1440         nonplain = 0;
1441         for (r = regions; r != nil; r = r->next)
1442                                 if (r->type != RGN_PLAIN)
1443                                                 nonplain++;
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) {
1449                                 free(dpts);
1450                                 return(nil);
1451         }
1452         dpts->npts = ndpts;
1453         dpts->next = nil;
1454
1455         /* Pick out dominant points. */
1456
1457                 /* Record start point. */
1458                 dp = 0;
1459                 previx = 0;
1460                 dpts->pts[dp++] = pts->pts[previx];
1461
1462                 for (curr = regions; curr != nil; curr = curr->next)
1463                         if (curr->type != RGN_PLAIN) {
1464                                                 int max_v = 0;
1465                                                 int min_v = 0x7fffffff; /* maxint */
1466                                                 int max_ix = -1;
1467                                                 int min_ix = -1;
1468                                                 int i;
1469
1470                                                 for (i = curr->start; i <= curr->end; i++) {
1471                                                         int v = cas[i];
1472                                                         if (v > max_v) { max_v = v; max_ix = i; }
1473                                                         if (v < min_v) { min_v = v; min_ix = i; }
1474                                                         if (lidebug > 1)
1475                                                                 fprint(2, "  %d\n", v);
1476                                                 }
1477
1478                                                 currix = (curr->type == RGN_CONVEX ? max_ix : min_ix);
1479
1480                                                 /* Record midpoint. */
1481                                                 dpts->pts[dp++] = pts->pts[previx + (currix - previx) / 2];
1482
1483                                                 /* Record extreme point. */
1484                                                 dpts->pts[dp++] = pts->pts[currix];
1485
1486                                                 previx = currix;
1487                         }
1488
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];
1493
1494         /* Compute chain-code. */
1495         lialg_compute_chain_code(dpts);
1496
1497         free(cas);
1498         return(dpts);
1499 }
1500
1501
1502 static int *lialg_compute_contour_angle_set(point_list *pts,
1503                                                                                            region_list *regions) {
1504                 int *V;
1505                 region_list *curr_reg;
1506                 int i;
1507
1508                 V = malloc(pts->npts*sizeof(int));
1509
1510                 V[0] = 18000;
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) {
1514                                                                 V[i] = 18000;
1515                                         } else {
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);
1522                                         }
1523                                 }
1524         }
1525         V[pts->npts - 1] = 18000;
1526
1527         return(V);
1528 }
1529
1530
1531 /*
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.
1536  *
1537  */
1538 static void lialg_score_stroke(point_list *input_dompts, point_list *curr_dompts, int *sim, int *dist) {
1539         *sim = MIN_SIM;
1540         *dist = MAX_DIST;
1541
1542         *sim = lialg_compute_similarity(input_dompts, curr_dompts);
1543         if (*sim < SIM_THLD) goto done;
1544
1545         *dist = lialg_compute_distance(input_dompts, curr_dompts);
1546
1547 done:
1548         if (lidebug)
1549                 fprint(2, "%d, %d\n", *sim, *dist);
1550 }
1551
1552
1553 static int lialg_compute_similarity(point_list *input_dompts, point_list *curr_dompts) {
1554         int sim;
1555         point_list *A, *B;
1556         int N, M;
1557         int **G;
1558         int *junk;
1559         int i, j;
1560
1561         /* A is the     longer sequence, length N. */
1562         /* B is the shorter sequence, length M. */
1563                 if (input_dompts->npts >= curr_dompts->npts) {
1564                                 A = input_dompts;
1565                                 N = input_dompts->npts;
1566                                 B = curr_dompts;
1567                                 M = curr_dompts->npts;
1568         } else {
1569                                 A = curr_dompts;
1570                                 N = curr_dompts->npts;
1571                                 B = input_dompts;
1572                                 M = input_dompts->npts;
1573                 }
1574
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));
1583
1584                 for (i = 1; i < M; i++) {
1585                         int bval = B->pts[i-1].chaincode;
1586
1587                         /* Source column. */
1588                         G[i][0] = 0;
1589
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;
1594
1595                                 G[i][j] = (diff == 0)
1596                                   ? 10
1597                                   : (diff == 1)
1598                                   ? 6
1599                                   : 0;
1600                         }
1601
1602                         /* Sink column. */
1603                         G[i][N] = 0;
1604                 }
1605
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];
1613
1614                                 if (i < (M - 1)) {
1615                                         int tmp = G[i + 1][j + 1];
1616                                         if (tmp > max) max = tmp;
1617                                 }
1618
1619                                 G[i][j] += max;
1620                 }
1621
1622                 sim = (10 * G[1][0] + (N - 1) / 2) / (N - 1);
1623
1624                 if (G != nil)
1625                                 free(G);
1626                 if (junk != nil)
1627                                 free(junk);
1628                 return(sim);
1629 }
1630
1631
1632 static int lialg_compute_distance(point_list *input_dompts,
1633                                                                    point_list *curr_dompts) {
1634                 int dist;
1635                 point_list *A, *B;
1636                 int N, M;
1637                 int **C;
1638                 int *junk;
1639                 int *BE;
1640                 int *TE;
1641                 int i, j;
1642
1643                 /* A is the     longer sequence, length N. */
1644                 /* B is the shorter sequence, length M. */
1645                 if (input_dompts->npts >= curr_dompts->npts) {
1646                                 A = input_dompts;
1647                                 N = input_dompts->npts;
1648                                 B = curr_dompts;
1649                                 M = curr_dompts->npts;
1650                 }
1651                 else {
1652                                 A = curr_dompts;
1653                                 N = curr_dompts->npts;
1654                                 B = input_dompts;
1655                                 M = input_dompts->npts;
1656                 }
1657
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));
1662
1663                 for (j = 1; j <= N; j++) {
1664                                 int bot, top;
1665
1666                                 bot = j + (M - DP_BAND);
1667                                 if (bot > M) bot = M;
1668                                 BE[j] = bot;
1669
1670                                 top = j - (N - DP_BAND);
1671                                 if (top < 1) top = 1;
1672                                 TE[j] = top;
1673                 }
1674
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));
1684
1685                 for (i = 1; i <= M; i++) {
1686                                 int bx = B->pts[i-1].x;
1687                                 int by = B->pts[i-1].y;
1688
1689                                 for (j = 1; j <= N; j++) {
1690                                                 int ax = A->pts[j-1].x;
1691                                                 int ay = A->pts[j-1].y;
1692                                                 int dx = bx - ax;
1693                                                 int dy = by - ay;
1694                                                 int dist = isqrt(10000 * (dx * dx + dy * dy));
1695
1696                                                 C[i][j] = dist;
1697                                 }
1698                 }
1699
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--) {
1705                                                 int min = MAX_DIST;
1706
1707                                                 if (i > BE[j] || i < TE[j] || (j == N && i == M))
1708                                                                 continue;
1709
1710                                                 if (j < N) {
1711                                                                 if (i >= TE[j+1]) {
1712                                                                                 int tmp = C[i][j+1];
1713                                                                                 if (tmp < min)
1714                                                                                                 min = tmp;
1715                                                                 }
1716
1717                                                                 if (i < M) {
1718                                                                                 int tmp = C[i+1][j+1];
1719                                                                                 if (tmp < min)
1720                                                                                                 min = tmp;
1721                                                                 }
1722                                                 }
1723
1724                                                 if (i < BE[j]) {
1725                                                                 int tmp = C[i+1][j];
1726                                                                 if (tmp < min) min = tmp;
1727                                                 }
1728
1729                                                 C[i][j] += min;
1730                                 }
1731
1732                 dist = (C[1][1] + N / 2) / N;
1733
1734                 if (C != nil) free(C);
1735                 if (junk != nil) free(junk);
1736                 if (BE != nil) free(BE);
1737                 if (TE != nil) free(TE);
1738                 return(dist);
1739 }
1740
1741
1742 /*************************************************************
1743
1744   Digest-processing routines
1745
1746  *************************************************************/
1747
1748 static int lialg_read_classifier_digest(rClassifier *rec) {
1749         int nclasses;
1750         FILE *fp;
1751
1752         /* Try to open the corresponding digest file. */
1753         {
1754                                 char *clx_path;
1755                                 char *dot;
1756
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);
1761
1762                                 /* Truncate the path after the last dot. */
1763                                 dot = strrchr(clx_path, '.');
1764                                 if (dot == nil) { free(clx_path); return(-1); }
1765                                 *(dot + 1) = 0;
1766
1767                                 /* Append the classifier-digest extension. */
1768                                 strcat(clx_path, "clx");
1769
1770                                 fp = fopen(clx_path, "r");
1771                                 if (fp == nil) {
1772                                                 free(clx_path);
1773                                                 return(-1);
1774                                 }
1775
1776                                 free(clx_path);
1777         }
1778
1779         /* Read-in the name and dominant points for each class. */
1780         for (nclasses = 0; !feof(fp); nclasses++) {
1781                 point_list *dpts = nil;
1782                 char class[BUFSIZ];
1783                 int npts;
1784                 int j;
1785
1786                 if (fscanf(fp, "%s %d", class, &npts) != 2) {
1787                         if (feof(fp)) break;
1788
1789                         goto failed;
1790                 }
1791                 rec->cnames[nclasses] = strdup(class);
1792
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;
1798                 dpts->npts = npts;
1799                 dpts->next = nil;
1800
1801                 /* Read in each point. */
1802                 for (j = 0; j < npts; j++) {
1803                         int x, y;
1804
1805                         if (fscanf(fp, "%d %d", &x, &y) != 2) goto failed;
1806                         dpts->pts[j].x = x;
1807                         dpts->pts[j].y = y;
1808                 }
1809
1810                 /* Compute the chain-code. */
1811                 lialg_compute_chain_code(dpts);
1812
1813                 /* Store the list in the rec data structure. */
1814                 rec->dompts[nclasses] = dpts;
1815
1816                 continue;
1817
1818 failed:
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;
1824                         }
1825                         if (rec->dompts[nclasses] != nil) {
1826                                 delete_examples(rec->dompts[nclasses]);
1827                                 rec->dompts[nclasses] = nil;
1828                         }
1829                 }
1830                 if (dpts != nil)
1831                         delete_examples(dpts);
1832                 fclose(fp);
1833                 return(-1);
1834         }
1835
1836         fclose(fp);
1837         return(0);
1838 }
1839
1840
1841 /*************************************************************
1842
1843   Canonicalization routines
1844
1845  *************************************************************/
1846
1847 static int lialg_canonicalize_examples(rClassifier *rec) {
1848         int i;
1849         int nclasses;
1850
1851         if (lidebug) {
1852                 fprint(2, "lialg_canonicalize_examples working on %s\n",
1853                                 rec->file_name);
1854         }
1855         /* Initialize canonical-example arrays. */
1856         for (i = 0; i < MAXSCLASSES; i++) {
1857                 rec->canonex[i] = nil;
1858         }
1859
1860         /* Figure out number of classes. */
1861         for (nclasses = 0;
1862                   nclasses < MAXSCLASSES && rec->cnames[nclasses] != nil;
1863                   nclasses++)
1864                 ;
1865
1866         /* Canonicalize the examples for each class. */
1867         for (i = 0; i < nclasses; i++) {
1868                 int j, k;
1869                 int nex;
1870                 point_list *pts, *tmp, *avg;
1871                 int maxxrange, maxyrange;
1872                 int minx, miny, maxx, maxy;
1873                 int avgxrange, avgyrange, avgxoff, avgyoff, avgscale;
1874
1875
1876                 if (lidebug) {
1877                         fprint(2, "lialg_canonicalize_examples working on class %s\n",
1878                                         rec->cnames[i]);
1879                 }
1880                 /* Make a copy of the examples. */
1881                 pts = nil;
1882                 tmp = rec->ex[i];
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);
1886                                 return(-1);
1887                         }
1888                 }
1889
1890                 /* Canonicalize each example, and derive the max x and y ranges. */
1891                 maxxrange = 0;
1892                 maxyrange = 0;
1893                 for (j = 0, tmp = pts; j < nex; j++, tmp = tmp->next) {
1894                         if (lialg_canonicalize_example_stroke(tmp) != 0) {
1895                 if (lidebug) {
1896                                         fprint(2, "lialg_canonicalize_example_stroke returned error\n");
1897                                 }
1898                                 return(-1);
1899                         }
1900
1901                         if (tmp->xrange > maxxrange) maxxrange = tmp->xrange;
1902                         if (tmp->yrange > maxyrange) maxyrange = tmp->yrange;
1903                 }
1904
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;
1910                 }
1911                 else {
1912                         maxxrange = (maxxrange * CANONICAL_Y + maxyrange / 2) / maxyrange;
1913                         maxyrange = CANONICAL_Y;
1914                 }
1915
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);
1922                                 return(-1);
1923                         }
1924                 }
1925
1926                 /* Average the examples; leave average in first example. */
1927                 avg = pts;                              /* careful aliasing!! */
1928                 for (k = 0; k < NCANONICAL; k++) {
1929                         int xsum = 0;
1930                         int ysum = 0;
1931
1932                         for (j = 0, tmp = pts; j < nex; j++, tmp = tmp->next) {
1933                                                 xsum += tmp->pts[k].x;
1934                                                 ysum += tmp->pts[k].y;
1935                         }
1936
1937                         avg->pts[k].x = (xsum + j / 2) / j;
1938                         avg->pts[k].y = (ysum + j / 2) / j;
1939                 }
1940
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);
1951                         return(-1);
1952                 }
1953
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);
1962                         return(-1);
1963                 }
1964
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);
1968                         return(-1);
1969                 }
1970                 (rec->canonex[i])->xrange = maxx - minx;
1971                 (rec->canonex[i])->yrange = maxy - miny;
1972
1973                 if (lidebug) {
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);
1977                         }
1978                 }
1979
1980                 /* Compute dominant points of canonical representation. */
1981                 rec->dompts[i] = lialg_compute_dominant_points(avg);
1982
1983                 /* Clean up. */
1984                 delete_examples(pts);
1985         }
1986
1987         /* Sanity check. */
1988         for (i = 0; i < nclasses; i++) {
1989                 char *best_name = lialg_recognize_stroke(rec, rec->canonex[i]);
1990
1991                 if (best_name != rec->cnames[i])
1992                         fprint(2, "%s, best = %s\n", rec->cnames[i], best_name);
1993         }
1994
1995         return(0);
1996 }
1997
1998
1999 static int lialg_canonicalize_example_stroke(point_list *points) {
2000         int minx, miny, maxx, maxy, xrange, yrange, scale;
2001
2002         /* Filter out points that are too close. */
2003         if (lialg_filter_points(points) != 0) return(-1);
2004
2005         /* Must be at least two points! */
2006         if (points->npts < 2) {
2007                 if (lidebug) {
2008                         fprint(2, "lialg_canonicalize_example_stroke: npts=%d\n",
2009                                         points->npts);
2010                 }
2011                 return(-1);
2012         }
2013
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) {
2023                 if (lidebug) {
2024                         fprint(2, "lialg_translate_points (minx=%d,miny=%d,scale=%d) returned error\n", minx, miny, scale);
2025                 }
2026                 return(-1);
2027         }
2028
2029         /* Compute an equivalent stroke with equi-distant points. */
2030         if (lialg_compute_equipoints(points) != 0) return(-1);
2031
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) {
2035                 if (lidebug) {
2036                         fprint(2, "lialg_translate_points (minx=%d,miny=%d) returned error\n", minx, miny);
2037                 }
2038                 return(-1);
2039         }
2040
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;
2046
2047         if (lidebug) {
2048                 int i;
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);
2052                 fflush(stderr);
2053         }
2054
2055         return(0);
2056 }
2057
2058
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);
2064         int i;
2065         int dist_since_last_eqpt;
2066         int remaining_seglen;
2067         int dist_to_next_eqpt;
2068
2069         if (equipoints == nil) {
2070                 fprint(2, "can't allocate memory in lialg_compute_equipoints");
2071                 return(-1);
2072         }
2073
2074         if (lidebug) {
2075                 fprint(2, "compute_equipoints:  npts = %d, pathlen = %d, equidist = %d\n",
2076                                 points->npts, pathlen, equidist);
2077                 fflush(stderr);
2078         }
2079
2080         /* First original point is an equipoint. */
2081         equipoints[0] = points->pts[0];
2082         nequipoints++;
2083         dist_since_last_eqpt = 0;
2084
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;
2092
2093                 while (remaining_seglen >= dist_to_next_eqpt) {
2094                         if (dx1 == 0) {
2095                                 /* x-coordinate stays the same */
2096                                 if (dy1 >= 0)
2097                                         endy += dist_to_next_eqpt;
2098                                 else
2099                                         endy -= dist_to_next_eqpt;
2100                         }
2101                         else {
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;
2106
2107                                 if (dy < 0) dy = -dy;
2108                                 if (dx1 >= 0)
2109                                         endx += dx;
2110                                 else
2111                                         endx -= dx;
2112                                 if (dy1 >= 0)
2113                                         endy += dy;
2114                                 else
2115                                         endy -= dy;
2116                         }
2117
2118                         equipoints[nequipoints].x = (endx + 50) / 100;
2119                         equipoints[nequipoints].y = (endy + 50) / 100;
2120                         nequipoints++;
2121 /*          assert(nequipoints <= NCANONICAL);*/
2122                         dist_since_last_eqpt = 0;
2123                         remaining_seglen -= dist_to_next_eqpt;
2124                         dist_to_next_eqpt = equidist;
2125                 }
2126
2127                 dist_since_last_eqpt += remaining_seglen;
2128         }
2129
2130         /* Take care of last equipoint. */
2131         if (nequipoints == NCANONICAL) {
2132                 /* Good. */
2133         } else if (nequipoints == (NCANONICAL - 1)) {
2134                 /* Make last original point the last equipoint. */
2135                 equipoints[nequipoints] = points->pts[points->npts - 1];
2136         } else {
2137           if (lidebug) {
2138                 fprint(2,"lialg_compute_equipoints: nequipoints = %d\n",
2139                                 nequipoints);
2140           }
2141 /*      assert(false);*/
2142                 return(-1);
2143         }
2144
2145         points->npts = NCANONICAL;
2146         delete_pen_point_array(points->pts);
2147         points->pts = equipoints;
2148         return(0);
2149 }
2150
2151
2152 /*************************************************************
2153
2154   Utility routines
2155
2156  *************************************************************/
2157
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));
2161 }
2162
2163
2164 /* Result is x 100. */
2165 static int lialg_compute_pathlen_subset(point_list *points,
2166                                                                                    int start, int end) {
2167         int pathlen;
2168         int i;
2169
2170         pathlen = 0;
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));
2175                 pathlen += dist;
2176         }
2177
2178         return(pathlen);
2179 }
2180
2181
2182 /* Note that this does NOT update points->xrange and points->yrange! */
2183 static int lialg_filter_points(point_list *points) {
2184         int filtered_npts;
2185         pen_point *filtered_pts = mallocz(points->npts*sizeof(pen_point), 1);
2186         int i;
2187
2188         if (filtered_pts == nil) {
2189                 fprint(2, "can't allocate memory in lialg_filter_points");
2190                 return(-1);
2191         }
2192
2193         filtered_pts[0] = points->pts[0];
2194         filtered_npts = 1;
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;
2200
2201                 if (magsq >= DIST_SQ_THRESHOLD) {
2202                         filtered_pts[filtered_npts] = points->pts[i];
2203                         filtered_npts++;
2204                 }
2205         }
2206
2207         points->npts = filtered_npts;
2208         delete_pen_point_array(points->pts);
2209         points->pts = filtered_pts;
2210         return(0);
2211 }
2212
2213
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,
2217                                                                    int minx, int miny,
2218                                                                    int scalex, int scaley) {
2219         int i;
2220
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;
2224         }
2225
2226         return(0);
2227 }
2228
2229
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;
2234         int i;
2235
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;
2244         }
2245
2246         *pminx = minx;
2247         *pminy = miny;
2248         *pmaxx = maxx;
2249         *pmaxy = maxy;
2250 }
2251
2252
2253 int wtvals[] = {100, 104, 117, 143, 189, 271, 422};
2254
2255 static void lialg_compute_lpf_parameters(void) {
2256         int i;
2257
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);
2262                                 int wt = wtvals[i];
2263                                 lialg_lpfwts[LP_FILTER_WIDTH - i] = wt;
2264                                 lialg_lpfwts[LP_FILTER_WIDTH + i] = wt;
2265                 }
2266                 lialg_lpfconst = 0;
2267                 for (i = 0; i < (2 * LP_FILTER_WIDTH + 1); i++) {
2268                                 lialg_lpfconst += lialg_lpfwts[i];
2269                 }
2270 }
2271
2272
2273 /* Code from Joseph Hall (jnhall@sat.mot.com). */
2274 static int isqrt(int n) {
2275                 register int i;
2276                 register long k0, k1, nn;
2277
2278                 for (nn = i = n, k0 = 2; i > 0; i >>= 2, k0 <<= 1)
2279                                 ;
2280                 nn <<= 2;
2281         for (;;) {
2282                                 k1 = (nn / k0 + k0) >> 1;
2283                                 if (((k0 ^ k1) & ~1) == 0)
2284                                 break;
2285                                 k0 = k1;
2286                 }
2287                 return (int) ((k1 + 1) >> 1);
2288 }
2289
2290
2291 /* Helper routines from Mark Hayter. */
2292 static int likeatan(int tantop, int tanbot) {
2293                 int t;
2294                 /* Use tan(theta)=top/bot --> order for t */
2295                 /* t in range 0..0x40000 */
2296
2297                 if ((tantop == 0) && (tanbot == 0))
2298                                 t = 0;
2299         else
2300                 {
2301                                 t = (tantop << 16) / (abs(tantop) + abs(tanbot));
2302                                 if (tanbot < 0)
2303                                                 t = 0x20000 - t;
2304                                 else
2305                                                 if (tantop < 0) t = 0x40000 + t;
2306                 }
2307                 return t;
2308 }
2309
2310
2311 static int quadr(int t) {
2312         return (8 - (((t + 0x4000) >> 15) & 7)) & 7;
2313 }