packages/services/gfx/mw/v2_0/src/demos/nxscribble/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
  20 #include <sys/types.h>
  21 #include <stdio.h>
  22 #include <string.h>
  23 #ifndef ELX
  24 #include <stdlib.h>
  25 #endif
  26 #include <math.h>
  27 #include <locale.h>
  28 #include <hre_internal.h>
  29 #include <setjmp.h>
  30 #include "util.h"
  31 #include "matrix.h"
  32 #include "sc.h"
  33 #include "li_recognizer.h"
  34 #include "li_recognizer_internal.h"
  35
  36
  37 int lidebug = 0;
  38
  39 /*LI Magic Number.*/
  40
  41 #define LI_MAGIC 0xACCBADDD
  42
  43 #define CHECK_LI_MAGIC(_a) \
  44   ((_a) != NULL && ((li_recognizer*)(_a))->li_magic == LI_MAGIC)
  45
  46
  47 static void lialg_initialize(rClassifier *);
  48 static int lialg_read_classifier_digest(rClassifier *);
  49 static int lialg_canonicalize_examples(rClassifier *);
  50 static char *lialg_recognize_stroke(rClassifier *, point_list *);
  51
  52
  53 char* li_err_msg = NULL;
  54 char _zdebug_flag[128];
  55
  56 #ifndef __ECOS
  57 // This is standard - defined in <stdlib.h>
  58 #define bcopy(s1,s2,n) memcpy(s2,s1,n)
  59 #endif
  60
  61 #if 0 /* was #ifdef mips*/
  62 char *strdup(char* from) {
  63    char* to;
  64    int len = strlen(from) + 1;
  65
  66    /*   to = (char *) safe_malloc( len * sizeof(char) );*/
  67    to = allocate(len, char);
  68    memcpy(to, from, len);
  69    return to;
  70 }
  71 #endif
  72
  73
  74 /*Freeing classifier*/
  75
  76 static void
  77 free_rClassifier(rClassifier* rc);
  78
  79 /*
  80  * Point List Support
  81 */
  82
  83 static point_list*
  84 add_example(point_list* l,int npts,pen_point* pts)
  85 {
  86     pen_point* lpts = make_pen_point_array(npts);
  87     /*    point_list* p = (point_list*)safe_malloc(sizeof(point_list));*/
  88     point_list *p = allocate(1, point_list);
  89
  90     p->npts = npts;
  91     p->pts = lpts;
  92     p->next = l;       /*Order doesn't matter, so we stick on end.*/
  93
  94     /*Copy points.*/
  95
  96     bcopy(pts,lpts,npts * sizeof(pen_point));
  97
  98     return(p);
  99 }
 100
 101
 102 static void
 103 delete_examples(point_list* l)
 104 {
 105     point_list* p;
 106
 107     for( ; l != NULL; l = p ) {
 108         p = l->next;
 109         free(l->pts);
 110         free(l);
 111     }
 112 }
 113
 114 /*
 115  * Local functions
 116  */
 117
 118 /*
 119  * recognize_internal-Form Vector, use Classifier to classify, return char.
 120  */
 121
 122 static char*
 123   recognize_internal(rClassifier* rec,pen_stroke* str,int* rconf)
 124 {
 125     char *res = NULL;
 126     point_list *stroke = NULL;
 127
 128     stroke = add_example(NULL, str->ps_npts, str->ps_pts);
 129     if (stroke == NULL) return(NULL);
 130
 131     res = lialg_recognize_stroke(rec, stroke);
 132
 133     delete_examples(stroke);
 134     return(res);
 135 }
 136
 137 /*
 138  * file_path-Construct pathname, check for proper extension.
 139  */
 140
 141 static int
 142   file_path(char* dir,char* filename,char* pathname)
 143 {
 144     char* dot;
 145
 146     /*Check for proper extension on file name.*/
 147
 148     dot = strrchr(filename,'.');
 149
 150     if( dot == NULL ) {
 151         return(-1);
 152     }
 153
 154     /*Determine whether a gesture or character classifier.*/
 155
 156     if( strcmp(dot,LI_CLASSIFIER_EXTENSION) != 0 ) {
 157         return(-1);
 158     }
 159
 160     /*Concatenate directory and filename into pathname.*/
 161
 162     strcpy(pathname,dir);
 163     strcat(pathname,"/");
 164     strcat(pathname,filename);
 165
 166     return(0);
 167 }
 168
 169 /*read_classifier_points-Read points so classifier can be extended.*/
 170
 171 static int
 172 read_classifier_points(FILE* fd,int nclss,point_list** ex,char** cnames)
 173 {
 174     int i,j,k;
 175     char buf[BUFSIZ];
 176     int nex = 0;
 177     char* names[MAXSCLASSES];
 178     point_list* examples[MAXSCLASSES];
 179     pen_point* pts;
 180     int npts;
 181
 182     /*Initialize*/
 183
 184     for( i = 0; i < MAXSCLASSES; i++ ) {
 185         names[i] = NULL;
 186         examples[i] = NULL;
 187     }
 188
 189     /*Go thru classes.*/
 190
 191 /* ari */
 192 /* fprintf(stderr, "Classes: [ "); */
 193
 194     for( k = 0; k < nclss; k++ ) {
 195
 196         /*Read class name and number of examples.*/
 197
 198         if( fscanf(fd,"%d %s",&nex,buf) != 2 ) {
 199             printf("%s *FAILED* - line: %d\n", __FUNCTION__, __LINE__);
 200             goto unallocate;
 201         }
 202
 203         /*Save class name*/
 204
 205         names[k] = strdup(buf);
 206 /* ari */
 207 /* fprintf(stderr, "%s ", buf); */
 208
 209         /*Read examples.*/
 210
 211         for( i = 0; i < nex; i++ ) {
 212
 213             /*Read number of points.*/
 214
 215             if( fscanf(fd,"%d",&npts) != 1 ) {
 216                 printf("%s *FAILED* - line: %d\n", __FUNCTION__, __LINE__);
 217                 goto unallocate; /*Boy would I like exceptions!*/
 218             }
 219
 220             /*Allocate array for points.*/
 221
 222             if( (pts = make_pen_point_array(npts)) == NULL ) {
 223                 printf("%s *FAILED* - line: %d\n", __FUNCTION__, __LINE__);
 224                 goto unallocate;
 225             }
 226
 227             /*Read in points.*/
 228
 229             for( j = 0; j < npts; j++ ) {
 230                 int x,y;
 231                 int jj;
 232                 if( fscanf(fd,"%d %d",&x,&y) != 2 ) {
 233                     delete_pen_point_array(pts);
 234                     printf("%s *FAILED* - line: %d\n", __FUNCTION__, __LINE__);
 235                     printf("class = %d/%d/%s, ex = %d/%d, pt: %d/%d\n",
 236                                 k, nclss, names[k], i, nex, j, npts);
 237                     for (jj = 0;  jj < j;  jj++) {
 238                         printf("pts[%d] = %d/%d\n", jj, pts[jj].x, pts[jj].y);
 239                     }
 240                     goto unallocate;
 241                 }
 242                 pts[j].x = x;
 243                 pts[j].y = y;
 244             }
 245
 246             /*Add example*/
 247
 248             if( (examples[k] = add_example(examples[k],npts,pts)) == NULL ) {
 249                 delete_pen_point_array(pts);
 250                 printf("%s *FAILED* - line: %d\n", __FUNCTION__, __LINE__);
 251                 goto unallocate;
 252             }
 253
 254             delete_pen_point_array(pts);
 255
 256           }
 257       }
 258
 259 /* ari -- end of list of classes */
 260 /* fprintf(stderr, "]\n"); */
 261
 262     /*Transfer to recognizer.*/
 263
 264     bcopy(examples,ex,sizeof(examples));
 265     bcopy(names,cnames,sizeof(names));
 266
 267     return(0);
 268
 269     /*Error. Deallocate memory and return.*/
 270
 271   unallocate:
 272
 273     for( ; k >= 0; k-- ) {
 274         delete_examples(examples[k]);
 275         free(names[k]);
 276     }
 277
 278     error("Error in reading example points from classifier file");
 279     return(-1);
 280 }
 281
 282 /*read_classifier-Read a classifier file.*/
 283
 284 static int read_classifier(FILE* fd,rClassifier* rc)
 285 {
 286     sClassifier sc;
 287
 288     li_err_msg = NULL;
 289
 290     /*Read in classifier file.*/
 291
 292     if( (sc = sRead(fd)) == NULL ) {
 293         return(-1);
 294     }
 295
 296     /*Read in the example points, so classifier can be extended.*/
 297
 298     if( read_classifier_points(fd,sc->nclasses,rc->ex,rc->cnames) != 0 ) {
 299         sFreeClassifier(sc);
 300         return(-1);
 301     }
 302
 303     /*Transfer sClassifier to the rClassifier*/
 304
 305     rc->sc = sc;
 306
 307     return(0);
 308 }
 309
 310 /*
 311  * Extension Functions
 312 */
 313
 314 /* getClasses and clearState are by Ari */
 315
 316 static int
 317 recognizer_getClasses (recognizer r, char ***list, int *nc)
 318 {
 319     int i, nclasses;
 320     li_recognizer* rec;
 321     sClassifier sc;
 322     char **ret;
 323
 324     rec = (li_recognizer*)r->recognizer_specific;
 325
 326     /*Check for LI recognizer.*/
 327
 328     if( !CHECK_LI_MAGIC(rec) ) {
 329         li_err_msg = "Not a LI recognizer";
 330         return(-1);
 331     }
 332
 333     sc = rec->li_rc.sc;
 334     *nc = nclasses = sc->nclasses;
 335     /*    ret = (char **) safe_malloc (nclasses * sizeof(char*));*/
 336     ret = allocate(nclasses, char*);
 337
 338     for (i = 0; i < nclasses; i++) {
 339       ret[i] = rec->li_rc.cnames[i];   /* only the 1st char of the cname */
 340     }
 341     *list = ret;
 342     return 0;
 343 }
 344
 345 static int
 346 recognizer_clearState (recognizer r)
 347 {
 348   /*This operation isn't supported by the LI recognizer.*/
 349
 350   li_err_msg = "Clearing state is not supported by the LI recognizer";
 351
 352   return(-1);
 353 }
 354
 355 static bool isa_li(recognizer r)
 356 { return(CHECK_LI_MAGIC(r)); }
 357
 358 static int
 359 recognizer_train(recognizer r,rc* rec_xt,u_int nstrokes,
 360                     pen_stroke* strokes,rec_element* re,
 361                     bool replace_p)
 362 {
 363   /*This operation isn't supported by the LI recognizer.*/
 364
 365   li_err_msg = "Training is not supported by the LI recognizer";
 366
 367   return(-1);
 368 }
 369
 370 int
 371 li_recognizer_get_example (recognizer   r,
 372                            int          class,
 373                            int          instance,
 374                            char         **name,
 375                            pen_point    **points,
 376                            int          *npts)
 377 {
 378     li_recognizer   *rec = (li_recognizer*)r->recognizer_specific;
 379     sClassifier     sc = rec->li_rc.sc;
 380     point_list      *pl;
 381
 382     if( !CHECK_LI_MAGIC(rec) ) {
 383         li_err_msg = "Not a LI recognizer";
 384         return(-1);
 385     }
 386     if (class > sc->nclasses)
 387         return -1;
 388     pl = rec->li_rc.canonex[class];
 389     while (instance && pl)
 390     {
 391         pl = pl->next;
 392         instance--;
 393     }
 394     if (!pl)
 395         return -1;
 396     *name = rec->li_rc.cnames[class];
 397     *points = pl->pts;
 398     *npts = pl->npts;
 399     return 0;
 400 }
 401
 402 /*
 403  * API Functions
 404  */
 405
 406
 407 /*li_recognizer_load-Load a classifier file.*/
 408
 409 static int li_recognizer_load(recognizer r,char* dir,char* filename)
 410 {
 411     FILE *fd;
 412     char* pathname;
 413     li_recognizer* rec;
 414     rClassifier* rc;
 415
 416     rec = (li_recognizer*)r->recognizer_specific;
 417
 418     /*Make sure recognizer's OK*/
 419
 420     if( !CHECK_LI_MAGIC(rec) ) {
 421         li_err_msg = "Not a LI recognizer";
 422         return(-1);
 423     }
 424
 425     rc = &(rec->li_rc);
 426
 427     /*Check parameters.*/
 428
 429     if( filename == NULL ) {
 430         li_err_msg = "Invalid parameters";
 431         return(-1);
 432     }
 433
 434     /*We let the directory be null.*/
 435
 436     if( dir == NULL || (int)strlen(dir) <= 0 ) {
 437         dir = ".";
 438     }
 439
 440     /*Make full pathname and check filename*/
 441
 442     /* pathname = (char*)safe_malloc(strlen(dir) + strlen(filename) + 2)); */
 443
 444     pathname = allocate(strlen(dir) + strlen(filename) + 2, char);
 445     if( file_path(dir,filename,pathname) == -1 ) {
 446         free(pathname);
 447         li_err_msg = "Not a LI recognizer classifier file";
 448         return(-1);
 449     }
 450
 451     /* Try to short-circuit the full classifier-file processing. */
 452     rc->file_name = pathname;
 453     if (lialg_read_classifier_digest(rc) == 0)
 454         return(0);
 455     rc->file_name = NULL;
 456
 457     /*Open the file*/
 458
 459     if( (fd = fopen(pathname,"r")) == NULL ) {
 460         free(pathname);
 461         li_err_msg = "Can't open classifier file";
 462 /* ari */
 463         /* fprintf(stderr, "Trying to open %s.\n", pathname); */
 464         return(-1);
 465
 466     }
 467
 468     /*If rClassifier is OK, then delete it first.*/
 469
 470     if( rc->file_name != NULL ) {
 471       free_rClassifier(rc);
 472     }
 473
 474     /*Read classifier.*/
 475
 476     if( read_classifier(fd,rc) < 0 ) {
 477         free(pathname);
 478         return(-1);
 479     }
 480
 481     /*Close file.*/
 482
 483     fclose(fd);
 484
 485     /*Add classifier name.*/
 486
 487     rc->file_name = pathname;
 488
 489     /* Canonicalize examples. */
 490     if (lialg_canonicalize_examples(rc) != 0) {
 491         free(pathname);
 492         rc->file_name = NULL;
 493         return(-1);
 494     }
 495
 496     return(0);
 497 }
 498
 499 /*li_recognizer_save-Save a classifier file.*/
 500
 501 static int li_recognizer_save(recognizer r,char* dir,char* filename)
 502 {
 503   /*This operation isn't supported by the LI recognizer.*/
 504
 505   li_err_msg = "Saving is not supported by the LI recognizer";
 506
 507   return(-1);
 508 }
 509
 510 static wordset
 511 li_recognizer_load_dictionary(recognizer rec,char* directory,char* name)
 512 {
 513   /*This operation isn't supported by the LI recognizer.*/
 514
 515   li_err_msg = "Dictionaries are not supported by the LI recognizer";
 516
 517   return(NULL);
 518
 519 }
 520
 521 static int
 522 li_recognizer_save_dictionary(recognizer rec,
 523                                char* directory,
 524                                char* name,
 525                                wordset dict)
 526 {
 527   /*This operation isn't supported by the LI recognizer.*/
 528
 529   li_err_msg = "Dictionaries are not supported by the LI recognizer";
 530
 531   return(-1);
 532
 533 }
 534
 535 static int
 536 li_recognizer_free_dictionary(recognizer rec,wordset dict)
 537 {
 538   /*This operation isn't supported by the LI recognizer.*/
 539
 540   li_err_msg = "Dictionaries are not supported by the LI recognizer";
 541
 542   return(-1);
 543
 544 }
 545
 546 static int
 547 li_recognizer_add_to_dictionary(recognizer rec,letterset* word,wordset dict)
 548 {
 549   /*This operation isn't supported by the LI recognizer.*/
 550
 551   li_err_msg = "Dictionaries are not supported by the LI recognizer";
 552
 553   return(-1);
 554
 555 }
 556
 557 static int
 558 li_recognizer_delete_from_dictionary(recognizer rec,
 559                                       letterset* word,
 560                                       wordset dict)
 561 {
 562   /*This operation isn't supported by the LI recognizer.*/
 563
 564   li_err_msg = "Dictionaries are not supported by the LI recognizer";
 565
 566   return(-1);
 567
 568 }
 569
 570 static char*
 571 li_recognizer_error(recognizer rec)
 572 {
 573     char* ret = li_err_msg;
 574
 575     /*Check for LI recognizer.*/
 576
 577     if( !CHECK_LI_MAGIC(rec->recognizer_specific) ) {
 578         li_err_msg = "Not a LI recognizer";
 579         return(NULL);
 580     }
 581
 582     li_err_msg = NULL;
 583
 584     return(ret);
 585 }
 586
 587 static int
 588 li_recognizer_clear(recognizer r,bool delete_points_p)
 589 {
 590     li_recognizer* rec;
 591
 592     rec = (li_recognizer*)r->recognizer_specific;
 593
 594     /*Check for LI recognizer.*/
 595
 596     if( !CHECK_LI_MAGIC(rec) ) {
 597         li_err_msg = "Not a LI recognizer";
 598         return(0);
 599     }
 600
 601     return(0);
 602 }
 603
 604 static int
 605 li_recognizer_set_context(recognizer r,rc* rec_xt)
 606 {
 607   /*This operation isn't supported by the LI recognizer.*/
 608
 609   li_err_msg = "Contexts are not supported by the LI recognizer";
 610
 611   return(-1);
 612 }
 613
 614 static rc*
 615 li_recognizer_get_context(recognizer r)
 616 {
 617   /*This operation isn't supported by the LI recognizer.*/
 618
 619   li_err_msg = "Contexts are not supported by the LI recognizer";
 620
 621   return(NULL);
 622 }
 623
 624 static int
 625 li_recognizer_get_buffer(recognizer r, u_int* nstrokes,pen_stroke** strokes)
 626 {
 627   /*This operation isn't supported by the LI recognizer.*/
 628
 629   li_err_msg = "Buffer get/set are not supported by the LI recognizer";
 630
 631   return(-1);
 632 }
 633
 634 static int
 635 li_recognizer_set_buffer(recognizer r,u_int nstrokes,pen_stroke* strokes)
 636 {
 637   /*This operation isn't supported by the LI recognizer.*/
 638
 639   li_err_msg = "Buffer get/set are not supported by the LI recognizer";
 640
 641   return(-1);
 642 }
 643
 644 static int
 645 li_recognizer_translate(recognizer r,
 646                          u_int ncs,
 647                          pen_stroke* tps,
 648                          bool correlate_p,
 649                          int* nret,
 650                          rec_alternative** ret)
 651 {
 652     char* clss = NULL;
 653     li_recognizer* rec;
 654     int conf;
 655     rClassifier* rc;
 656
 657     rec = (li_recognizer*)r->recognizer_specific;
 658
 659     *nret = 0;
 660     *ret = NULL;
 661
 662     /*Check for LI recognizer.*/
 663
 664     if( !CHECK_LI_MAGIC(rec) ) {
 665         li_err_msg = "Not a LI recognizer";
 666         return(-1);
 667     }
 668
 669    rc = &(rec->li_rc);
 670
 671     /*Check for valid parameters.*/
 672     if (ncs < 1) {
 673         li_err_msg = "Invalid parameters: ncs";
 674         return(-1);
 675     }
 676     if( tps == NULL) {
 677         li_err_msg = "Invalid parameters: tps";
 678         return(-1);
 679     }
 680     if( nret == NULL) {
 681         li_err_msg = "Invalid parameters: nret";
 682         return(-1);
 683     }
 684     if( ret == NULL) {
 685         li_err_msg = "Invalid parameters: ret";
 686         return(-1);
 687     }
 688
 689 /*    if( ncs < 1 || tps == NULL || nret == NULL || ret == NULL) {
 690         li_err_msg = "Invalid parameters";
 691         return(-1);
 692     }
 693 */
 694
 695     /*Check for null classifier. It must have at least one.*/
 696 /*
 697     if( rec->li_rc.sc == NULL ) {
 698         li_err_msg = "No classifier";
 699         return(-1);
 700     }
 701 */
 702
 703     /*
 704      * Go through the stroke array and recognize. Since this is a single
 705      *   stroke recognizer, each stroke is treated as a separate
 706      *   character or gesture. We allow only characters or gestures
 707      *   to be recognized at one time, since otherwise, handling
 708      *   the display of segmentation would be difficult.
 709     */
 710     clss = recognize_internal(rc,tps,&conf);
 711     if (clss == NULL) {
 712 /*
 713         li_err_msg = "unrecognized character";
 714         return(-1);
 715 */
 716         *nret = 1;
 717         return(0);
 718     }
 719
 720     /*Return values.*/
 721     *nret = 1;
 722     return(*clss);
 723 }
 724
 725
 726 static rec_fn*
 727 li_recognizer_get_extension_functions(recognizer rec)
 728 {
 729     rec_fn* ret;
 730
 731     /*Check for LI recognizer.*/
 732
 733     if( !CHECK_LI_MAGIC(rec->recognizer_specific) ) {
 734         li_err_msg = "Not a LI recognizer";
 735         return(NULL);
 736     }
 737
 738     ret = make_rec_fn_array(LI_NUM_EX_FNS);
 739
 740 /* ari -- clearState & getClasses are mine */
 741     ret[LI_GET_CLASSES] = (rec_fn)recognizer_getClasses;
 742     ret[LI_CLEAR] = (rec_fn)recognizer_clearState;
 743     ret[LI_ISA_LI] = (rec_fn)isa_li;
 744     ret[LI_TRAIN] = (rec_fn)recognizer_train;
 745
 746     return(ret);
 747 }
 748
 749 static char**
 750 li_recognizer_get_gesture_names(recognizer r)
 751 {
 752   /*This operation isn't supported by the LI recognizer.*/
 753
 754   li_err_msg = "Gestures are not supported by the LI recognizer";
 755
 756   return(NULL);
 757 }
 758
 759 static xgesture
 760 li_recognizer_set_gesture_action(recognizer r,
 761                                   char* name,
 762                                   xgesture fn,
 763                                   void* wsinfo)
 764 {
 765   /*This operation isn't supported by the LI recognizer.*/
 766
 767   li_err_msg = "Gestures are not supported by the LI recognizer";
 768
 769   return(NULL);
 770 }
 771
 772 /*
 773  * Exported Functions
 774 */
 775
 776 /*RECOGNIZER_INITIALIZE-Initialize the recognizer.*/
 777
 778 /* note from ari:  this expands via pre-processor to
 779  *
 780  * recognizer __recognizer_internal_initialize(rec_info* ri)
 781 */
 782
 783 RECOGNIZER_INITIALIZE(ri)
 784 {
 785     recognizer r;
 786     li_recognizer* rec;
 787     int i;
 788
 789     /*Check that locale matches.*/
 790
 791     if( strcmp(ri->ri_locale,LI_SUPPORTED_LOCALE) != 0 ) {
 792         li_err_msg = "Not a supported locale";
 793 fprintf(stderr, "Locale error.\n");
 794 #if 0
 795         return(NULL);
 796 #endif
 797     }
 798
 799     /*
 800      * Check that character sets match. Note that this is only approximate,
 801      * since the classifier file will have more information.
 802     */
 803
 804     if( ri->ri_subset != NULL ) {
 805       for(i = 0; ri->ri_subset[i] != NULL; i++ ) {
 806
 807         if( strcmp(ri->ri_subset[i],UPPERCASE) != 0 &&
 808             strcmp(ri->ri_subset[i],LOWERCASE) != 0  &&
 809             strcmp(ri->ri_subset[i],DIGITS) != 0 &&
 810             strcmp(ri->ri_subset[i],GESTURE) != 0 ) {
 811           li_err_msg = "Not a supported character set";
 812 fprintf(stderr, "charset error.\n");
 813
 814           return(NULL);
 815         }
 816       }
 817     }
 818
 819 /* ari */
 820     r = make_recognizer(ri);
 821     /*fprintf(stderr, "past make_recognizer.\n");*/
 822
 823     if( r == NULL ) {
 824         li_err_msg = "Can't allocate storage";
 825
 826         return(NULL);
 827     }
 828
 829     /*Make a LI recognizer structure.*/
 830
 831
 832     /* rec = (li_recognizer*)safe_malloc(sizeof(li_recognizer))) == NULL ); */
 833
 834     rec = allocate(1, li_recognizer);
 835
 836     r->recognizer_specific = rec;
 837
 838     rec->li_rc.file_name = NULL;
 839     rec->li_rc.sc = NULL;
 840
 841     /*Initialize the recognizer struct.*/
 842
 843     r->recognizer_load_state = li_recognizer_load;
 844     r->recognizer_save_state = li_recognizer_save;
 845     r->recognizer_load_dictionary = li_recognizer_load_dictionary;
 846     r->recognizer_save_dictionary = li_recognizer_save_dictionary;
 847     r->recognizer_free_dictionary = li_recognizer_free_dictionary;
 848     r->recognizer_add_to_dictionary = li_recognizer_add_to_dictionary;
 849     r->recognizer_delete_from_dictionary = li_recognizer_delete_from_dictionary;
 850     r->recognizer_error = li_recognizer_error;
 851     r->recognizer_translate = li_recognizer_translate;
 852     r->recognizer_get_context = li_recognizer_get_context;
 853     r->recognizer_set_context = li_recognizer_set_context;
 854     r->recognizer_get_buffer = li_recognizer_get_buffer;
 855     r->recognizer_set_buffer = li_recognizer_set_buffer;
 856     r->recognizer_clear = li_recognizer_clear;
 857     r->recognizer_get_extension_functions =
 858       li_recognizer_get_extension_functions;
 859     r->recognizer_get_gesture_names = li_recognizer_get_gesture_names;
 860     r->recognizer_set_gesture_action =
 861       li_recognizer_set_gesture_action;
 862
 863     /*Initialize LI Magic Number.*/
 864
 865     rec->li_magic = LI_MAGIC;
 866
 867     /*Initialize rClassifier.*/
 868
 869     rec->li_rc.file_name = NULL;
 870
 871     for( i = 0; i < MAXSCLASSES; i++ ) {
 872         rec->li_rc.ex[i] = NULL;
 873         rec->li_rc.cnames[i] = NULL;
 874     }
 875
 876     lialg_initialize(&rec->li_rc);
 877
 878     /*Get rid of error message. Not needed here.*/
 879     li_err_msg = NULL;
 880
 881     return(r);
 882 }
 883
 884 /*free_rClassifier-Free the rClassifier.*/
 885
 886 static void
 887 free_rClassifier(rClassifier* rc)
 888 {
 889     int i;
 890
 891     if( rc->file_name != NULL) {
 892         free(rc->file_name);
 893     }
 894
 895     for( i = 0; rc->ex[i] != NULL; i++) {
 896         delete_examples(rc->ex[i]);
 897         free(rc->cnames[i]);
 898     }
 899
 900     if(rc->sc != NULL ) {
 901         sFreeClassifier(rc->sc);
 902     }
 903 }
 904
 905 /*RECOGNIZER_FINALIZE-Deallocate the recognizer, finalize.*/
 906
 907 RECOGNIZER_FINALIZE(r)
 908 {
 909     li_recognizer* rec = (li_recognizer*)r->recognizer_specific;
 910
 911     /*Make sure this is a li_recognizer first*/
 912
 913     if( !CHECK_LI_MAGIC(rec) ) {
 914         li_err_msg = "Not a LI recognizer";
 915         return(-1);
 916     }
 917
 918     /*Deallocate rClassifier resources.*/
 919
 920     free_rClassifier(&(rec->li_rc));
 921
 922     /*Deallocate the li_recognizer struct.*/
 923
 924     free(rec);
 925
 926     /*Deallocate the recognizer*/
 927
 928     delete_recognizer(r);
 929
 930     return(0);
 931 }
 932
 933
 934 /* **************************************************
 935
 936   Implementation of the Li/Yeung recognition algorithm
 937
 938 ************************************************** */
 939
 940 /*#include <assert.h>*/
 941 #ifdef __ECOS
 942 #define MAXINT 0x7FFFFFFF
 943 #else
 944 #include <values.h>
 945 #endif
 946 #include <sys/time.h>
 947
 948 #ifdef  __ultrix
 949 /* Ultrix doesn't have these declarations in math.h! */
 950 extern double rint(double);
 951 extern float expf(float);
 952 #endif
 953
 954 #ifdef  ELX
 955 extern double rint (double);
 956 extern float expf (float);      /* N.B.  exp() appears to be broken on ELX! */
 957 #endif
 958
 959 #define WORST_SCORE     MAXINT
 960
 961 /* Dynamic programming parameters */
 962 #define DP_BAND         3
 963 #define MIN_SIM         0
 964 #define MAX_DIST        MAXINT
 965 #define SIM_THLD        60      /* x 100 */
 966 #define DIST_THLD       3200    /* x 100 */
 967
 968 /* Low-pass filter parameters -- empirically derived */
 969 #define LP_FILTER_WIDTH 6
 970 #define LP_FILTER_ITERS 8
 971 #define LP_FILTER_THLD  250     /* x 100 */
 972 #define LP_FILTER_MIN   5
 973
 974 /* Pseudo-extrema parameters -- empirically derived */
 975 #define PE_AL_THLD      1500    /* x 100 */
 976 #define PE_ATCR_THLD    135     /* x 100 */
 977
 978 /* Contour-angle derivation parameters */
 979 #define T_ONE           1
 980 #define T_TWO           20
 981
 982 /* Pre-processing and canonicalization parameters */
 983 #define CANONICAL_X     108
 984 #define CANONICAL_Y     128
 985 #define DIST_SQ_THRESHOLD   (3*3)       /* copied from fv.h */
 986 #define NCANONICAL      50
 987
 988 /* Tap-handling parameters */
 989 #define TAP_CHAR        "."
 990 #define TAP_TIME_THLD   150         /* msec */
 991 #define TAP_DIST_THLD   75          /* dx * dx + dy * dy */
 992 #define TAP_PATHLEN     1000        /* x 100 */
 993
 994
 995 /* Overload the time field of the pen_point struct with the chain-code. */
 996 #define chaincode   time
 997
 998 /* region types */
 999 #define RGN_CONVEX  0
1000 #define RGN_CONCAVE 1
1001 #define RGN_PLAIN   2
1002 #define RGN_PSEUDO  3
1003
1004
1005 typedef struct RegionList {
1006     int start;
1007     int end;
1008     int type;
1009     struct RegionList *next;
1010 } region_list;
1011
1012
1013 /* direction-code table; indexed by dx, dy */
1014 static int lialg_dctbl[3][3] = {{1, 0, 7}, {2, 0x7FFFFFFF, 6}, {3, 4, 5}};
1015
1016 /* low-pass filter weights */
1017 static int lialg_lpfwts[2 * LP_FILTER_WIDTH + 1];
1018 static int lialg_lpfconst = -1;
1019
1020
1021 static int lialg_preprocess_stroke(point_list *);
1022 static point_list *lialg_compute_dominant_points(point_list *);
1023 static point_list *lialg_interpolate_points(point_list *);
1024 static void lialg_bresline(pen_point *, pen_point *, point_list *, int *);
1025 static void lialg_compute_chain_code(point_list *);
1026 static void lialg_compute_unit_chain_code(point_list *);
1027 static region_list *lialg_compute_regions(point_list *);
1028 static point_list *lialg_compute_dompts(point_list *, region_list *);
1029 static int *lialg_compute_contour_angle_set(point_list *, region_list *);
1030 static void lialg_score_stroke(point_list *, point_list *, int *, int *);
1031 static int lialg_compute_similarity(point_list *, point_list *);
1032 static int lialg_compute_distance(point_list *, point_list *);
1033
1034 static int lialg_read_classifier_digest(rClassifier *);
1035
1036 static int lialg_canonicalize_examples(rClassifier *);
1037 static int lialg_canonicalize_example_stroke(point_list *);
1038 static int lialg_compute_equipoints(point_list *);
1039
1040 static int lialg_compute_pathlen(point_list *);
1041 static int lialg_compute_pathlen_subset(point_list *, int, int);
1042 static int lialg_filter_points(point_list *);
1043 static int lialg_translate_points(point_list *, int, int, int, int);
1044 static void lialg_get_bounding_box(point_list *, int *, int *, int *, int *);
1045 static void lialg_compute_lpf_parameters();
1046 static int isqrt(int);
1047 static int likeatan(int, int);
1048 static int quadr(int);
1049
1050
1051 /*************************************************************
1052
1053   Core routines for the Li/Yeung recognition algorithm
1054
1055  *************************************************************/
1056
1057 static void lialg_initialize(rClassifier *rec) {
1058     int i;
1059
1060     /* Initialize the dompts arrays. */
1061     for (i = 0; i < MAXSCLASSES; i++) {
1062         rec->dompts[i] = NULL;
1063     }
1064 }
1065
1066
1067 /*
1068  *  Main recognition routine -- called by HRE API.
1069  */
1070 static char *lialg_recognize_stroke(rClassifier *rec, point_list *stroke) {
1071     int i;
1072     char *name = NULL;
1073     point_list *input_dompts = NULL;
1074     char *best_name = NULL;
1075     int best_score = WORST_SCORE;
1076     char *curr_name;
1077     point_list *curr_dompts = NULL;
1078     /*struct timeval stv, etv;
1079     int elapsed;*/
1080
1081     /*    (void)gettimeofday(&stv, NULL);*/
1082
1083     if (stroke->npts < 1) goto done;
1084
1085     /* Check for tap. */
1086     {
1087 /*
1088         pen_point *startpt = &stroke->pts[0];
1089         pen_point *endpt = &stroke->pts[stroke->npts - 1];
1090         int dt = endpt->time - startpt->time;
1091         int dx = endpt->x - startpt->x;
1092         int dy = endpt->y - startpt->y;
1093         int magsq = dx * dx + dy * dy;
1094 */
1095
1096         /* First thing is to filter out ``close points.'' */
1097         if (lialg_filter_points(stroke) != 0) return(NULL);
1098
1099         /* Unfortunately, we don't have the actual time that each point */
1100         /* was recorded (i.e., dt is invalid).  Hence, we have to use a */
1101         /* heuristic based on total distance and the number of points. */
1102         if (stroke->npts == 1 || lialg_compute_pathlen(stroke) < TAP_PATHLEN)
1103             return(TAP_CHAR);
1104     }
1105
1106     /* Pre-process input stroke. */
1107     if (lialg_preprocess_stroke(stroke) != 0) goto done;
1108
1109     /* Compute its dominant points. */
1110     input_dompts = lialg_compute_dominant_points(stroke);
1111     if (input_dompts == NULL) goto done;
1112
1113     /* Score input stroke against every class in classifier. */
1114     for (i = 0, curr_name = rec->cnames[i], curr_dompts = rec->dompts[i];
1115           i < MAXSCLASSES && curr_name != NULL && curr_dompts != NULL;
1116           i++, curr_name = rec->cnames[i], curr_dompts = rec->dompts[i]) {
1117         int sim;
1118         int dist;
1119         int curr_score;
1120
1121         lialg_score_stroke(input_dompts, curr_dompts, &sim, &dist);
1122         curr_score = dist;
1123
1124         if (lidebug && curr_score < DIST_THLD)
1125             fprintf(stderr, "(%s, %d, %d) ", curr_name, sim, dist);
1126
1127         /* Is it the best so far? */
1128         if (curr_score < best_score && curr_score <= DIST_THLD) {
1129             best_score = curr_score;
1130             best_name = curr_name;
1131         }
1132     }
1133
1134     if (lidebug)
1135         fprintf(stderr, "\n");
1136
1137     /* No errors. */
1138     name = best_name;
1139
1140 done:
1141     delete_examples(input_dompts);
1142     /*    (void)gettimeofday(&etv, NULL);
1143           elapsed = (1000 * (etv.tv_sec - stv.tv_sec)) + ((etv.tv_usec - stv.tv_usec + 500) / 1000);
1144           fprintf(stderr, "elapsed = %d\n", elapsed);
1145      */
1146     return(name);
1147 }
1148
1149
1150 static int lialg_preprocess_stroke(point_list *points) {
1151     int minx, miny, maxx, maxy, xrange, yrange, scale, xoff, yoff;
1152
1153     /* Filter out points that are too close. */
1154     /* We did this earlier, when we checked for a tap. */
1155 /*
1156     if (lialg_filter_points(points) != 0) return(-1);
1157 */
1158
1159 /*    assert(points->npts > 0);*/
1160
1161     /* Scale up to avoid conversion errors. */
1162     lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
1163     xrange = maxx - minx;
1164     yrange = maxy - miny;
1165     scale = ( ((100 * xrange + CANONICAL_X / 2) / CANONICAL_X) >
1166               ((100 * yrange + CANONICAL_Y / 2) / CANONICAL_Y))
1167       ? (100 * CANONICAL_X + xrange / 2) / xrange
1168       : (100 * CANONICAL_Y + yrange / 2) / yrange;
1169     if (lialg_translate_points(points, minx, miny, scale, scale) != 0)
1170         return(-1);
1171
1172     /* Center the stroke. */
1173     lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
1174     xrange = maxx - minx;
1175     yrange = maxy - miny;
1176     xoff = -((CANONICAL_X - xrange + 1) / 2);
1177     yoff = -((CANONICAL_Y - yrange + 1) / 2);
1178     if (lialg_translate_points(points, xoff, yoff, 100, 100) != 0) return(-1);
1179
1180     /* Store the x and y ranges in the point list. */
1181     xrange = maxx - minx;
1182     yrange = maxy - miny;
1183     points->xrange = xrange;
1184     points->yrange = yrange;
1185
1186     if (lidebug) {
1187         int i;
1188         fprintf(stderr, "After pre-processing:   %d %d %d %d\n",
1189                 minx, miny, maxx, maxy);
1190         for (i = 0; i < points->npts; i++)
1191             fprintf(stderr, "      (%d %d)\n",
1192                     points->pts[i].x, points->pts[i].y);
1193         fflush(stderr);
1194     }
1195
1196     return(0);
1197 }
1198
1199
1200 static point_list *lialg_compute_dominant_points(point_list *points) {
1201     point_list *ipts = NULL;
1202     region_list *regions = NULL;
1203     point_list *dpts = NULL;
1204
1205     /* Interpolate points. */
1206     ipts = lialg_interpolate_points(points);
1207     if (ipts == NULL) return(NULL);
1208     if (lidebug) {
1209         int j;
1210         fprintf(stderr, "After interpolation:  %d ipts\n", ipts->npts);
1211         for (j = 0; j < ipts->npts; j++) {
1212             fprintf(stderr, "  (%d, %d), %ld\n",
1213                     ipts->pts[j].x, ipts->pts[j].y, ipts->pts[j].chaincode);
1214         }
1215         fflush(stderr);
1216     }
1217
1218     /* Compute regions. */
1219     regions = lialg_compute_regions(ipts);
1220 /*    assert(regions != NULL);*/
1221
1222     /* Compute dominant points. */
1223     dpts = lialg_compute_dompts(ipts, regions);
1224     if (lidebug) {
1225         int j;
1226         fprintf(stderr, "Dominant points:  ");
1227         for (j = 0; j < dpts->npts; j++) {
1228             fprintf(stderr, "%d %d (%ld)  ",
1229                     dpts->pts[j].x, dpts->pts[j].y, dpts->pts[j].chaincode);
1230         }
1231         fprintf(stderr, "\n");
1232         fflush(stderr);
1233     }
1234
1235     /* Delete region data structure. */
1236     {
1237         region_list *curr, *next;
1238         for (curr = regions; curr != NULL; ) {
1239             next = curr->next;
1240             free(curr);
1241             curr = next;
1242         }
1243     }
1244     delete_examples(ipts);
1245     return(dpts);
1246 }
1247
1248
1249 /* Input points are assumed to be integer-valued! */
1250 static point_list *lialg_interpolate_points(point_list *points) {
1251     int i, j;
1252     int maxpts;
1253     point_list *newpts;
1254
1255     /* Compute an upper-bound on the number of interpolated points. */
1256     maxpts = 0;
1257     for (i = 0; i < (points->npts - 1); i++) {
1258         pen_point *pta = &(points->pts[i]);
1259         pen_point *ptb = &(points->pts[i+1]);
1260         maxpts += abs(pta->x - ptb->x) + abs(pta->y - ptb->y);
1261     }
1262
1263     /* Allocate an array of the requisite size. */
1264     maxpts += points->npts;
1265     /* newpts = (point_list *)safe_malloc(sizeof(point_list)); */
1266     newpts = allocate(1, point_list);
1267     newpts->pts = make_pen_point_array(maxpts);
1268     if (newpts->pts == NULL) {
1269         free(newpts);
1270         return(NULL);
1271     }
1272     newpts->npts = maxpts;
1273     newpts->next = NULL;
1274
1275     /* Interpolate each of the segments. */
1276     j = 0;
1277     for (i = 0; i < (points->npts - 1); i++) {
1278         pen_point *startpt = &(points->pts[i]);
1279         pen_point *endpt = &(points->pts[i+1]);
1280
1281         lialg_bresline(startpt, endpt, newpts, &j);
1282
1283         j--;    /* end point gets recorded as start point of next segment! */
1284     }
1285
1286     /* Add-in last point. */
1287     newpts->pts[j++] = points->pts[points->npts - 1];
1288     newpts->npts = j;
1289
1290     /* Compute the chain code for P (the list of points). */
1291     lialg_compute_unit_chain_code(newpts);
1292
1293     return(newpts);
1294 }
1295
1296
1297 /* This implementation is due to Kenny Hoff. */
1298 static void lialg_bresline(pen_point *startpt, pen_point *endpt,
1299                             point_list *newpts, int *j) {
1300     int Ax, Ay, Bx, By, dX, dY, Xincr, Yincr;
1301
1302     Ax = startpt->x;
1303     Ay = startpt->y;
1304     Bx = endpt->x;
1305     By = endpt->y;
1306
1307     /* INITIALIZE THE COMPONENTS OF THE ALGORITHM THAT ARE NOT AFFECTED */
1308     /* BY THE SLOPE OR DIRECTION OF THE LINE */
1309     dX = abs(Bx-Ax);    /* store the change in X and Y of the line endpoints */
1310     dY = abs(By-Ay);
1311
1312     /* DETERMINE "DIRECTIONS" TO INCREMENT X AND Y (REGARDLESS OF DECISION) */
1313     if (Ax > Bx) { Xincr=-1; } else { Xincr=1; }    /* which direction in X? */
1314     if (Ay > By) { Yincr=-1; } else { Yincr=1; }    /* which direction in Y? */
1315
1316     /* DETERMINE INDEPENDENT VARIABLE (ONE THAT ALWAYS INCREMENTS BY 1 (OR -1) ) */
1317     /* AND INITIATE APPROPRIATE LINE DRAWING ROUTINE (BASED ON FIRST OCTANT */
1318     /* ALWAYS). THE X AND Y'S MAY BE FLIPPED IF Y IS THE INDEPENDENT VARIABLE. */
1319     if (dX >= dY) {         /* if X is the independent variable */
1320         int dPr = dY<<1;            /* amount to increment decision if right is chosen (always) */
1321         int dPru = dPr - (dX<<1);   /* amount to increment decision if up is chosen */
1322         int P = dPr - dX;           /* decision variable start value */
1323
1324         /* process each point in the line one at a time (just use dX) */
1325         for (; dX>=0; dX--) {
1326             newpts->pts[*j].x = Ax;
1327             newpts->pts[*j].y = Ay;
1328             (*j)++;
1329
1330             if (P > 0) {        /* is the pixel going right AND up? */
1331                 Ax+=Xincr;      /* increment independent variable */
1332                 Ay+=Yincr;      /* increment dependent variable */
1333                 P+=dPru;        /* increment decision (for up) */
1334             }
1335             else {              /* is the pixel just going right? */
1336                 Ax+=Xincr;      /* increment independent variable */
1337                 P+=dPr;         /* increment decision (for right) */
1338             }
1339         }
1340     }
1341     else {                  /* if Y is the independent variable */
1342         int dPr = dX<<1;            /* amount to increment decision if right is chosen (always) */
1343         int dPru = dPr - (dY<<1);   /* amount to increment decision if up is chosen */
1344         int P  = dPr - dY;          /* decision variable start value */
1345
1346         /* process each point in the line one at a time (just use dY) */
1347         for (; dY>=0; dY--) {
1348             newpts->pts[*j].x = Ax;
1349             newpts->pts[*j].y = Ay;
1350             (*j)++;
1351
1352             if (P > 0) {        /* is the pixel going up AND right? */
1353                 Ax+=Xincr;      /* increment dependent variable */
1354                 Ay+=Yincr;      /* increment independent variable */
1355                 P+=dPru;        /* increment decision (for up) */
1356             }
1357             else {              /* is the pixel just going up? */
1358                 Ay+=Yincr;      /* increment independent variable */
1359                 P+=dPr;         /* increment decision (for right) */
1360             }
1361         }
1362     }
1363 }
1364
1365
1366 static void lialg_compute_chain_code(point_list *pts) {
1367     int i;
1368
1369     for (i = 0; i < (pts->npts - 1); i++) {
1370         pen_point *startpt = &(pts->pts[i]);
1371         pen_point *endpt = &(pts->pts[i+1]);
1372         int dx = endpt->x - startpt->x;
1373         int dy = endpt->y - startpt->y;
1374 /*
1375         int tmp = rint(4.0 * atan2((double)dx, (double)dy) / M_PI);
1376         int dircode = (10 + tmp) % 8;
1377 */
1378         int tmp = quadr(likeatan(dy, dx));
1379         int dircode = (12 - tmp) % 8;
1380
1381         startpt->chaincode = dircode;
1382     }
1383 }
1384
1385
1386 static void lialg_compute_unit_chain_code(point_list *pts) {
1387     int i;
1388
1389     for (i = 0; i < (pts->npts - 1); i++) {
1390         pen_point *startpt = &(pts->pts[i]);
1391         pen_point *endpt = &(pts->pts[i+1]);
1392         int dx = endpt->x - startpt->x;
1393         int dy = endpt->y - startpt->y;
1394         int dircode = lialg_dctbl[dx+1][dy+1];
1395
1396 /*      assert(dircode < 8);*/
1397         startpt->chaincode = dircode;
1398     }
1399 }
1400
1401
1402 static region_list *lialg_compute_regions(point_list *pts) {
1403     region_list *regions = NULL;
1404     region_list *curr_reg = NULL;
1405     int *R[2 + LP_FILTER_ITERS];
1406     int *junk;
1407     int *curr, *next;
1408     int i, j;
1409
1410     /* Initialize low-pass filter parameters if necessary. */
1411     if (lialg_lpfconst == -1)
1412         lialg_compute_lpf_parameters();
1413
1414     /* Allocate a 2 x pts->npts array for use in computing the (filtered) Angle set, A_n. */
1415     /*    junk = (int *)safe_malloc((2 + LP_FILTER_ITERS) * pts->npts * sizeof(int)); */
1416     junk = allocate((2 + LP_FILTER_ITERS) * pts->npts, int);
1417     for (i = 0; i < (2 + LP_FILTER_ITERS); i++)
1418         R[i] = junk + (i * pts->npts);
1419     curr = R[0];
1420
1421     /* Compute the Angle set, A, in the first element of array R. */
1422     /* Values in R are in degrees, x 100. */
1423     curr[0] = 18000;                            /* a_0 */
1424     for (i = 1; i < (pts->npts - 1); i++) {
1425         int d_i = pts->pts[i].chaincode;
1426         int d_iminusone = pts->pts[i-1].chaincode;
1427         int a_i;
1428
1429         if (d_iminusone < d_i)
1430             d_iminusone += 8;
1431
1432         a_i = (d_iminusone - d_i) % 8;
1433
1434         /* convert to degrees, x 100 */
1435         curr[i] = ((12 - a_i) % 8) * 45 * 100;
1436     }
1437     curr[pts->npts - 1] = 18000;                /* a_L-1 */
1438
1439     /* Perform a number of filtering iterations. */
1440     next = R[1];
1441     for (j = 0; j < LP_FILTER_ITERS; j++, curr = R[j], next = R[j+1]) {
1442         for (i = 0; i < pts->npts; i++) {
1443             int k;
1444
1445             next[i] = 0;
1446
1447             for (k = i - LP_FILTER_WIDTH; k <= i + LP_FILTER_WIDTH; k++) {
1448                 int oldval = (k < 0 || k >= pts->npts) ? 18000 : curr[k];
1449                 next[i] += oldval * lialg_lpfwts[k - (i - LP_FILTER_WIDTH)];    /* overflow? */
1450             }
1451
1452             next[i] /= lialg_lpfconst;
1453         }
1454     }
1455
1456     /* Do final thresholding around PI. */
1457     /* curr and next are set-up correctly at end of previous loop! */
1458     for (i = 0; i < pts->npts; i++) {
1459         next[i] = (abs(curr[i] - 18000) < LP_FILTER_THLD)
1460           ? 18000
1461           : curr[i];
1462     }
1463     curr = next;
1464
1465     /* Debugging. */
1466     if (lidebug > 1) {
1467         for (i = 0; i < pts->npts; i++) {
1468             fprintf(stderr, "%3d:  (%3d, %3d)  %ld  ",
1469                     i, pts->pts[i].x, pts->pts[i].y, pts->pts[i].chaincode);
1470             for (j = 0; j < 2 + LP_FILTER_ITERS; j++)
1471                 fprintf(stderr, "%d  ", R[j][i]);
1472             fprintf(stderr, "\n");
1473         }
1474     }
1475
1476     /* Do the region segmentation. */
1477     {
1478         int start, end;
1479         int currtype;
1480
1481 #define RGN_TYPE(val)\
1482     (((val) == 18000)\
1483         ? RGN_PLAIN\
1484         : ((val) < 18000 ? RGN_CONCAVE : RGN_CONVEX))
1485
1486         start = 0;
1487         currtype = RGN_TYPE(curr[0]);
1488         /*      regions = (region_list *)safe_malloc(sizeof(region_list));*/
1489         regions = allocate(1, region_list);
1490         curr_reg = regions;
1491         curr_reg->start = start;
1492         curr_reg->end = 0;
1493         curr_reg->type = currtype;
1494         curr_reg->next = NULL;
1495         for (i = 1; i < pts->npts; i++) {
1496             int nexttype = RGN_TYPE(curr[i]);
1497
1498             if (nexttype != currtype) {
1499                 region_list *next_reg = NULL;
1500
1501                 end = i - 1;
1502                 curr_reg->end = end;
1503                 if (lidebug > 1)
1504                     fprintf(stderr, "  (%d, %d), %d\n", start, end, currtype);
1505
1506                 start = i;
1507                 currtype = nexttype;
1508                 /* next_reg = (region_list *)safe_malloc(sizeof(region_list));*/
1509                 next_reg = allocate(1, region_list);
1510                 next_reg->start = start;
1511                 next_reg->end = 0;
1512                 next_reg->type = nexttype;
1513                 next_reg->next = NULL;
1514
1515                 curr_reg->next = next_reg;
1516                 curr_reg = next_reg;
1517             }
1518         }
1519         end = i - 1;
1520         curr_reg->end = end;
1521         if (lidebug > 1)
1522             fprintf(stderr, "  (%d, %d), %d\n", start, end, currtype);
1523
1524         /* Filter out convex/concave regions that are too short. */
1525         for (curr_reg = regions; curr_reg; curr_reg = curr_reg->next)
1526             if (curr_reg->type == RGN_PLAIN) {
1527                 region_list *next_reg;
1528
1529                 for (next_reg = curr_reg->next;
1530                      next_reg != NULL &&
1531                        (next_reg->end - next_reg->start) < LP_FILTER_MIN;
1532                      next_reg = curr_reg->next) {
1533                     /* next_reg must not be plain, and it must be followed by a plain */
1534                     /* assert(next_reg->type != RGN_PLAIN); */
1535                     /* assert(next_reg->next != NULL && (next_reg->next)->type == RGN_PLAIN); */
1536
1537                     curr_reg->next = (next_reg->next)->next;
1538                     curr_reg->end = (next_reg->next)->end;
1539
1540                     free(next_reg->next);
1541                     free(next_reg);
1542                 }
1543             }
1544
1545         /* Add-in pseudo-extremes. */
1546         {
1547             region_list *tmp, *prev_reg;
1548
1549             tmp = regions;
1550             regions = NULL;
1551             prev_reg = NULL;
1552             for (curr_reg = tmp; curr_reg; curr_reg = curr_reg->next) {
1553                 if (curr_reg->type == RGN_PLAIN) {
1554                     int arclen = lialg_compute_pathlen_subset(pts,
1555                                                                 curr_reg->start,
1556                                                                 curr_reg->end);
1557                     int dx = pts->pts[curr_reg->end].x -
1558                       pts->pts[curr_reg->start].x;
1559                     int dy = pts->pts[curr_reg->end].y -
1560                       pts->pts[curr_reg->start].y;
1561                     int chordlen = isqrt(10000 * (dx * dx + dy * dy));
1562                     int atcr = (chordlen == 0) ? 0 : (100 * arclen + chordlen / 2) / chordlen;
1563
1564                     if (lidebug)
1565                         fprintf(stderr, "%d, %d, %d\n", arclen, chordlen, atcr);
1566
1567                     /* Split region if necessary. */
1568                     if (arclen >= PE_AL_THLD && atcr >= PE_ATCR_THLD) {
1569                         int mid = curr_reg->start + (curr_reg->end - curr_reg->start) / 2;
1570                         int end = curr_reg->end;
1571                         region_list *saved_next = curr_reg->next;
1572
1573                         curr_reg->end = mid - 1;
1574                         if (prev_reg == NULL)
1575                             regions = curr_reg;
1576                         else
1577                             prev_reg->next = curr_reg;
1578                         prev_reg = curr_reg;
1579
1580                         /* curr_reg = (region_list *)safe_malloc(sizeof(region_list));*/
1581                         curr_reg = allocate(1, region_list);
1582                         curr_reg->start = mid;
1583                         curr_reg->end = mid;
1584                         curr_reg->type = RGN_PSEUDO;
1585                         curr_reg->next = NULL;
1586                         prev_reg->next = curr_reg;
1587                         prev_reg = curr_reg;
1588
1589                         /* curr_reg = (region_list *)safe_malloc(sizeof(region_list)); */
1590                         curr_reg = allocate(1, region_list);
1591                         curr_reg->start = mid + 1;
1592                         curr_reg->end = end;
1593                         curr_reg->type = RGN_PLAIN;
1594                         curr_reg->next = NULL;
1595                         prev_reg->next = curr_reg;
1596                         prev_reg = curr_reg;
1597
1598                         curr_reg->next = saved_next;
1599                         continue;
1600                     }
1601                 }
1602
1603                 if (prev_reg == NULL)
1604                     regions = curr_reg;
1605                 else
1606                     prev_reg->next = curr_reg;
1607                 prev_reg = curr_reg;
1608             }
1609         }
1610     }
1611
1612     free(junk);
1613     return(regions);
1614 }
1615
1616
1617 static point_list *lialg_compute_dompts(point_list *pts, region_list *regions) {
1618     point_list *dpts = NULL;
1619     int ndpts;
1620     int *cas = NULL;
1621     int nonplain;
1622     region_list *r;
1623
1624     /* Compute contour angle set. */
1625     cas = lialg_compute_contour_angle_set(pts, regions);
1626 /*    assert(cas != NULL);*/
1627
1628     /* Dominant points include:  start_pt, end_pt, extrema_of_non_plain_regions, midpts of the preceding. */
1629     nonplain = 0;
1630     for (r = regions; r != NULL; r = r->next)
1631         if (r->type != RGN_PLAIN) nonplain++;
1632     ndpts = 2 * (2 + nonplain) - 1;
1633     /* dpts = (point_list *)safe_malloc(sizeof(point_list)); */
1634     dpts = allocate(1, point_list);
1635     dpts->pts = make_pen_point_array(ndpts);
1636     if (dpts->pts == NULL) {
1637         free(dpts);
1638         return(NULL);
1639     }
1640     dpts->npts = ndpts;
1641     dpts->next = NULL;
1642
1643     /* Pick out dominant points. */
1644     {
1645         region_list *curr;
1646         int dp;
1647         int previx;
1648         int currix;
1649
1650         /* Record start point. */
1651         dp = 0;
1652         previx = 0;
1653         dpts->pts[dp++] = pts->pts[previx];
1654
1655         for (curr = regions; curr != NULL; curr = curr->next)
1656             if (curr->type != RGN_PLAIN) {
1657                 int max_v = 0;
1658                 int min_v = MAXINT;
1659                 int max_ix = -1;
1660                 int min_ix = -1;
1661                 int i;
1662
1663                 for (i = curr->start; i <= curr->end; i++) {
1664                     int v = cas[i];
1665                     if (v > max_v) { max_v = v; max_ix = i; }
1666                     if (v < min_v) { min_v = v; min_ix = i; }
1667                     if (lidebug > 1)
1668                         fprintf(stderr, "  %d\n", v);
1669                 }
1670
1671                 currix = (curr->type == RGN_CONVEX ? max_ix : min_ix);
1672
1673                 /* Record midpoint. */
1674                 dpts->pts[dp++] = pts->pts[previx + (currix - previx) / 2];
1675
1676                 /* Record extreme point. */
1677                 dpts->pts[dp++] = pts->pts[currix];
1678
1679                 previx = currix;
1680             }
1681
1682         /* Record last mid-point and end point. */
1683         currix = pts->npts - 1;
1684         dpts->pts[dp++] = pts->pts[previx + (currix - previx) / 2];
1685         dpts->pts[dp++] = pts->pts[currix];
1686     }
1687
1688     /* Compute chain-code. */
1689     lialg_compute_chain_code(dpts);
1690
1691     free(cas);
1692     return(dpts);
1693 }
1694
1695
1696 static int *lialg_compute_contour_angle_set(point_list *pts,
1697                                                region_list *regions) {
1698     int *V = NULL;
1699     region_list *curr_reg, *prev_reg;
1700     int i;
1701
1702     /*    V = (int *)safe_malloc(pts->npts * sizeof(int));*/
1703     V = allocate(pts->npts, int);
1704
1705     V[0] = 18000;
1706     for (curr_reg = regions; curr_reg != NULL;
1707             prev_reg = curr_reg, curr_reg = curr_reg->next) {
1708         for (i = curr_reg->start; i <= curr_reg->end; i++) {
1709             if (curr_reg->type == RGN_PLAIN) {
1710                 V[i] = 18000;
1711             }
1712             else {
1713 #ifdef  notdef
1714                 /* XXX - eliminate floating point */
1715                 region_list *next_reg = curr_reg->next;
1716                 int b = curr_reg->start;
1717                 int h = prev_reg->start;
1718                 int t = next_reg->end;
1719                 int pts_before = i - h;
1720                 int pts_after = t - i;
1721                 int min_pts = (pts_before < pts_after)
1722                   ? pts_before
1723                   : pts_after;
1724                 int k = (min_pts < T_ONE)
1725                   ? T_ONE
1726                   : (min_pts > T_TWO)
1727                   ? T_TWO
1728                   : min_pts;
1729                 float sum = 0.0;
1730
1731                 for (j = 1; j <= k; j++) {
1732                     int ptA = i - j;
1733                     int ptB = i + j - 1;
1734                     int d_A = pts->pts[ptA].chaincode;
1735                     int d_B = pts->pts[ptB].chaincode;
1736                     int a_i;
1737
1738                     if (d_A < d_B)
1739                         d_A += 8;
1740
1741                     a_i = (d_A - d_B) % 8;
1742
1743                     /* convert to radians */
1744                     if (a_i == 4 && curr_reg->type == RGN_CONVEX)
1745                         sum += M_2_PI;
1746                     else
1747                         sum += (float)((12 - a_i) % 8) / 4.0 * M_PI;
1748                 }
1749                 V[i] = sum / (float)k;
1750 #else
1751                 /* For now, simply choose the mid-point. */
1752                 int isMidPt = (i == (curr_reg->start +
1753                                      (curr_reg->end - curr_reg->start) / 2));
1754                 V[i] = (curr_reg->type == RGN_CONVEX)
1755                   ? (isMidPt ? 18000 : 0)
1756                   : (isMidPt ? 0 : 18000);
1757 #endif
1758             }
1759         }
1760     }
1761     V[pts->npts - 1] = 18000;
1762
1763     return(V);
1764 }
1765
1766
1767 /*
1768  *  First compute the similarity between the two strings.
1769  *  If it's above a threshold, compute the distance between
1770  *  the two and return it as the ``score.''
1771  *  Otherwise, return the constant WORST_SCORE.
1772  *
1773  */
1774 static void lialg_score_stroke(point_list *input_dompts, point_list *curr_dompts, int *sim, int *dist) {
1775     *sim = MIN_SIM;
1776     *dist = MAX_DIST;
1777
1778     *sim = lialg_compute_similarity(input_dompts, curr_dompts);
1779     if (*sim < SIM_THLD) goto done;
1780
1781     *dist = lialg_compute_distance(input_dompts, curr_dompts);
1782
1783 done:
1784     if (lidebug)
1785         fprintf(stderr, "%d, %d\n", *sim, *dist);
1786 }
1787
1788
1789 static int lialg_compute_similarity(point_list *input_dompts,
1790                                      point_list *curr_dompts) {
1791     int sim = 0;
1792     point_list *A, *B;
1793     int N, M;
1794     int **G = NULL;
1795     int *junk = NULL;
1796     int i, j;
1797
1798     /* A is the longer sequence, length N. */
1799     /* B is the shorter sequence, length M. */
1800     if (input_dompts->npts >= curr_dompts->npts) {
1801         A = input_dompts;
1802         N = input_dompts->npts;
1803         B = curr_dompts;
1804         M = curr_dompts->npts;
1805     }
1806     else {
1807         A = curr_dompts;
1808         N = curr_dompts->npts;
1809         B = input_dompts;
1810         M = input_dompts->npts;
1811     }
1812
1813     /* Allocate and initialize the Gain matrix, G. */
1814     /* The size of G is M x (N + 1). */
1815     /* Note that row 0 is unused. */
1816     /* Similarities are x 10. */
1817     {
1818         /*      G = (int **)safe_malloc(M * sizeof(int *));*/
1819         G = allocate(M, int *);
1820         /*      junk = (int *)safe_malloc(M * (N + 1) * sizeof(int)); */
1821         junk = allocate(M * (N + 1), int);
1822         for (i = 0; i < M; i++)
1823             G[i] = junk + (i * (N + 1));
1824
1825         for (i = 1; i < M; i++) {
1826             int bval = B->pts[i-1].chaincode;
1827
1828             /* Source column. */
1829             G[i][0] = 0;
1830
1831             for (j = 1; j < N; j++) {
1832                 int aval = A->pts[j-1].chaincode;
1833                 int diff = abs(bval - aval);
1834                 if (diff > 4) diff = 8 - diff;
1835
1836                 G[i][j] = (diff == 0)
1837                   ? 10
1838                   : (diff == 1)
1839                   ? 6
1840                   : 0;
1841             }
1842
1843             /* Sink column. */
1844             G[i][N] = 0;
1845         }
1846     }
1847
1848     /* Do the DP algorithm. */
1849     /* Proceed in column order, from highest column to the lowest. */
1850     /* Within each column, proceed from the highest row to the lowest. */
1851     /* Skip the highest column. */
1852     {
1853         for (j = N - 1; j >= 0; j--)
1854             for (i = M - 1; i > 0; i--) {
1855                 int max = G[i][j + 1];
1856
1857                 if (i < (M - 1)) {
1858                     int tmp = G[i + 1][j + 1];
1859                     if (tmp > max) max = tmp;
1860                 }
1861
1862                 G[i][j] += max;
1863             }
1864
1865         sim = (10 * G[1][0] + (N - 1) / 2) / (N - 1);
1866     }
1867
1868     if (G != NULL) free(G);
1869     if (junk != NULL) free(junk);
1870     return(sim);
1871 }
1872
1873
1874 static int lialg_compute_distance(point_list *input_dompts,
1875                                    point_list *curr_dompts) {
1876     int dist = MAX_DIST;
1877     point_list *A, *B;
1878     int N, M;
1879     int **C = NULL;
1880     int *junk = NULL;
1881     int *BE = NULL;
1882     int *TE = NULL;
1883     int i, j;
1884
1885     /* A is the longer sequence, length N. */
1886     /* B is the shorter sequence, length M. */
1887     if (input_dompts->npts >= curr_dompts->npts) {
1888         A = input_dompts;
1889         N = input_dompts->npts;
1890         B = curr_dompts;
1891         M = curr_dompts->npts;
1892     }
1893     else {
1894         A = curr_dompts;
1895         N = curr_dompts->npts;
1896         B = input_dompts;
1897         M = input_dompts->npts;
1898     }
1899
1900     /* Construct the helper vectors, BE and TE, which say for each column */
1901     /* what are the ``bottom'' and ``top'' rows of interest. */
1902     {
1903         /*      BE = (int *)safe_malloc((N + 1) * sizeof(int));*/
1904         BE = allocate((N + 1), int);
1905         /*      TE = (int *)safe_malloc((N + 1) * sizeof(int)); */
1906         TE = allocate((N + 1), int);
1907
1908         for (j = 1; j <= N; j++) {
1909             int bot, top;
1910
1911             bot = j + (M - DP_BAND);
1912             if (bot > M) bot = M;
1913             BE[j] = bot;
1914
1915             top = j - (N - DP_BAND);
1916             if (top < 1) top = 1;
1917             TE[j] = top;
1918         }
1919     }
1920
1921     /* Allocate and initialize the Cost matrix, C. */
1922     /* The size of C is (M + 1) x (N + 1). */
1923     /* Note that row and column 0 are unused. */
1924     /* Costs are x 100. */
1925     {
1926         /*      C = (int **)safe_malloc((M + 1) * sizeof(int *)); */
1927         C = allocate((M + 1), int *);
1928         /*      junk = (int *)safe_malloc((M + 1) * (N + 1) * sizeof(int)); */
1929         junk = allocate((M + 1) * (N + 1), int);
1930         for (i = 0; i <= M; i++)
1931             C[i] = junk + (i * (N + 1));
1932
1933         for (i = 1; i <= M; i++) {
1934             int bx = B->pts[i-1].x;
1935             int by = B->pts[i-1].y;
1936
1937             for (j = 1; j <= N; j++) {
1938                 int ax = A->pts[j-1].x;
1939                 int ay = A->pts[j-1].y;
1940                 int dx = bx - ax;
1941                 int dy = by - ay;
1942                 int dist = isqrt(10000 * (dx * dx + dy * dy));
1943
1944                 C[i][j] = dist;
1945             }
1946         }
1947     }
1948
1949     /* Do the DP algorithm. */
1950     /* Proceed in column order, from highest column to the lowest. */
1951     /* Within each column, proceed from the highest row to the lowest. */
1952     {
1953         for (j = N; j > 0; j--)
1954             for (i = M; i > 0; i--) {
1955                 int min = MAX_DIST;
1956
1957                 if (i > BE[j] || i < TE[j] || (j == N && i == M))
1958                     continue;
1959
1960                 if (j < N) {
1961                     if (i >= TE[j+1]) {
1962                         int tmp = C[i][j+1];
1963                         if (tmp < min) min = tmp;
1964                     }
1965
1966                     if (i < M) {
1967                         int tmp = C[i+1][j+1];
1968                         if (tmp < min) min = tmp;
1969                     }
1970                 }
1971
1972                 if (i < BE[j]) {
1973                     int tmp = C[i+1][j];
1974                     if (tmp < min) min = tmp;
1975                 }
1976
1977                 C[i][j] += min;
1978             }
1979
1980         dist = (C[1][1] + N / 2) / N;
1981     }
1982
1983     if (C != NULL) free(C);
1984     if (junk != NULL) free(junk);
1985     if (BE != NULL) free(BE);
1986     if (TE != NULL) free(TE);
1987     return(dist);
1988 }
1989
1990
1991 /*************************************************************
1992
1993   Digest-processing routines
1994
1995  *************************************************************/
1996
1997 static int lialg_read_classifier_digest(rClassifier *rec) {
1998     int nclasses;
1999     FILE *fp = NULL;
2000
2001     /* Try to open the corresponding digest file. */
2002     {
2003         char *clx_path;
2004         char *dot;
2005
2006         /* Get a copy of the filename, with some room on the end. */
2007         /*      clx_path = safe_malloc(strlen(rec->file_name) + 5); */
2008         clx_path = allocate(strlen(rec->file_name) + 5, char);
2009         strcpy(clx_path, rec->file_name);
2010
2011         /* Truncate the path after the last dot. */
2012         dot = strrchr(clx_path, '.');
2013         if (dot == NULL) { free(clx_path); return(-1); }
2014         *(dot + 1) = 0;
2015
2016         /* Append the classifier-digest extension. */
2017         strcat(clx_path, "clx");
2018
2019         fp = fopen(clx_path, "r");
2020         if (fp == NULL) { free(clx_path); return(-1); }
2021
2022         free(clx_path);
2023     }
2024
2025     /* Read-in the name and dominant points for each class. */
2026     for (nclasses = 0; !feof(fp); nclasses++) {
2027         point_list *dpts = NULL;
2028         char class[BUFSIZ];
2029         int npts;
2030         int j;
2031
2032         if (fscanf(fp, "%s %d", class, &npts) != 2) {
2033             if (feof(fp)) break;
2034
2035             goto failed;
2036         }
2037         rec->cnames[nclasses] = strdup(class);
2038
2039         /* Allocate a dominant-points list. */
2040         /* dpts = (point_list *)safe_malloc(sizeof(point_list)); */
2041         dpts = allocate(1, point_list);
2042         dpts->pts = make_pen_point_array(npts);
2043         if (dpts->pts == NULL) goto failed;
2044         dpts->npts = npts;
2045         dpts->next = NULL;
2046
2047         /* Read in each point. */
2048         for (j = 0; j < npts; j++) {
2049             int x, y;
2050
2051             if (fscanf(fp, "%d %d", &x, &y) != 2) goto failed;
2052             dpts->pts[j].x = x;
2053             dpts->pts[j].y = y;
2054         }
2055
2056         /* Compute the chain-code. */
2057         lialg_compute_chain_code(dpts);
2058
2059         /* Store the list in the rec data structure. */
2060         rec->dompts[nclasses] = dpts;
2061
2062         continue;
2063
2064 failed:
2065         fprintf(stderr, "read_classifier_digest failed...\n");
2066         for (; nclasses >= 0; nclasses--) {
2067             if (rec->cnames[nclasses] != NULL) {
2068                 free(rec->cnames[nclasses]);
2069                 rec->cnames[nclasses] = NULL;
2070             }
2071             if (rec->dompts[nclasses] != NULL) {
2072                 delete_examples(rec->dompts[nclasses]);
2073                 rec->dompts[nclasses] = NULL;
2074             }
2075         }
2076         if (dpts != NULL)
2077             delete_examples(dpts);
2078         fclose(fp);
2079         return(-1);
2080     }
2081
2082     fclose(fp);
2083     return(0);
2084 }
2085
2086
2087 /*************************************************************
2088
2089   Canonicalization routines
2090
2091  *************************************************************/
2092
2093 static int lialg_canonicalize_examples(rClassifier *rec) {
2094     int i;
2095     int nclasses;
2096
2097     if (lidebug) {
2098         fprintf(stderr, "lialg_canonicalize_examples working on %s\n",
2099                 rec->file_name);
2100     }
2101     /* Initialize canonical-example arrays. */
2102     for (i = 0; i < MAXSCLASSES; i++) {
2103         rec->canonex[i] = NULL;
2104     }
2105
2106     /* Figure out number of classes. */
2107     for (nclasses = 0;
2108           nclasses < MAXSCLASSES && rec->cnames[nclasses] != NULL;
2109           nclasses++)
2110         ;
2111
2112     /* Canonicalize the examples for each class. */
2113     for (i = 0; i < nclasses; i++) {
2114         int j, k;
2115         int nex;
2116         point_list *pts, *tmp, *avg;
2117         int maxxrange, maxyrange;
2118         int minx, miny, maxx, maxy;
2119         int avgxrange, avgyrange, avgxoff, avgyoff, avgscale;
2120
2121
2122         if (lidebug) {
2123             fprintf(stderr, "lialg_canonicalize_examples working on class %s\n",
2124                     rec->cnames[i]);
2125         }
2126         /* Make a copy of the examples. */
2127         pts = NULL;
2128         tmp = rec->ex[i];
2129         for (nex = 0; tmp != NULL; nex++, tmp = tmp->next) {
2130             if ((pts = add_example(pts, tmp->npts, tmp->pts)) == NULL) {
2131                 delete_examples(pts);
2132                 return(-1);
2133             }
2134         }
2135
2136         /* Canonicalize each example, and derive the max x and y ranges. */
2137         maxxrange = 0;
2138         maxyrange = 0;
2139         for (j = 0, tmp = pts; j < nex; j++, tmp = tmp->next) {
2140             if (lialg_canonicalize_example_stroke(tmp) != 0) {
2141                 if (lidebug) {
2142                     fprintf(stderr, "lialg_canonicalize_example_stroke returned error\n");
2143                 }
2144                 return(-1);
2145             }
2146
2147             if (tmp->xrange > maxxrange) maxxrange = tmp->xrange;
2148             if (tmp->yrange > maxyrange) maxyrange = tmp->yrange;
2149         }
2150
2151         /* Normalize max ranges. */
2152         if (((100 * maxxrange + CANONICAL_X / 2) / CANONICAL_X) >
2153             ((100 * maxyrange + CANONICAL_Y / 2) / CANONICAL_Y)) {
2154             maxyrange = (maxyrange * CANONICAL_X + maxxrange / 2) / maxxrange;
2155             maxxrange = CANONICAL_X;
2156         }
2157         else {
2158             maxxrange = (maxxrange * CANONICAL_Y + maxyrange / 2) / maxyrange;
2159             maxyrange = CANONICAL_Y;
2160         }
2161
2162         /* Re-scale each example to max ranges. */
2163         for (j = 0, tmp = pts; j < nex; j++, tmp = tmp->next) {
2164             int scalex = (tmp->xrange == 0) ? 100 : (100 * maxxrange + tmp->xrange / 2) / tmp->xrange;
2165             int scaley = (tmp->yrange == 0) ? 100 : (100 * maxyrange + tmp->yrange / 2) / tmp->yrange;
2166             if (lialg_translate_points(tmp, 0, 0, scalex, scaley) != 0) {
2167                 delete_examples(pts);
2168                 return(-1);
2169             }
2170         }
2171
2172         /* Average the examples; leave average in first example. */
2173         avg = pts;                              /* careful aliasing!! */
2174         for (k = 0; k < NCANONICAL; k++) {
2175             int xsum = 0;
2176             int ysum = 0;
2177
2178             for (j = 0, tmp = pts; j < nex; j++, tmp = tmp->next) {
2179                 xsum += tmp->pts[k].x;
2180                 ysum += tmp->pts[k].y;
2181             }
2182
2183             avg->pts[k].x = (xsum + j / 2) / j;
2184             avg->pts[k].y = (ysum + j / 2) / j;
2185         }
2186
2187         /* Compute BB of averaged stroke and re-scale. */
2188         lialg_get_bounding_box(avg, &minx, &miny, &maxx, &maxy);
2189         avgxrange = maxx - minx;
2190         avgyrange = maxy - miny;
2191         avgscale = (((100 * avgxrange + CANONICAL_X / 2) / CANONICAL_X) >
2192                     ((100 * avgyrange + CANONICAL_Y / 2) / CANONICAL_Y))
2193           ? (100 * CANONICAL_X + avgxrange / 2) / avgxrange
2194           : (100 * CANONICAL_Y + avgyrange / 2) / avgyrange;
2195         if (lialg_translate_points(avg, minx, miny, avgscale, avgscale) != 0) {
2196             delete_examples(pts);
2197             return(-1);
2198         }
2199
2200         /* Re-compute the x and y ranges and center the stroke. */
2201         lialg_get_bounding_box(avg, &minx, &miny, &maxx, &maxy);
2202         avgxrange = maxx - minx;
2203         avgyrange = maxy - miny;
2204         avgxoff = -((CANONICAL_X - avgxrange + 1) / 2);
2205         avgyoff = -((CANONICAL_Y - avgyrange + 1) / 2);
2206         if (lialg_translate_points(avg, avgxoff, avgyoff, 100, 100) != 0) {
2207             delete_examples(pts);
2208             return(-1);
2209         }
2210
2211         /* Create a point list to serve as the ``canonical representation. */
2212         if ((rec->canonex[i] = add_example(NULL, avg->npts, avg->pts)) == NULL) {
2213             delete_examples(pts);
2214             return(-1);
2215         }
2216         (rec->canonex[i])->xrange = maxx - minx;
2217         (rec->canonex[i])->yrange = maxy - miny;
2218
2219         if (lidebug) {
2220             fprintf(stderr, "%s, avgpts = %d\n", rec->cnames[i], avg->npts);
2221             for (j = 0; j < avg->npts; j++) {
2222                 fprintf(stderr, "  (%d, %d)\n",
2223                         avg->pts[j].x, avg->pts[j].y);
2224             }
2225         }
2226
2227         /* Compute dominant points of canonical representation. */
2228         rec->dompts[i] = lialg_compute_dominant_points(avg);
2229
2230         /* Clean up. */
2231         delete_examples(pts);
2232     }
2233
2234     /* Sanity check. */
2235     for (i = 0; i < nclasses; i++) {
2236         char *best_name = lialg_recognize_stroke(rec, rec->canonex[i]);
2237
2238         if (best_name != rec->cnames[i])
2239             fprintf(stderr, "%s, best = %s\n", rec->cnames[i], best_name);
2240     }
2241
2242     return(0);
2243 }
2244
2245
2246 static int lialg_canonicalize_example_stroke(point_list *points) {
2247     int minx, miny, maxx, maxy, xrange, yrange, scale;
2248
2249     /* Filter out points that are too close. */
2250     if (lialg_filter_points(points) != 0) return(-1);
2251
2252     /* Must be at least two points! */
2253     if (points->npts < 2) {
2254         if (lidebug) {
2255             fprintf(stderr, "lialg_canonicalize_example_stroke: npts=%d\n",
2256                     points->npts);
2257         }
2258         return(-1);
2259     }
2260
2261     /* Scale up to avoid conversion errors. */
2262     lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
2263     xrange = maxx - minx;
2264     yrange = maxy - miny;
2265     scale = (((100 * xrange + CANONICAL_X / 2) / CANONICAL_X) >
2266              ((100 * yrange + CANONICAL_Y / 2) / CANONICAL_Y))
2267       ? (100 * CANONICAL_X + xrange / 2) / xrange
2268       : (100 * CANONICAL_Y + yrange / 2) / yrange;
2269     if (lialg_translate_points(points, minx, miny, scale, scale) != 0) {
2270         if (lidebug) {
2271             fprintf(stderr, "lialg_translate_points (minx=%d,miny=%d,scale=%d) returned error\n", minx, miny, scale);
2272         }
2273         return(-1);
2274     }
2275
2276     /* Compute an equivalent stroke with equi-distant points. */
2277     if (lialg_compute_equipoints(points) != 0) return(-1);
2278
2279     /* Re-translate the points to the origin. */
2280     lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
2281     if (lialg_translate_points(points, minx, miny, 100, 100) != 0) {
2282         if (lidebug) {
2283             fprintf(stderr, "lialg_translate_points (minx=%d,miny=%d) returned error\n", minx, miny);
2284         }
2285         return(-1);
2286     }
2287
2288     /* Store the x and y ranges in the point list. */
2289     xrange = maxx - minx;
2290     yrange = maxy - miny;
2291     points->xrange = xrange;
2292     points->yrange = yrange;
2293
2294     if (lidebug) {
2295         int i;
2296         fprintf(stderr, "Canonicalized:   %d, %d, %d, %d\n", minx, miny, maxx, maxy);
2297         for (i = 0; i < points->npts; i++)
2298             fprintf(stderr, "      (%d %d)\n",
2299                     points->pts[i].x, points->pts[i].y);
2300         fflush(stderr);
2301     }
2302
2303     return(0);
2304 }
2305
2306
2307 static int lialg_compute_equipoints(point_list *points) {
2308     pen_point *equipoints = make_pen_point_array(NCANONICAL);
2309     int nequipoints = 0;
2310     int pathlen = lialg_compute_pathlen(points);
2311     int equidist = (pathlen + (NCANONICAL - 1) / 2) / (NCANONICAL - 1);
2312     int i;
2313     int dist_since_last_eqpt;
2314     int remaining_seglen;
2315     int dist_to_next_eqpt;
2316
2317     if (equipoints == NULL) {
2318         error("can't allocate memory in lialg_compute_equipoints");
2319         return(-1);
2320     }
2321
2322     if (lidebug) {
2323         fprintf(stderr, "compute_equipoints:  npts = %d, pathlen = %d, equidist = %d\n",
2324                 points->npts, pathlen, equidist);
2325         fflush(stderr);
2326     }
2327
2328     /* First original point is an equipoint. */
2329     equipoints[0] = points->pts[0];
2330     nequipoints++;
2331     dist_since_last_eqpt = 0;
2332
2333     for (i = 1; i < points->npts; i++) {
2334         int dx1 = points->pts[i].x - points->pts[i-1].x;
2335         int dy1 = points->pts[i].y - points->pts[i-1].y;
2336         int endx = 100 * points->pts[i-1].x;
2337         int endy = 100 * points->pts[i-1].y;
2338         remaining_seglen = isqrt(10000 * (dx1 * dx1 + dy1 * dy1));
2339         dist_to_next_eqpt = equidist - dist_since_last_eqpt;
2340
2341         while (remaining_seglen >= dist_to_next_eqpt) {
2342             if (dx1 == 0) {
2343                 /* x-coordinate stays the same */
2344                 if (dy1 >= 0)
2345                     endy += dist_to_next_eqpt;
2346                 else
2347                     endy -= dist_to_next_eqpt;
2348             }
2349             else {
2350                 int slope = (100 * dy1 + dx1 / 2) / dx1;
2351                 int tmp = isqrt(10000 + slope * slope);
2352                 int dx = (100 * dist_to_next_eqpt + tmp / 2) / tmp;
2353                 int dy = (slope * dx + 50) / 100;
2354
2355                 if (dy < 0) dy = -dy;
2356                 if (dx1 >= 0)
2357                     endx += dx;
2358                 else
2359                     endx -= dx;
2360                 if (dy1 >= 0)
2361                     endy += dy;
2362                 else
2363                     endy -= dy;
2364             }
2365
2366             equipoints[nequipoints].x = (endx + 50) / 100;
2367             equipoints[nequipoints].y = (endy + 50) / 100;
2368             nequipoints++;
2369 /*          assert(nequipoints <= NCANONICAL);*/
2370             dist_since_last_eqpt = 0;
2371             remaining_seglen -= dist_to_next_eqpt;
2372             dist_to_next_eqpt = equidist;
2373         }
2374
2375         dist_since_last_eqpt += remaining_seglen;
2376     }
2377
2378     /* Take care of last equipoint. */
2379     if (nequipoints == NCANONICAL) {
2380         /* Good. */
2381     } else if (nequipoints == (NCANONICAL - 1)) {
2382         /* Make last original point the last equipoint. */
2383         equipoints[nequipoints] = points->pts[points->npts - 1];
2384     } else {
2385       if (lidebug) {
2386         fprintf(stderr,"lialg_compute_equipoints: nequipoints = %d\n",
2387                 nequipoints);
2388       }
2389 /*      assert(false);*/
2390         return(-1);
2391     }
2392
2393     points->npts = NCANONICAL;
2394     delete_pen_point_array(points->pts);
2395     points->pts = equipoints;
2396     return(0);
2397 }
2398
2399
2400 /*************************************************************
2401
2402   Utility routines
2403
2404  *************************************************************/
2405
2406 /* Result is x 100. */
2407 static int lialg_compute_pathlen(point_list *points) {
2408     return(lialg_compute_pathlen_subset(points, 0, points->npts - 1));
2409 }
2410
2411
2412 /* Result is x 100. */
2413 static int lialg_compute_pathlen_subset(point_list *points,
2414                                            int start, int end) {
2415     int pathlen;
2416     int i;
2417
2418     pathlen = 0;
2419     for (i = start + 1; i <= end; i++) {
2420         int dx = points->pts[i].x - points->pts[i-1].x;
2421         int dy = points->pts[i].y - points->pts[i-1].y;
2422         int dist = isqrt(10000 * (dx * dx + dy * dy));
2423         pathlen += dist;
2424     }
2425
2426     return(pathlen);
2427 }
2428
2429
2430 /* Note that this does NOT update points->xrange and points->yrange! */
2431 static int lialg_filter_points(point_list *points) {
2432     int filtered_npts;
2433     pen_point *filtered_pts = make_pen_point_array(points->npts);
2434     int i;
2435
2436     if (filtered_pts == NULL) {
2437         error("can't allocate memory in lialg_filter_points");
2438         return(-1);
2439     }
2440
2441     filtered_pts[0] = points->pts[0];
2442     filtered_npts = 1;
2443     for (i = 1; i < points->npts; i++) {
2444         int j = filtered_npts - 1;
2445         int dx = points->pts[i].x - filtered_pts[j].x;
2446         int dy = points->pts[i].y - filtered_pts[j].y;
2447         int magsq = dx * dx + dy * dy;
2448
2449         if (magsq >= DIST_SQ_THRESHOLD) {
2450             filtered_pts[filtered_npts] = points->pts[i];
2451             filtered_npts++;
2452         }
2453     }
2454
2455     points->npts = filtered_npts;
2456     delete_pen_point_array(points->pts);
2457     points->pts = filtered_pts;
2458     return(0);
2459 }
2460
2461
2462 /* scalex and scaley are x 100. */
2463 /* Note that this does NOT update points->xrange and points->yrange! */
2464 static int lialg_translate_points(point_list *points,
2465                                    int minx, int miny,
2466                                    int scalex, int scaley) {
2467     int i;
2468
2469     for (i = 0; i < points->npts; i++) {
2470         points->pts[i].x = ((points->pts[i].x - minx) * scalex + 50) / 100;
2471         points->pts[i].y = ((points->pts[i].y - miny) * scaley + 50) / 100;
2472     }
2473
2474     return(0);
2475 }
2476
2477
2478 static void lialg_get_bounding_box(point_list *points,
2479                                     int *pminx, int *pminy,
2480                                     int *pmaxx, int *pmaxy) {
2481     int minx, miny, maxx, maxy;
2482     int i;
2483
2484     minx = maxx = points->pts[0].x;
2485     miny = maxy = points->pts[0].y;
2486     for (i = 1; i < points->npts; i++) {
2487         pen_point *pt = &(points->pts[i]);
2488         if (pt->x < minx) minx = pt->x;
2489         if (pt->x > maxx) maxx = pt->x;
2490         if (pt->y < miny) miny = pt->y;
2491         if (pt->y > maxy) maxy = pt->y;
2492     }
2493
2494     *pminx = minx;
2495     *pminy = miny;
2496     *pmaxx = maxx;
2497     *pmaxy = maxy;
2498 }
2499
2500 #ifdef __ECOS
2501 float
2502 expf(float x)
2503 {
2504     return exp((double)x);
2505 }
2506 #endif
2507
2508
2509 static void lialg_compute_lpf_parameters() {
2510     int i;
2511
2512     for (i = LP_FILTER_WIDTH; i >= 0; i--) {
2513         float x = 0.04 * (i * i);
2514 #if defined(ARM_LINUX) || !defined(__GLIBC__)
2515         double tmp = 100.0 * exp((double)x);
2516 #else
2517         float tmp = 100.0 * expf(x);
2518 #endif
2519         int wt = rint((double)tmp);
2520
2521         lialg_lpfwts[LP_FILTER_WIDTH - i] = wt;
2522         lialg_lpfwts[LP_FILTER_WIDTH + i] = wt;
2523     }
2524     lialg_lpfconst = 0;
2525     for (i = 0; i < (2 * LP_FILTER_WIDTH + 1); i++) {
2526         lialg_lpfconst += lialg_lpfwts[i];
2527     }
2528 }
2529
2530
2531 /* Code from Joseph Hall (jnhall@sat.mot.com). */
2532 static int isqrt(int n) {
2533     register int i;
2534     register long k0, k1, nn;
2535
2536     for (nn = i = n, k0 = 2; i > 0; i >>= 2, k0 <<= 1)
2537         ;
2538     nn <<= 2;
2539     for (;;) {
2540         k1 = (nn / k0 + k0) >> 1;
2541         if (((k0 ^ k1) & ~1) == 0)
2542             break;
2543         k0 = k1;
2544     }
2545     return (int) ((k1 + 1) >> 1);
2546 }
2547
2548
2549 /* Helper routines from Mark Hayter. */
2550 static int likeatan(int tantop, int tanbot) {
2551     int t;
2552     /* Use tan(theta)=top/bot --> order for t */
2553     /* t in range 0..0x40000 */
2554
2555     if ((tantop == 0) && (tanbot == 0))
2556         t = 0;
2557     else
2558     {
2559         t = (tantop << 16) / (abs(tantop) + abs(tanbot));
2560         if (tanbot < 0)
2561             t = 0x20000 - t;
2562         else
2563             if (tantop < 0) t = 0x40000 + t;
2564     }
2565     return t;
2566 }
2567
2568
2569 static int quadr(int t) {
2570     return (8 - (((t + 0x4000) >> 15) & 7)) & 7;
2571 }