Subversion Repositories planix.SVN

Rev

Blame | Last modification | View Log | RSS feed

/*
 *  li_recognizer.c
 *
 *      Copyright 2000 Compaq Computer Corporation.
 *      Copying or modifying this code for any purpose is permitted,
 *      provided that this copyright notice is preserved in its entirety
 *      in all copies or modifications.
 *      COMPAQ COMPUTER CORPORATION MAKES NO WARRANTIES, EXPRESSED OR
 *      IMPLIED, AS TO THE USEFULNESS OR CORRECTNESS OF THIS CODE OR
 *
 *
 *  Adapted from cmu_recognizer.c by Jay Kistler.
 *  
 *  Where is the CMU copyright???? Gotta track it down - Jim Gettys
 *
 *  Credit to Dean Rubine, Jim Kempf, and Ari Rapkin.
 */

#include <u.h>
#include <libc.h>
#include <stdio.h>
#include <draw.h>
#include <scribble.h>
#include "scribbleimpl.h"

#include "hre_internal.h"
#include "li_recognizer_internal.h"

int lidebug = 0;

#define LI_MAGIC 0xACCBADDD

#define CHECK_LI_MAGIC(_a) \
  ((_a) != nil && ((li_recognizer*)(_a))->li_magic == LI_MAGIC)


static void lialg_initialize(rClassifier *);
static int lialg_read_classifier_digest(rClassifier *);
static int lialg_canonicalize_examples(rClassifier *);
static char *lialg_recognize_stroke(rClassifier *, point_list *);
static void lialg_compute_lpf_parameters(void);


char* li_err_msg = nil;

#define bcopy(s1,s2,n) memmove(s2,s1,n)

/*Freeing classifier*/

static void
free_rClassifier(rClassifier* rc);

/*
 * Point List Support
*/

static point_list*
add_example(point_list* l,int npts,pen_point* pts)
{
        pen_point* lpts = mallocz(npts*sizeof(pen_point), 1);
        point_list *p = malloc(sizeof(point_list));

        p->npts = npts;
        p->pts = lpts;
        p->next = l;       /*Order doesn't matter, so we stick on end.*/

        /*Copy points.*/

        bcopy(pts, lpts, npts * sizeof(pen_point));

        return(p);
}
        

static void
delete_examples(point_list* l)
{
        point_list* p;

        for( ; l != nil; l = p ) {
                p = l->next;
                free(l->pts);
                free(l);
        }
}

/*
 * Local functions
 */

/*
 * recognize_internal-Form Vector, use Classifier to classify, return char.
 */

static char*
recognize_internal(rClassifier* rec, Stroke* str, int*)
{
        char *res;
        point_list *stroke;

        stroke = add_example(nil, str->npts, str->pts);
        if (stroke == nil) return(nil);

        res = lialg_recognize_stroke(rec, stroke);

        delete_examples(stroke);
        return(res);
}

/*
 * file_path-Construct pathname, check for proper extension.
 */

static int
  file_path(char* dir,char* filename,char* pathname)
{
        char* dot;
        
        /*Check for proper extension on file name.*/
        
        dot = strrchr(filename,'.');
        
        if( dot == nil ) {
                return(-1);
        }

        /*Determine whether a gesture or character classifier.*/

        if( strcmp(dot,LI_CLASSIFIER_EXTENSION) != 0 ) {
                return(-1);
        }

        /*Concatenate directory and filename into pathname.*/
        
        strcpy(pathname,dir);
        strcat(pathname,"/");
        strcat(pathname,filename);
        
        return(0);
}

/*read_classifier_points-Read points so classifier can be extended.*/

static int 
read_classifier_points(FILE* fd,int nclss,point_list** ex,char** cnames)
{
        int i,j,k;
        char buf[BUFSIZ];
        int nex = 0;
        char* names[MAXSCLASSES];
        point_list* examples[MAXSCLASSES];
        pen_point* pts;
        int npts;

        /*Initialize*/

                for( i = 0; i < MAXSCLASSES; i++ ) {
                                names[i] = nil;
                                examples[i] = nil;
                }

                /*Go thru classes.*/

                for( k = 0; k < nclss; k++ ) {

                                /*Read class name and number of examples.*/
                
                                if( fscanf(fd,"%d %s",&nex,buf) != 2 )
                                        goto unallocate;
                
                                /*Save class name*/
                
                                names[k] = strdup(buf);
                
                                /*Read examples.*/
                                
                                for( i = 0; i < nex; i++ ) {
                                        
                                        /*Read number of points.*/
                                        
                                        if( fscanf(fd,"%d",&npts) != 1 )
                                                                goto unallocate; /*Boy would I like exceptions!*/
                                        
                                        /*Allocate array for points.*/
                                        
                                        if( (pts = mallocz(npts*sizeof(pen_point), 1)) == nil )
                                                                goto unallocate;
                                        
                                        /*Read in points.*/
                                        
                                        for( j = 0; j < npts; j++ ) {
                                                                int x,y;
                                                                if( fscanf(fd,"%d %d",&x,&y) != 2 ) {
                                                                        delete_pen_point_array(pts);
                                                                        goto unallocate;
                                                                }
                                                                pts[j].Point = Pt(x, y);
                                        }
                                        
                                        /*Add example*/
                                        
                                        if( (examples[k] = add_example(examples[k],npts,pts)) == nil ) {
                                                                delete_pen_point_array(pts);
                                                                goto unallocate;
                                        }
                                        
                                        delete_pen_point_array(pts);
                        
                                }
                }

/* ari -- end of list of classes */
/* fprint(2, "]\n"); */

                /*Transfer to recognizer.*/

                bcopy(examples,ex,sizeof(examples));
                bcopy(names,cnames,sizeof(names));

                return(0);

                /*Error. Deallocate memory and return.*/

unallocate:
                for( ; k >= 0; k-- ) {
                                delete_examples(examples[k]);
                                free(names[k]);
                }

                return(-1);
}

/*read_classifier-Read a classifier file.*/

static int read_classifier(FILE* fd,rClassifier* rc)
{
        
        li_err_msg = nil;

        /*Read in classifier file.*/

        if(fscanf(fd, "%d", &rc->nclasses) != 1)
                return -1;

        /*Read in the example points, so classifier can be extended.*/

        if( read_classifier_points(fd,rc->nclasses,rc->ex,rc->cnames) != 0 )
                return -1;

        return(0);
}

/*
 * Extension Functions
*/

/* getClasses and clearState are by Ari */

static int
recognizer_getClasses (recognizer r, char ***list, int *nc)
{
        int i, nclasses;
        li_recognizer* rec;
        char **ret;

                rec = (li_recognizer*)r->recognizer_specific;

                /*Check for LI recognizer.*/

                if( !CHECK_LI_MAGIC(rec) ) {
                                li_err_msg = "Not a LI recognizer";
                                return(-1);
        }
        
        *nc = nclasses = rec->li_rc.nclasses;
        ret = malloc(nclasses*sizeof(char*));

        for (i = 0; i < nclasses; i++)
                                ret[i] = rec->li_rc.cnames[i];   /* only the 1st char of the cname */
        *list = ret;
                return 0;
}

static int
recognizer_clearState (recognizer)
{
  /*This operation isn't supported by the LI recognizer.*/

  li_err_msg = "Clearing state is not supported by the LI recognizer";

  return(-1);
}

static bool isa_li(recognizer r)
{ return(CHECK_LI_MAGIC(r)); }

static int
recognizer_train(recognizer, rc*, uint, Stroke*, rec_element*, bool)
{
  /*This operation isn't supported by the LI recognizer.*/

  li_err_msg = "Training is not supported by the LI recognizer";

  return(-1);
}

int
li_recognizer_get_example (recognizer r,
                                                   int          class, 
                                                   int          instance,
                                                   char         **name, 
                                                   pen_point    **points,
                                                   int          *npts)
{
        li_recognizer   *rec = (li_recognizer*)r->recognizer_specific;
        int nclasses = rec->li_rc.nclasses;
        point_list          *pl;
        
                if( !CHECK_LI_MAGIC(rec) ) {
                                li_err_msg = "Not a LI recognizer";
                                return(-1);
                }
                if (class > nclasses)
                                return -1;
                pl = rec->li_rc.canonex[class];
                while (instance && pl)
                {
                                pl = pl->next;
                                instance--;
                }
                if (!pl)
                                return -1;
                *name = rec->li_rc.cnames[class];
                *points = pl->pts;
                *npts = pl->npts;
                return pl->npts; /* I hope [sjm] */
}

/*
 * API Functions
 */


/*li_recognizer_load-Load a classifier file.*/

static int li_recognizer_load(recognizer r, char* dir, char* filename)
{
                FILE *fd;
                char* pathname;
                li_recognizer* rec;
                rClassifier* rc;
                
                rec = (li_recognizer*)r->recognizer_specific;

                /*Make sure recognizer's OK*/

                if( !CHECK_LI_MAGIC(rec) ) {
                                li_err_msg = "Not a LI recognizer";
                                return(-1);
                }

                rc = &(rec->li_rc);

                /*Check parameters.*/

                if( filename == nil ) {
                                li_err_msg = "Invalid parameters";
                                return(-1);
                }

                /*We let the directory be null.*/

                if( dir == nil || (int)strlen(dir) <= 0 ) {
                                dir = ".";
                }

                if(0)fprint(2, "dir = %s, filename = %s\n",
                                dir, filename);

                /*Make full pathname and check filename*/

                pathname = malloc((strlen(dir) + strlen(filename) + 2)*sizeof(char));
                if( file_path(dir, filename, pathname) == -1 ) {
                                free(pathname);
                                li_err_msg = "Not a LI recognizer classifier file";
                                return -1;
                }

                /* Try to short-circuit the full classifier-file processing. */
                rc->file_name = pathname;
                if (lialg_read_classifier_digest(rc) == 0)
                                return(0);
                rc->file_name = nil;

                /*Open the file*/

                if( (fd = fopen(pathname,"r")) == nil ) {
                                li_err_msg = "Can't open classifier file";
                                if(0)fprint(2, "Can't open %s.\n", pathname);
                                free(pathname);
                                return(-1);
                }

                /*If rClassifier is OK, then delete it first.*/

                if( rc->file_name != nil ) {
                                free_rClassifier(rc);
                }

                /*Read classifier.*/
                
                if( read_classifier(fd,rc) < 0 ) {
                                free(pathname);
                                return(-1);
                }

                /*Close file.*/

                fclose(fd);

                /*Add classifier name.*/

                rc->file_name = pathname;

                /* Canonicalize examples. */
                if (lialg_canonicalize_examples(rc) != 0) {
                                free(pathname);
                                rc->file_name = nil;
                                return -1;
                }

                return(0);
}

/*li_recognizer_save-Save a classifier file.*/

static int li_recognizer_save(recognizer, char*, char*)
{ 
                /*This operation isn't supported by the LI recognizer.*/

                li_err_msg = "Saving is not supported by the LI recognizer";
                return -1;
}

static wordset
li_recognizer_load_dictionary(recognizer, char*, char*)
{
                /*This operation isn't supported by the LI recognizer.*/

                li_err_msg = "Dictionaries are not supported by the LI recognizer";
                return nil;
}

static int
li_recognizer_save_dictionary(recognizer, char*, char*, wordset)
{
                /*This operation isn't supported by the LI recognizer.*/
                li_err_msg = "Dictionaries are not supported by the LI recognizer";
                return -1;
}

static int
li_recognizer_free_dictionary(recognizer, wordset)
{
  /*This operation isn't supported by the LI recognizer.*/

  li_err_msg = "Dictionaries are not supported by the LI recognizer";

  return -1;

}

static int
li_recognizer_add_to_dictionary(recognizer, letterset*, wordset)
{
                /*This operation isn't supported by the LI recognizer.*/
                li_err_msg = "Dictionaries are not supported by the LI recognizer";
                return -1;
}

static int
li_recognizer_delete_from_dictionary(recognizer, letterset*, wordset)
{
                /*This operation isn't supported by the LI recognizer.*/
                li_err_msg = "Dictionaries are not supported by the LI recognizer";
                return -1;
}

static char*
li_recognizer_error(recognizer rec)
{
        char* ret = li_err_msg;

        /*Check for LI recognizer.*/

        if( !CHECK_LI_MAGIC(rec->recognizer_specific) ) {
                                li_err_msg = "Not a LI recognizer";
                                return nil;
        }
        li_err_msg = nil;
        return ret;
}

static int 
li_recognizer_clear(recognizer r, bool)
{
                li_recognizer* rec; 

                rec = (li_recognizer*)r->recognizer_specific;
                /*Check for LI recognizer.*/
                if( !CHECK_LI_MAGIC(rec) ) {
                                li_err_msg = "Not a LI recognizer";
                                return 0;
                }
                return 0;
}

static int 
li_recognizer_set_context(recognizer, rc*)
{
                /*This operation isn't supported by the LI recognizer.*/
                li_err_msg = "Contexts are not supported by the LI recognizer";
                return -1;
}

static rc*
li_recognizer_get_context(recognizer)
{
                /*This operation isn't supported by the LI recognizer.*/
                li_err_msg = "Contexts are not supported by the LI recognizer";
                return nil;
}

static int 
li_recognizer_get_buffer(recognizer, uint*, Stroke**)
{
                /*This operation isn't supported by the LI recognizer.*/
                li_err_msg = "Buffer get/set are not supported by the LI recognizer";
                return -1;
}

static int 
li_recognizer_set_buffer(recognizer, uint, Stroke*)
{
                /*This operation isn't supported by the LI recognizer.*/
                li_err_msg = "Buffer get/set are not supported by the LI recognizer";
                return -1;
}

static int
li_recognizer_translate(recognizer r, uint ncs, Stroke* tps, bool, int* nret, rec_alternative** ret)
{
        char* clss;
        li_recognizer* rec; 
        int conf;
        rClassifier* rc;
          
        rec = (li_recognizer*)r->recognizer_specific;

        *nret = 0;
        *ret = nil;

        /*Check for LI recognizer.*/

        if( !CHECK_LI_MAGIC(rec) ) {
                                li_err_msg = "Not a LI recognizer";
                                return(-1);
        }

                rc = &(rec->li_rc);

        /*Check for valid parameters.*/
        if (ncs < 1) {
                                li_err_msg = "Invalid parameters: ncs";
                                return(-1);
        }
        if( tps == nil) {
                                li_err_msg = "Invalid parameters: tps";
                                return(-1);
        }
        if( nret == nil) {
                                li_err_msg = "Invalid parameters: nret";
                                return(-1);
        }
        if( ret == nil) {
                                li_err_msg = "Invalid parameters: ret";
                                return(-1);
        }

        /*
         * Go through the stroke array and recognize. Since this is a single
         *   stroke recognizer, each stroke is treated as a separate
         *   character or gesture. We allow only characters or gestures
         *   to be recognized at one time, since otherwise, handling
         *   the display of segmentation would be difficult. 
        */
        clss = recognize_internal(rc,tps,&conf);
        if (clss == nil) {
                                *nret = 1;
                                return(0);
        }

        /*Return values.*/
        *nret = 1;
        return(*clss);
}


static rec_fn*
li_recognizer_get_extension_functions(recognizer rec)
{
        rec_fn* ret;

        /*Check for LI recognizer.*/

        if( !CHECK_LI_MAGIC(rec->recognizer_specific) ) {
                                li_err_msg = "Not a LI recognizer";
                                return(nil);
        }

        ret = make_rec_fn_array(LI_NUM_EX_FNS);

/* ari -- clearState & getClasses are mine */
        ret[LI_GET_CLASSES] =   (rec_fn)recognizer_getClasses;
        ret[LI_CLEAR] =                 (rec_fn)recognizer_clearState;
        ret[LI_ISA_LI] =                (rec_fn)isa_li;
        ret[LI_TRAIN] =                 (rec_fn)recognizer_train;

        return(ret);
}

static char**
li_recognizer_get_gesture_names(recognizer)
{
                /*This operation isn't supported by the LI recognizer.*/
                li_err_msg = "Gestures are not supported by the LI recognizer";
                return nil;
}

static xgesture
li_recognizer_set_gesture_action(recognizer, char*, xgesture, void*)
{
                /*This operation isn't supported by the LI recognizer.*/
                li_err_msg = "Gestures are not supported by the LI recognizer";
                return nil;
}

/*
 * Exported Functions
*/

/*RECOGNIZER_INITIALIZE-Initialize the recognizer.*/

/* note from ari:  this expands via pre-processor to
 *
 * recognizer __recognizer_internal_initialize(rec_info* ri)
 */

RECOGNIZER_INITIALIZE(ri)
{
        recognizer r;
        li_recognizer* rec;
        int i;

        /*Check that locale matches.*/

        if( strcmp(ri->ri_locale,LI_SUPPORTED_LOCALE) != 0 ) {
                li_err_msg = "Not a supported locale";
/* fprint(2, "Locale error.\n");*/
                return(nil);
        }

        /*
         * Check that character sets match. Note that this is only approximate,
         * since the classifier file will have more information.
        */

        if( ri->ri_subset != nil ) {
          for(i = 0; ri->ri_subset[i] != nil; i++ ) {

                if( strcmp(ri->ri_subset[i],UPPERCASE) != 0 &&
                        strcmp(ri->ri_subset[i],LOWERCASE) != 0  &&
                        strcmp(ri->ri_subset[i],DIGITS) != 0 &&
                        strcmp(ri->ri_subset[i],GESTURE) != 0 ) {
                  li_err_msg = "Not a supported character set";
/* fprint(2, "charset error.\n"); */

                  return(nil);
                }
          }
        }
                         
/* ari */
        r = make_recognizer(ri);
/* fprint(2, "past make_recognizer.\n"); */

        if( r == nil ) {
                li_err_msg = "Can't allocate storage";

                return(nil);
        }

        /*Make a LI recognizer structure.*/


        /* rec = (li_recognizer*)safe_malloc(sizeof(li_recognizer))) == nil ); */

        rec = malloc(sizeof(li_recognizer));

        r->recognizer_specific = rec;

        rec->li_rc.file_name = nil;
        rec->li_rc.nclasses = 0;

        /*Initialize the recognizer struct.*/

        r->recognizer_load_state = li_recognizer_load;
        r->recognizer_save_state = li_recognizer_save;
        r->recognizer_load_dictionary = li_recognizer_load_dictionary;
        r->recognizer_save_dictionary = li_recognizer_save_dictionary;
        r->recognizer_free_dictionary = li_recognizer_free_dictionary;
        r->recognizer_add_to_dictionary = li_recognizer_add_to_dictionary;
        r->recognizer_delete_from_dictionary = li_recognizer_delete_from_dictionary;
        r->recognizer_error = li_recognizer_error;
        r->recognizer_translate = li_recognizer_translate;
        r->recognizer_get_context = li_recognizer_get_context;
        r->recognizer_set_context = li_recognizer_set_context;
        r->recognizer_get_buffer = li_recognizer_get_buffer;
        r->recognizer_set_buffer = li_recognizer_set_buffer;
        r->recognizer_clear = li_recognizer_clear;
        r->recognizer_get_extension_functions = 
          li_recognizer_get_extension_functions;
        r->recognizer_get_gesture_names = li_recognizer_get_gesture_names;
        r->recognizer_set_gesture_action = 
          li_recognizer_set_gesture_action;

        /*Initialize LI Magic Number.*/

        rec->li_magic = LI_MAGIC;

        /*Initialize rClassifier.*/

        rec->li_rc.file_name = nil;

        for( i = 0; i < MAXSCLASSES; i++ ) {
                rec->li_rc.ex[i] = nil;
                rec->li_rc.cnames[i] = nil;
        }

        lialg_initialize(&rec->li_rc);

        /*Get rid of error message. Not needed here.*/
        li_err_msg = nil;

        return(r);
}

/*free_rClassifier-Free the rClassifier.*/

static void
free_rClassifier(rClassifier* rc)
{
        int i;

        if( rc->file_name != nil) {
                free(rc->file_name);
        }

        for( i = 0; rc->ex[i] != nil; i++) {
                delete_examples(rc->ex[i]);
                free(rc->cnames[i]);
        }

}

/*RECOGNIZER_FINALIZE-Deallocate the recognizer, finalize.*/

RECOGNIZER_FINALIZE(r)
{
                li_recognizer* rec = (li_recognizer*)r->recognizer_specific;

                /*Make sure this is a li_recognizer first*/
                if( !CHECK_LI_MAGIC(rec) ) {
                                li_err_msg = "Not a LI recognizer";
                                return(-1);
        }

                /*Deallocate rClassifier resources.*/

                free_rClassifier(&(rec->li_rc));

                /*Deallocate the li_recognizer struct.*/
                free(rec);

                /*Deallocate the recognizer*/
                delete_recognizer(r);

                return(0);
}


/* **************************************************

  Implementation of the Li/Yeung recognition algorithm

************************************************** */

#define WORST_SCORE     0x7fffffff

/* Dynamic programming parameters */
#define DP_BAND         3
#define MIN_SIM         0
#define MAX_DIST        0x7fffffff
#define SIM_THLD        60      /* x 100 */
#define DIST_THLD       3200    /* x 100 */

/* Low-pass filter parameters -- empirically derived */
#define LP_FILTER_WIDTH 6
#define LP_FILTER_ITERS 8
#define LP_FILTER_THLD  250     /* x 100 */
#define LP_FILTER_MIN   5

/* Pseudo-extrema parameters -- empirically derived */
#define PE_AL_THLD      1500    /* x 100 */
#define PE_ATCR_THLD    135     /* x 100 */

/* Contour-angle derivation parameters */
#define T_ONE           1
#define T_TWO           20

/* Pre-processing and canonicalization parameters */
#define CANONICAL_X     108
#define CANONICAL_Y     128
#define DIST_SQ_THRESHOLD   (3*3)       /* copied from fv.h */
#define NCANONICAL      50

/* Tap-handling parameters */
#define TAP_CHAR        "."
#define TAP_TIME_THLD   150         /* msec */
#define TAP_DIST_THLD   75          /* dx * dx + dy * dy */
#define TAP_PATHLEN     1000        /* x 100 */


/* region types */
#define RGN_CONVEX  0
#define RGN_CONCAVE 1
#define RGN_PLAIN   2
#define RGN_PSEUDO  3


typedef struct RegionList {
        int start;
        int end;
        int type;
        struct RegionList *next;
} region_list;


/* direction-code table; indexed by dx, dy */
static int lialg_dctbl[3][3] = {{1, 0, 7}, {2, 0x7FFFFFFF, 6}, {3, 4, 5}};

/* low-pass filter weights */
static int lialg_lpfwts[2 * LP_FILTER_WIDTH + 1];
static int lialg_lpfconst = -1;


static int lialg_preprocess_stroke(point_list *);
static point_list *lialg_compute_dominant_points(point_list *);
static point_list *lialg_interpolate_points(point_list *);
static void lialg_bresline(pen_point *, pen_point *, point_list *, int *);
static void lialg_compute_chain_code(point_list *);
static void lialg_compute_unit_chain_code(point_list *);
static region_list *lialg_compute_regions(point_list *);
static point_list *lialg_compute_dompts(point_list *, region_list *);
static int *lialg_compute_contour_angle_set(point_list *, region_list *);
static void lialg_score_stroke(point_list *, point_list *, int *, int *);
static int lialg_compute_similarity(point_list *, point_list *);
static int lialg_compute_distance(point_list *, point_list *);

static int lialg_read_classifier_digest(rClassifier *);

static int lialg_canonicalize_examples(rClassifier *);
static int lialg_canonicalize_example_stroke(point_list *);
static int lialg_compute_equipoints(point_list *);

static int lialg_compute_pathlen(point_list *);
static int lialg_compute_pathlen_subset(point_list *, int, int);
static int lialg_filter_points(point_list *);
static int lialg_translate_points(point_list *, int, int, int, int);
static void lialg_get_bounding_box(point_list *, int *, int *, int *, int *);
static void lialg_compute_lpf_parameters();
static int isqrt(int);
static int likeatan(int, int);
static int quadr(int);


/*************************************************************

  Core routines for the Li/Yeung recognition algorithm

 *************************************************************/

static void lialg_initialize(rClassifier *rec) {
        int i;

        /* Initialize the dompts arrays. */
        for (i = 0; i < MAXSCLASSES; i++) {
                rec->dompts[i] = nil;
        }
}


/*
 *  Main recognition routine -- called by HRE API.
 */
static char *lialg_recognize_stroke(rClassifier *rec, point_list *stroke) {
                int i;
                char *name = nil;
                point_list *input_dompts = nil;
                char *best_name = nil;
                int best_score = WORST_SCORE;
                char *curr_name;
                point_list *curr_dompts;

                /*    (void)gettimeofday(&stv, nil);*/

                if (stroke->npts < 1) goto done;

                /* Check for tap. */

                /* First thing is to filter out ``close points.'' */
                if (lialg_filter_points(stroke) != 0) return(nil);

                /* Unfortunately, we don't have the actual time that each point */
                /* was recorded (i.e., dt is invalid).  Hence, we have to use a */
                /* heuristic based on total distance and the number of points. */
                if (stroke->npts == 1 || lialg_compute_pathlen(stroke) < TAP_PATHLEN)
                        return(TAP_CHAR);

                /* Pre-process input stroke. */
                if (lialg_preprocess_stroke(stroke) != 0) goto done;

                /* Compute its dominant points. */
                input_dompts = lialg_compute_dominant_points(stroke);
                if (input_dompts == nil) goto done;

                /* Score input stroke against every class in classifier. */
                for (i = 0, curr_name = rec->cnames[i], curr_dompts = rec->dompts[i];
                                i < MAXSCLASSES && curr_name != nil && curr_dompts != nil;
                                i++, curr_name = rec->cnames[i], curr_dompts = rec->dompts[i]) {
                                int sim;
                                int dist;
                                int curr_score;

                                lialg_score_stroke(input_dompts, curr_dompts, &sim, &dist);
                                curr_score = dist;

                                if (lidebug && curr_score < DIST_THLD)
                                fprint(2, "(%s, %d, %d) ", curr_name, sim, dist);

                                /* Is it the best so far? */
                                if (curr_score < best_score && curr_score <= DIST_THLD) {
                                best_score = curr_score;
                                best_name = curr_name;
                                }
                }

                if (lidebug)
                                fprint(2, "\n");

                /* No errors. */
                name = best_name;

done:
                delete_examples(input_dompts);
                return(name);
}


static int lialg_preprocess_stroke(point_list *points) {
        int minx, miny, maxx, maxy, xrange, yrange, scale, xoff, yoff;

        /* Filter out points that are too close. */
        /* We did this earlier, when we checked for a tap. */
/*
        if (lialg_filter_points(points) != 0) return(-1);
*/

/*    assert(points->npts > 0);*/

        /* Scale up to avoid conversion errors. */
        lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
        xrange = maxx - minx;
        yrange = maxy - miny;
        scale = ( ((100 * xrange + CANONICAL_X / 2) / CANONICAL_X) > 
                          ((100 * yrange + CANONICAL_Y / 2) / CANONICAL_Y))
          ? (100 * CANONICAL_X + xrange / 2) / xrange
          : (100 * CANONICAL_Y + yrange / 2) / yrange;
        if (lialg_translate_points(points, minx, miny, scale, scale) != 0)
                return(-1);

        /* Center the stroke. */
        lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
        xrange = maxx - minx;
        yrange = maxy - miny;
        xoff = -((CANONICAL_X - xrange + 1) / 2);
        yoff = -((CANONICAL_Y - yrange + 1) / 2);
        if (lialg_translate_points(points, xoff, yoff, 100, 100) != 0) return(-1);

        /* Store the x and y ranges in the point list. */
        xrange = maxx - minx;
        yrange = maxy - miny;
        points->xrange = xrange;
        points->yrange = yrange;

        if (lidebug) {
                int i;
                fprint(2, "After pre-processing:   %d %d %d %d\n",
                                minx, miny, maxx, maxy);
                for (i = 0; i < points->npts; i++)
                        fprint(2, "      (%P)\n", points->pts[i].Point);
                fflush(stderr);
        }

        return(0);
}


static point_list *lialg_compute_dominant_points(point_list *points) {
        point_list *ipts;
        region_list *regions;
        point_list *dpts;

        /* Interpolate points. */
        ipts = lialg_interpolate_points(points);
        if (ipts == nil) return(nil);
        if (lidebug) {
                                int j;
                                fprint(2, "After interpolation:  %d ipts\n", ipts->npts);
                                for (j = 0; j < ipts->npts; j++) {
                                        fprint(2, "  (%P), %lud\n", ipts->pts[j].Point, ipts->pts[j].chaincode);
                                }
                                fflush(stderr);
        }

        /* Compute regions. */
        regions = lialg_compute_regions(ipts);
/*    assert(regions != nil);*/

        /* Compute dominant points. */
        dpts = lialg_compute_dompts(ipts, regions);
        if (lidebug) {
                                int j;
                                fprint(2, "Dominant points:  ");
                                for (j = 0; j < dpts->npts; j++) {
                                        fprint(2, "%P (%lud)  ", dpts->pts[j].Point, dpts->pts[j].chaincode);
                                }
                                fprint(2, "\n");
                                fflush(stderr);
        }

        /* Delete region data structure. */
        {
                                region_list *curr, *next;
                                for (curr = regions; curr != nil; ) {
                                        next = curr->next;
                                        free(curr);
                                        curr = next;
                                }
        }
        delete_examples(ipts);
        return(dpts);
}

/* Input points are assumed to be integer-valued! */
static point_list *lialg_interpolate_points(point_list *points) {
        int i, j;
        int maxpts;
        point_list *newpts;

        /* Compute an upper-bound on the number of interpolated points. */
        maxpts = 0;
        for (i = 0; i < (points->npts - 1); i++) {
                                pen_point *pta = &(points->pts[i]);
                                pen_point *ptb = &(points->pts[i+1]);
                                maxpts += abs(pta->x - ptb->x) + abs(pta->y - ptb->y);
        }

        /* Allocate an array of the requisite size. */
        maxpts += points->npts;
        /* newpts = (point_list *)safe_malloc(sizeof(point_list)); */
        newpts = malloc(sizeof(point_list));
        newpts->pts = mallocz(maxpts*sizeof(pen_point), 1);
        if (newpts->pts == nil) {
                                free(newpts);
                                return(nil);
        }
        newpts->npts = maxpts;
        newpts->next = nil;

        /* Interpolate each of the segments. */
        j = 0;
        for (i = 0; i < (points->npts - 1); i++) {
                                pen_point *startpt = &(points->pts[i]);
                                pen_point *endpt = &(points->pts[i+1]);

                                lialg_bresline(startpt, endpt, newpts, &j);

                                j--;    /* end point gets recorded as start point of next segment! */
        }

        /* Add-in last point. */
        newpts->pts[j++] = points->pts[points->npts - 1];
        newpts->npts = j;

        /* Compute the chain code for P (the list of points). */
        lialg_compute_unit_chain_code(newpts);

        return(newpts);
}


/* This implementation is due to Kenny Hoff. */
static void lialg_bresline(pen_point *startpt, pen_point *endpt,
                                                        point_list *newpts, int *j) {
        int Ax, Ay, Bx, By, dX, dY, Xincr, Yincr;

        Ax = startpt->x;
        Ay = startpt->y;
        Bx = endpt->x;
        By = endpt->y;

        /* INITIALIZE THE COMPONENTS OF THE ALGORITHM THAT ARE NOT AFFECTED */
        /* BY THE SLOPE OR DIRECTION OF THE     LINE */
        dX = abs(Bx-Ax);        /* store the change in X and Y of the line endpoints */
        dY = abs(By-Ay);

        /* DETERMINE "DIRECTIONS" TO INCREMENT X AND Y (REGARDLESS OF DECISION) */
        if (Ax > Bx) { Xincr=-1; } else { Xincr=1; }    /* which direction in X? */
        if (Ay > By) { Yincr=-1; } else { Yincr=1; }    /* which direction in Y? */

        /* DETERMINE INDEPENDENT VARIABLE (ONE THAT ALWAYS INCREMENTS BY 1 (OR -1) ) */
        /* AND INITIATE APPROPRIATE LINE DRAWING ROUTINE (BASED ON FIRST OCTANT */
        /* ALWAYS). THE X AND Y'S MAY BE FLIPPED IF Y IS THE INDEPENDENT VARIABLE. */
        if (dX >= dY) {     /* if X is the independent variable */
                                int dPr = dY<<1;            /* amount to increment decision if right is chosen (always) */
                                int dPru = dPr - (dX<<1);   /* amount to increment decision if up is chosen */
                                int P = dPr - dX;           /* decision variable start value */
                
                                /* process each point in the line one at a time (just use dX) */
                                for (; dX>=0; dX--) {
                                        newpts->pts[*j].x = Ax;
                                        newpts->pts[*j].y = Ay;
                                        (*j)++;
                
                                        if (P > 0) {    /* is the pixel going right AND up? */
                                                                Ax+=Xincr;      /* increment independent variable */
                                                                Ay+=Yincr;      /* increment dependent variable */
                                                                P+=dPru;        /* increment decision (for up) */
                                        } else {                /* is the pixel just going right? */
                                                                Ax+=Xincr;      /* increment independent variable */
                                                                P+=dPr;         /* increment decision (for right) */
                                        }
                                }
        } else {                    /* if Y is the independent variable */
                                int dPr = dX<<1;            /* amount to increment decision if right is chosen (always) */
                                int dPru = dPr - (dY<<1);   /* amount to increment decision if up is chosen */
                                int P  = dPr - dY;          /* decision variable start value */
                
                                /* process each point in the line one at a time (just use dY) */
                                for (; dY>=0; dY--) {
                                        newpts->pts[*j].x = Ax;
                                        newpts->pts[*j].y = Ay;
                                        (*j)++;
                
                                        if (P > 0) {    /* is the pixel going up AND right? */
                                                                Ax+=Xincr;      /* increment dependent variable */
                                                                Ay+=Yincr;      /* increment independent variable */
                                                                P+=dPru;        /* increment decision (for up) */
                                        } else {                /* is the pixel just going up? */
                                                                Ay+=Yincr;      /* increment independent variable */
                                                                P+=dPr;         /* increment decision (for right) */
                                        }
                                }
        }
}

static void lialg_compute_chain_code(point_list *pts) {
        int i;

        for (i = 0; i < (pts->npts - 1); i++) {
                                pen_point *startpt = &(pts->pts[i]);
                                pen_point *endpt = &(pts->pts[i+1]);
                                int dx = endpt->x - startpt->x;
                                int dy = endpt->y - startpt->y;
                                int tmp = quadr(likeatan(dy, dx));
                                int dircode = (12 - tmp) % 8;
                
                                startpt->chaincode = dircode;
        }
}


static void lialg_compute_unit_chain_code(point_list *pts) {
        int i;

        for (i = 0; i < (pts->npts - 1); i++) {
                                pen_point *startpt = &(pts->pts[i]);
                                pen_point *endpt = &(pts->pts[i+1]);
                                int dx = endpt->x - startpt->x;
                                int dy = endpt->y - startpt->y;
                                int dircode = lialg_dctbl[dx+1][dy+1];

                                startpt->chaincode = dircode;
        }
}


static region_list *lialg_compute_regions(point_list *pts) {
                region_list *regions;
                region_list *curr_reg;
                int *R[2 + LP_FILTER_ITERS];
                int *junk;
                int *curr, *next;
                int i, j;

                /* Initialize low-pass filter parameters if necessary. */
                if (lialg_lpfconst == -1)
                                lialg_compute_lpf_parameters();

        /* Allocate a 2 x pts->npts array for use in computing the (filtered) Angle set, A_n. */
        junk = malloc((2 + LP_FILTER_ITERS) * pts->npts*sizeof(int));
        for (i = 0; i < (2 + LP_FILTER_ITERS); i++)
                R[i] = junk + (i * pts->npts);
        curr = R[0];

        /* Compute the Angle set, A, in the first element of array R. */
        /* Values in R are in degrees, x 100. */
        curr[0] = 18000;                                /* a_0 */
        for (i = 1; i < (pts->npts - 1); i++) {
                                int d_i = pts->pts[i].chaincode;
                                int d_iminusone = pts->pts[i-1].chaincode;
                                int a_i;
                
                                if (d_iminusone < d_i)
                                        d_iminusone += 8;
                
                                a_i = (d_iminusone - d_i) % 8;
                
                                /* convert to degrees, x 100 */
                                curr[i] = ((12 - a_i) % 8) * 45 * 100;
        }
        curr[pts->npts - 1]     = 18000;                /* a_L-1 */

        /* Perform a number of filtering iterations. */
        next = R[1];
        for (j = 0; j < LP_FILTER_ITERS; j++, curr = R[j], next = R[j+1]) {
                                for (i = 0; i < pts->npts; i++) {
                                        int k;
                
                                        next[i] = 0;
                
                                        for (k = i - LP_FILTER_WIDTH; k <= i + LP_FILTER_WIDTH; k++) {
                                                int oldval = (k < 0 || k >= pts->npts) ? 18000 : curr[k];
                                                next[i] += oldval * lialg_lpfwts[k - (i - LP_FILTER_WIDTH)];    /* overflow? */
                                        }
                
                                        next[i] /= lialg_lpfconst;
                                }
        }

        /* Do final thresholding around PI. */
        /* curr and next are set-up correctly at end of previous loop! */
        for (i = 0; i < pts->npts; i++)
                                next[i] = (abs(curr[i] - 18000) < LP_FILTER_THLD) ? 18000 : curr[i];
        curr = next;

        /* Debugging. */
        if (lidebug > 1) {
                                for (i = 0; i < pts->npts; i++) {
                                        fprint(2, "%3d:  (%P)  %lud  ",
                                                                i, pts->pts[i].Point, pts->pts[i].chaincode);
                                        for (j = 0; j < 2 + LP_FILTER_ITERS; j++)
                                                fprint(2, "%d  ", R[j][i]);
                                        fprint(2, "\n");
                                }
        }

        /* Do the region segmentation. */
        {
                                int start, end;
                                int currtype;
                
#define RGN_TYPE(val) (((val)==18000)?RGN_PLAIN:((val)<18000?RGN_CONCAVE:RGN_CONVEX))
                
                                start = 0;
                                currtype = RGN_TYPE(curr[0]);
                                regions = malloc(sizeof(region_list));
                                curr_reg = regions;
                                curr_reg->start = start;
                                curr_reg->end = 0;
                                curr_reg->type = currtype;
                                curr_reg->next = nil;
                                for (i = 1; i < pts->npts; i++) {
                                        int nexttype = RGN_TYPE(curr[i]);
                
                                        if (nexttype != currtype) {
                                                                region_list *next_reg;
                                
                                                                end = i - 1;
                                                                curr_reg->end = end;
                                                                if (lidebug > 1)
                                                                        fprint(2, "  (%d, %d), %d\n", start, end, currtype);
                                
                                                                start = i;
                                                                currtype = nexttype;
                                                                next_reg = malloc(sizeof(region_list));
                                                                next_reg->start = start;
                                                                next_reg->end = 0;
                                                                next_reg->type = nexttype;
                                                                next_reg->next = nil;
                                
                                                                curr_reg->next = next_reg;
                                                                curr_reg = next_reg;
                                        }
                                }
                                end = i - 1;
                                curr_reg->end = end;
                                if (lidebug > 1)
                                        fprint(2, "  (%d, %d), %d\n", start, end, currtype);

                                /* Filter out convex/concave regions that are too short. */
                                for (curr_reg = regions; curr_reg; curr_reg = curr_reg->next)
                                        if (curr_reg->type == RGN_PLAIN) {
                                                                region_list *next_reg;
                                
                                                                for (next_reg = curr_reg->next;
                                                                         next_reg != nil &&
                                                                           (next_reg->end - next_reg->start) < LP_FILTER_MIN;
                                                                         next_reg = curr_reg->next) {
                                                                        /* next_reg must not be plain, and it must be followed by a plain */
                                                                        /* assert(next_reg->type != RGN_PLAIN); */
                                                                        /* assert(next_reg->next != nil && (next_reg->next)->type == RGN_PLAIN); */
                                
                                                                        curr_reg->next = (next_reg->next)->next;
                                                                        curr_reg->end = (next_reg->next)->end;
                                
                                                                        free(next_reg->next);
                                                                        free(next_reg);
                                                                }
                                        }
                
                                /* Add-in pseudo-extremes. */
                                {
                                        region_list *tmp, *prev_reg;
                
                                        tmp = regions;
                                        regions = nil;
                                        prev_reg = nil;
                                        for (curr_reg = tmp; curr_reg; curr_reg = curr_reg->next) {
                                                if (curr_reg->type == RGN_PLAIN) {
                                                        int arclen = lialg_compute_pathlen_subset(pts,
                                                                                                                                                curr_reg->start,
                                                                                                                                                curr_reg->end);
                                                        int dx = pts->pts[curr_reg->end].x -
                                                          pts->pts[curr_reg->start].x;
                                                        int dy = pts->pts[curr_reg->end].y -
                                                          pts->pts[curr_reg->start].y;
                                                        int chordlen = isqrt(10000 * (dx * dx + dy * dy));
                                                        int atcr = (chordlen == 0) ? 0 : (100 * arclen + chordlen / 2) / chordlen;

                                                        if (lidebug)
                                                                fprint(2, "%d, %d, %d\n", arclen, chordlen, atcr);
                
                                                        /* Split region if necessary. */
                                                        if (arclen >= PE_AL_THLD && atcr >= PE_ATCR_THLD) {
                                                                int mid = curr_reg->start + (curr_reg->end - curr_reg->start) / 2;
                                                                int end = curr_reg->end;
                                                                region_list *saved_next = curr_reg->next;
                
                                                                curr_reg->end = mid - 1;
                                                                if (prev_reg == nil)
                                                                        regions = curr_reg;
                                                                else
                                                                        prev_reg->next = curr_reg;
                                                                prev_reg = curr_reg;
                
                                                                /* curr_reg = (region_list *)safe_malloc(sizeof(region_list));*/
                                                                curr_reg = malloc(sizeof(region_list));
                                                                curr_reg->start = mid;
                                                                curr_reg->end = mid;
                                                                curr_reg->type = RGN_PSEUDO;
                                                                curr_reg->next = nil;
                                                                prev_reg->next = curr_reg;
                                                                prev_reg = curr_reg;
                
                                                                /* curr_reg = (region_list *)malloc(sizeof(region_list)); */
                                                                curr_reg = malloc(sizeof(region_list));
                                                                curr_reg->start = mid + 1;
                                                                curr_reg->end = end;
                                                                curr_reg->type = RGN_PLAIN;
                                                                curr_reg->next = nil;
                                                                prev_reg->next = curr_reg;
                                                                prev_reg = curr_reg;
                
                                                                curr_reg->next = saved_next;
                                                                continue;
                                                        }
                                                }
                
                                                if (prev_reg == nil)
                                                        regions = curr_reg;
                                                else
                                                        prev_reg->next = curr_reg;
                                                prev_reg = curr_reg;
                                        }
                                }
                }

        free(junk);
        return(regions);
}


static point_list *lialg_compute_dompts(point_list *pts, region_list *regions) {
        point_list *dpts;
        int ndpts;
        int *cas;
        int nonplain;
        region_list *r;
                region_list *curr;
                int dp;
                int previx;
                int currix;

        /* Compute contour angle set. */
        cas = lialg_compute_contour_angle_set(pts, regions);

        /* Dominant points include:  start_pt, end_pt, extrema_of_non_plain_regions, midpts of the preceding. */
        nonplain = 0;
        for (r = regions; r != nil; r = r->next)
                                if (r->type != RGN_PLAIN)
                                                nonplain++;
        ndpts = 2 * (2 + nonplain) - 1;
        /* dpts = (point_list *)safe_malloc(sizeof(point_list)); */
        dpts = malloc(sizeof(point_list));
        dpts->pts = mallocz(ndpts*sizeof(pen_point), 1);
        if (dpts->pts == nil) {
                                free(dpts);
                                return(nil);
        }
        dpts->npts = ndpts;
        dpts->next = nil;

        /* Pick out dominant points. */

                /* Record start point. */
                dp = 0;
                previx = 0;
                dpts->pts[dp++] = pts->pts[previx];

                for (curr = regions; curr != nil; curr = curr->next)
                        if (curr->type != RGN_PLAIN) {
                                                int max_v = 0;
                                                int min_v = 0x7fffffff; /* maxint */
                                                int max_ix = -1;
                                                int min_ix = -1;
                                                int i;
                
                                                for (i = curr->start; i <= curr->end; i++) {
                                                        int v = cas[i];
                                                        if (v > max_v) { max_v = v; max_ix = i; }
                                                        if (v < min_v) { min_v = v; min_ix = i; }
                                                        if (lidebug > 1)
                                                                fprint(2, "  %d\n", v);
                                                }
                
                                                currix = (curr->type == RGN_CONVEX ? max_ix : min_ix);
                
                                                /* Record midpoint. */
                                                dpts->pts[dp++] = pts->pts[previx + (currix - previx) / 2];
                
                                                /* Record extreme point. */
                                                dpts->pts[dp++] = pts->pts[currix];
                
                                                previx = currix;
                        }

                /* Record last mid-point and end point. */
                currix = pts->npts - 1;
                dpts->pts[dp++] = pts->pts[previx + (currix - previx) / 2];
                dpts->pts[dp] = pts->pts[currix];

        /* Compute chain-code. */
        lialg_compute_chain_code(dpts);

        free(cas);
        return(dpts);
}


static int *lialg_compute_contour_angle_set(point_list *pts,
                                                                                           region_list *regions) {
                int *V;
                region_list *curr_reg;
                int i;

                V = malloc(pts->npts*sizeof(int));

                V[0] = 18000;
                for (curr_reg = regions; curr_reg != nil; curr_reg = curr_reg->next) {
                                for (i = curr_reg->start; i <= curr_reg->end; i++) {
                                        if (curr_reg->type == RGN_PLAIN) {
                                                                V[i] = 18000;
                                        } else {
                                                                /* For now, simply choose the mid-point. */
                                                                int isMidPt = i == (curr_reg->start +
                                                                                                         (curr_reg->end - curr_reg->start) / 2);
                                                                V[i] = (curr_reg->type == RGN_CONVEX)
                                                                  ? (isMidPt ? 18000 : 0)
                                                                  : (isMidPt ? 0 : 18000);
                                        }
                                }
        }
        V[pts->npts - 1] = 18000;

        return(V);
}


/*
 *  First compute the similarity between the two strings.
 *  If it's above a threshold, compute the distance between
 *  the two and return it as the ``score.''
 *  Otherwise, return the constant WORST_SCORE.
 *
 */
static void lialg_score_stroke(point_list *input_dompts, point_list *curr_dompts, int *sim, int *dist) {
        *sim = MIN_SIM;
        *dist = MAX_DIST;

        *sim = lialg_compute_similarity(input_dompts, curr_dompts);
        if (*sim < SIM_THLD) goto done;

        *dist = lialg_compute_distance(input_dompts, curr_dompts);

done:
        if (lidebug)
                fprint(2, "%d, %d\n", *sim, *dist);
}


static int lialg_compute_similarity(point_list *input_dompts, point_list *curr_dompts) {
        int sim;
        point_list *A, *B;
        int N, M;
        int **G;
        int *junk;
        int i, j;

        /* A is the     longer sequence, length N. */
        /* B is the shorter sequence, length M. */
                if (input_dompts->npts >= curr_dompts->npts) {
                                A = input_dompts;
                                N = input_dompts->npts;
                                B = curr_dompts;
                                M = curr_dompts->npts;
        } else {
                                A = curr_dompts;
                                N = curr_dompts->npts;
                                B = input_dompts;
                                M = input_dompts->npts;
                }

                /* Allocate and initialize the Gain matrix, G. */
                /* The size of G is M x (N + 1). */
                /* Note that row 0 is unused. */
                /* Similarities are x 10. */
                G = malloc(M*sizeof(int *));
                junk = malloc(M * (N + 1) * sizeof(int));
                for (i = 0; i < M; i++)
                        G[i] = junk + (i * (N + 1));

                for (i = 1; i < M; i++) {
                        int bval = B->pts[i-1].chaincode;

                        /* Source column. */
                        G[i][0] = 0;

                        for (j = 1; j < N; j++) {
                                int aval = A->pts[j-1].chaincode;
                                int diff = abs(bval - aval);
                                if (diff > 4) diff = 8 - diff;

                                G[i][j] = (diff == 0)
                                  ? 10
                                  : (diff == 1)
                                  ? 6
                                  : 0;
                        }

                        /* Sink column. */
                        G[i][N] = 0;
                }

        /* Do the DP algorithm. */
        /* Proceed in column order, from highest column to the lowest. */
        /* Within each column, proceed from the highest row to the lowest. */
        /* Skip the highest column. */
                for (j = N - 1; j >= 0; j--)
                        for (i = M - 1; i > 0; i--) {
                                int max = G[i][j + 1];

                                if (i < (M - 1)) {
                                        int tmp = G[i + 1][j + 1];
                                        if (tmp > max) max = tmp;
                                }

                                G[i][j] += max;
                }

                sim = (10 * G[1][0] + (N - 1) / 2) / (N - 1);

                if (G != nil)
                                free(G);
                if (junk != nil)
                                free(junk);
                return(sim);
}


static int lialg_compute_distance(point_list *input_dompts,
                                                                   point_list *curr_dompts) {
                int dist;
                point_list *A, *B;
                int N, M;
                int **C;
                int *junk;
                int *BE;
                int *TE;
                int i, j;

                /* A is the     longer sequence, length N. */
                /* B is the shorter sequence, length M. */
                if (input_dompts->npts >= curr_dompts->npts) {
                                A = input_dompts;
                                N = input_dompts->npts;
                                B = curr_dompts;
                                M = curr_dompts->npts;
                }
                else {
                                A = curr_dompts;
                                N = curr_dompts->npts;
                                B = input_dompts;
                                M = input_dompts->npts;
                }

                /* Construct the helper vectors, BE and TE, which say for each column */
                /* what are the ``bottom'' and ``top'' rows of interest. */
                BE = malloc((N + 1)*sizeof(int));
                TE = malloc((N + 1)*sizeof(int));

                for (j = 1; j <= N; j++) {
                                int bot, top;

                                bot = j + (M - DP_BAND);
                                if (bot > M) bot = M;
                                BE[j] = bot;

                                top = j - (N - DP_BAND);
                                if (top < 1) top = 1;
                                TE[j] = top;
                }

                /* Allocate and initialize the Cost matrix, C. */
                /* The size of C is (M + 1) x (N + 1). */
                /* Note that row and column 0 are unused. */
                /* Costs are x 100. */
                /*      C = (int **)safe_malloc((M + 1) * sizeof(int *)); */
                C = malloc((M + 1)*sizeof( int *));
                junk = malloc((M + 1) * (N + 1)*sizeof(int));
                for (i = 0; i <= M; i++)
                                C[i] = junk + (i * (N + 1));

                for (i = 1; i <= M; i++) {
                                int bx = B->pts[i-1].x;
                                int by = B->pts[i-1].y;

                                for (j = 1; j <= N; j++) {
                                                int ax = A->pts[j-1].x;
                                                int ay = A->pts[j-1].y;
                                                int dx = bx - ax;
                                                int dy = by - ay;
                                                int dist = isqrt(10000 * (dx * dx + dy * dy));
                                
                                                C[i][j] = dist;
                                }
                }

                /* Do the DP algorithm. */
                /* Proceed in column order, from highest column to the lowest. */
                /* Within each column, proceed from the highest row to the lowest. */
                for (j = N; j > 0; j--)
                                for (i = M; i > 0; i--) {
                                                int min = MAX_DIST;
                
                                                if (i > BE[j] || i < TE[j] || (j == N && i == M))
                                                                continue;
                
                                                if (j < N) {
                                                                if (i >= TE[j+1]) {
                                                                                int tmp = C[i][j+1];
                                                                                if (tmp < min)
                                                                                                min = tmp;
                                                                }
                
                                                                if (i < M) {
                                                                                int tmp = C[i+1][j+1];
                                                                                if (tmp < min)
                                                                                                min = tmp;
                                                                }
                                                }
                
                                                if (i < BE[j]) {
                                                                int tmp = C[i+1][j];
                                                                if (tmp < min) min = tmp;
                                                }
                
                                                C[i][j] += min;
                                }

                dist = (C[1][1] + N / 2) / N;

                if (C != nil) free(C);
                if (junk != nil) free(junk);
                if (BE != nil) free(BE);
                if (TE != nil) free(TE);
                return(dist);
}


/*************************************************************

  Digest-processing routines

 *************************************************************/

static int lialg_read_classifier_digest(rClassifier *rec) {
        int nclasses;
        FILE *fp;

        /* Try to open the corresponding digest file. */
        {
                                char *clx_path;
                                char *dot;
                
                                /* Get a copy of the filename, with some room on the end. */
                                /*      clx_path = safe_malloc(strlen(rec->file_name) + 5); */
                                clx_path = malloc((strlen(rec->file_name) + 5) *sizeof(char));
                                strcpy(clx_path, rec->file_name);
                
                                /* Truncate the path after the last dot. */
                                dot = strrchr(clx_path, '.');
                                if (dot == nil) { free(clx_path); return(-1); }
                                *(dot + 1) = 0;
                
                                /* Append the classifier-digest extension. */
                                strcat(clx_path, "clx");
                
                                fp = fopen(clx_path, "r");
                                if (fp == nil) {
                                                free(clx_path);
                                                return(-1);
                                }
                
                                free(clx_path);
        }

        /* Read-in the name and dominant points for each class. */
        for (nclasses = 0; !feof(fp); nclasses++) {
                point_list *dpts = nil;
                char class[BUFSIZ];
                int npts;
                int j;

                if (fscanf(fp, "%s %d", class, &npts) != 2) {
                        if (feof(fp)) break;

                        goto failed;
                }
                rec->cnames[nclasses] = strdup(class);

                /* Allocate a dominant-points list. */
                /* dpts = (point_list *)safe_malloc(sizeof(point_list)); */
                dpts = malloc(sizeof(point_list));
                dpts->pts = mallocz(npts*sizeof(pen_point), 1);
                if (dpts->pts == nil) goto failed;
                dpts->npts = npts;
                dpts->next = nil;

                /* Read in each point. */
                for (j = 0; j < npts; j++) {
                        int x, y;

                        if (fscanf(fp, "%d %d", &x, &y) != 2) goto failed;
                        dpts->pts[j].x = x;
                        dpts->pts[j].y = y;
                }

                /* Compute the chain-code. */
                lialg_compute_chain_code(dpts);

                /* Store the list in the rec data structure. */
                rec->dompts[nclasses] = dpts;

                continue;

failed:
                fprint(2, "read_classifier_digest failed...\n");
                for (; nclasses >= 0; nclasses--) {
                        if (rec->cnames[nclasses] != nil) {
                                free(rec->cnames[nclasses]);
                                rec->cnames[nclasses] = nil;
                        }
                        if (rec->dompts[nclasses] != nil) {
                                delete_examples(rec->dompts[nclasses]);
                                rec->dompts[nclasses] = nil;
                        }
                }
                if (dpts != nil)
                        delete_examples(dpts);
                fclose(fp);
                return(-1);
        }

        fclose(fp);
        return(0);
}


/*************************************************************

  Canonicalization routines

 *************************************************************/

static int lialg_canonicalize_examples(rClassifier *rec) {
        int i;
        int nclasses;

        if (lidebug) {
                fprint(2, "lialg_canonicalize_examples working on %s\n",
                                rec->file_name);
        }
        /* Initialize canonical-example arrays. */
        for (i = 0; i < MAXSCLASSES; i++) {
                rec->canonex[i] = nil;
        }

        /* Figure out number of classes. */
        for (nclasses = 0;
                  nclasses < MAXSCLASSES && rec->cnames[nclasses] != nil;
                  nclasses++)
                ;

        /* Canonicalize the examples for each class. */
        for (i = 0; i < nclasses; i++) {
                int j, k;
                int nex;
                point_list *pts, *tmp, *avg;
                int maxxrange, maxyrange;
                int minx, miny, maxx, maxy;
                int avgxrange, avgyrange, avgxoff, avgyoff, avgscale;

                
                if (lidebug) {
                        fprint(2, "lialg_canonicalize_examples working on class %s\n",
                                        rec->cnames[i]);
                }
                /* Make a copy of the examples. */
                pts = nil;
                tmp = rec->ex[i];
                for (nex = 0; tmp != nil; nex++, tmp = tmp->next) {
                        if ((pts = add_example(pts, tmp->npts, tmp->pts)) == nil) {
                                delete_examples(pts);
                                return(-1);
                        }
                }

                /* Canonicalize each example, and derive the max x and y ranges. */
                maxxrange = 0;
                maxyrange = 0;
                for (j = 0, tmp = pts; j < nex; j++, tmp = tmp->next) {
                        if (lialg_canonicalize_example_stroke(tmp) != 0) {
                if (lidebug) {
                                        fprint(2, "lialg_canonicalize_example_stroke returned error\n");
                                }
                                return(-1);
                        }

                        if (tmp->xrange > maxxrange) maxxrange = tmp->xrange;
                        if (tmp->yrange > maxyrange) maxyrange = tmp->yrange;
                }

                /* Normalize max ranges. */
                if (((100 * maxxrange + CANONICAL_X / 2) / CANONICAL_X) >
                        ((100 * maxyrange + CANONICAL_Y / 2) / CANONICAL_Y)) {
                        maxyrange = (maxyrange * CANONICAL_X + maxxrange / 2) / maxxrange;
                        maxxrange = CANONICAL_X;
                }
                else {
                        maxxrange = (maxxrange * CANONICAL_Y + maxyrange / 2) / maxyrange;
                        maxyrange = CANONICAL_Y;
                }

                /* Re-scale each example to max ranges. */
                for (j = 0, tmp = pts; j < nex; j++, tmp = tmp->next) {
                        int scalex = (tmp->xrange == 0) ? 100 : (100 * maxxrange + tmp->xrange / 2) / tmp->xrange;
                        int scaley = (tmp->yrange == 0) ? 100 : (100 * maxyrange + tmp->yrange / 2) / tmp->yrange;
                        if (lialg_translate_points(tmp, 0, 0, scalex, scaley) != 0) {
                                delete_examples(pts);
                                return(-1);
                        }
                }

                /* Average the examples; leave average in first example. */
                avg = pts;                              /* careful aliasing!! */
                for (k = 0; k < NCANONICAL; k++) {
                        int xsum = 0;
                        int ysum = 0;

                        for (j = 0, tmp = pts; j < nex; j++, tmp = tmp->next) {
                                                xsum += tmp->pts[k].x;
                                                ysum += tmp->pts[k].y;
                        }

                        avg->pts[k].x = (xsum + j / 2) / j;
                        avg->pts[k].y = (ysum + j / 2) / j;
                }

                /* Compute BB of averaged stroke and re-scale. */
                lialg_get_bounding_box(avg, &minx, &miny, &maxx, &maxy);
                avgxrange = maxx - minx;
                avgyrange = maxy - miny;
                avgscale = (((100 * avgxrange + CANONICAL_X / 2) / CANONICAL_X) >
                                        ((100 * avgyrange + CANONICAL_Y / 2) / CANONICAL_Y))
                  ? (100 * CANONICAL_X + avgxrange / 2) / avgxrange
                  : (100 * CANONICAL_Y + avgyrange / 2) / avgyrange;
                if (lialg_translate_points(avg, minx, miny, avgscale, avgscale) != 0) {
                        delete_examples(pts);
                        return(-1);
                }

                /* Re-compute the x and y ranges and center the stroke. */
                lialg_get_bounding_box(avg, &minx, &miny, &maxx, &maxy);
                avgxrange = maxx - minx;
                avgyrange = maxy - miny;
                avgxoff = -((CANONICAL_X - avgxrange + 1) / 2);
                avgyoff = -((CANONICAL_Y - avgyrange + 1) / 2);
                if (lialg_translate_points(avg, avgxoff, avgyoff, 100, 100) != 0) {
                        delete_examples(pts);
                        return(-1);
                }

                /* Create a point list to serve as the ``canonical representation. */
                if ((rec->canonex[i] = add_example(nil, avg->npts, avg->pts)) == nil) {
                        delete_examples(pts);
                        return(-1);
                }
                (rec->canonex[i])->xrange = maxx - minx;
                (rec->canonex[i])->yrange = maxy - miny;

                if (lidebug) {
                        fprint(2, "%s, avgpts = %d\n", rec->cnames[i], avg->npts);
                        for (j = 0; j < avg->npts; j++) {
                                                fprint(2, "  (%P)\n", avg->pts[j].Point);
                        }
                }

                /* Compute dominant points of canonical representation. */
                rec->dompts[i] = lialg_compute_dominant_points(avg);

                /* Clean up. */
                delete_examples(pts);
        }

        /* Sanity check. */
        for (i = 0; i < nclasses; i++) {
                char *best_name = lialg_recognize_stroke(rec, rec->canonex[i]);

                if (best_name != rec->cnames[i])
                        fprint(2, "%s, best = %s\n", rec->cnames[i], best_name);
        }

        return(0);
}


static int lialg_canonicalize_example_stroke(point_list *points) {
        int minx, miny, maxx, maxy, xrange, yrange, scale;

        /* Filter out points that are too close. */
        if (lialg_filter_points(points) != 0) return(-1);

        /* Must be at least two points! */
        if (points->npts < 2) {
                if (lidebug) {
                        fprint(2, "lialg_canonicalize_example_stroke: npts=%d\n",
                                        points->npts);
                }
                return(-1);
        }

        /* Scale up to avoid conversion errors. */
        lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
        xrange = maxx - minx;
        yrange = maxy - miny;
        scale = (((100 * xrange + CANONICAL_X / 2) / CANONICAL_X) >
                         ((100 * yrange + CANONICAL_Y / 2) / CANONICAL_Y))
          ? (100 * CANONICAL_X + xrange / 2) / xrange
          : (100 * CANONICAL_Y + yrange / 2) / yrange;
        if (lialg_translate_points(points, minx, miny, scale, scale) != 0) {
                if (lidebug) {
                        fprint(2, "lialg_translate_points (minx=%d,miny=%d,scale=%d) returned error\n", minx, miny, scale);
                }
                return(-1);
        }

        /* Compute an equivalent stroke with equi-distant points. */
        if (lialg_compute_equipoints(points) != 0) return(-1);

        /* Re-translate the points to the origin. */
        lialg_get_bounding_box(points, &minx, &miny, &maxx, &maxy);
        if (lialg_translate_points(points, minx, miny, 100, 100) != 0) {
                if (lidebug) {
                        fprint(2, "lialg_translate_points (minx=%d,miny=%d) returned error\n", minx, miny);
                }
                return(-1);
        }

        /* Store the x and y ranges in the point list. */
        xrange = maxx - minx;
        yrange = maxy - miny;
        points->xrange = xrange;
        points->yrange = yrange;

        if (lidebug) {
                int i;
                fprint(2, "Canonicalized:   %d, %d, %d, %d\n", minx, miny, maxx, maxy);
                for (i = 0; i < points->npts; i++)
                        fprint(2, "      (%P)\n", points->pts[i].Point);
                fflush(stderr);
        }

        return(0);
}


static int lialg_compute_equipoints(point_list *points) {
        pen_point *equipoints = mallocz(NCANONICAL*sizeof(pen_point), 1);
        int nequipoints = 0;
        int pathlen = lialg_compute_pathlen(points);
        int equidist = (pathlen + (NCANONICAL - 1) / 2) / (NCANONICAL - 1);
        int i;
        int dist_since_last_eqpt;
        int remaining_seglen;
        int dist_to_next_eqpt;

        if (equipoints == nil) {
                fprint(2, "can't allocate memory in lialg_compute_equipoints");
                return(-1);
        }

        if (lidebug) {
                fprint(2, "compute_equipoints:  npts = %d, pathlen = %d, equidist = %d\n",
                                points->npts, pathlen, equidist);
                fflush(stderr);
        }

        /* First original point is an equipoint. */
        equipoints[0] = points->pts[0];
        nequipoints++;
        dist_since_last_eqpt = 0;

        for (i = 1; i < points->npts; i++) {
                int dx1 = points->pts[i].x - points->pts[i-1].x;
                int dy1 = points->pts[i].y - points->pts[i-1].y;
                int endx = 100 * points->pts[i-1].x;
                int endy = 100 * points->pts[i-1].y;
                remaining_seglen = isqrt(10000 * (dx1 * dx1 + dy1 * dy1));
                dist_to_next_eqpt = equidist - dist_since_last_eqpt;

                while (remaining_seglen >= dist_to_next_eqpt) {
                        if (dx1 == 0) {
                                /* x-coordinate stays the same */
                                if (dy1 >= 0)
                                        endy += dist_to_next_eqpt;
                                else
                                        endy -= dist_to_next_eqpt;
                        }
                        else {
                                int slope = (100 * dy1 + dx1 / 2) / dx1;
                                int tmp = isqrt(10000 + slope * slope);
                                int dx = (100 * dist_to_next_eqpt + tmp / 2) / tmp;
                                int dy = (slope * dx + 50) / 100;

                                if (dy < 0) dy = -dy;
                                if (dx1 >= 0)
                                        endx += dx;
                                else
                                        endx -= dx;
                                if (dy1 >= 0)
                                        endy += dy;
                                else
                                        endy -= dy;
                        }

                        equipoints[nequipoints].x = (endx + 50) / 100;
                        equipoints[nequipoints].y = (endy + 50) / 100;
                        nequipoints++;
/*          assert(nequipoints <= NCANONICAL);*/
                        dist_since_last_eqpt = 0;
                        remaining_seglen -= dist_to_next_eqpt;
                        dist_to_next_eqpt = equidist;
                }

                dist_since_last_eqpt += remaining_seglen;
        }

        /* Take care of last equipoint. */
        if (nequipoints == NCANONICAL) {
                /* Good. */
        } else if (nequipoints == (NCANONICAL - 1)) {
                /* Make last original point the last equipoint. */
                equipoints[nequipoints] = points->pts[points->npts - 1];
        } else {
          if (lidebug) {
                fprint(2,"lialg_compute_equipoints: nequipoints = %d\n", 
                                nequipoints);
          }
/*      assert(false);*/
                return(-1);
        }

        points->npts = NCANONICAL;
        delete_pen_point_array(points->pts);
        points->pts = equipoints;
        return(0);
}


/*************************************************************

  Utility routines

 *************************************************************/

/* Result is x 100. */
static int lialg_compute_pathlen(point_list *points) {
        return(lialg_compute_pathlen_subset(points, 0, points->npts - 1));
}


/* Result is x 100. */
static int lialg_compute_pathlen_subset(point_list *points,
                                                                                   int start, int end) {
        int pathlen;
        int i;

        pathlen = 0;
        for (i = start + 1; i <= end; i++) {
                int dx = points->pts[i].x - points->pts[i-1].x;
                int dy = points->pts[i].y - points->pts[i-1].y;
                int dist = isqrt(10000 * (dx * dx + dy * dy));
                pathlen += dist;
        }

        return(pathlen);
}


/* Note that this does NOT update points->xrange and points->yrange! */
static int lialg_filter_points(point_list *points) {
        int filtered_npts;
        pen_point *filtered_pts = mallocz(points->npts*sizeof(pen_point), 1);
        int i;

        if (filtered_pts == nil) {
                fprint(2, "can't allocate memory in lialg_filter_points");
                return(-1);
        }

        filtered_pts[0] = points->pts[0];
        filtered_npts = 1;
        for (i = 1; i < points->npts; i++) {
                int j = filtered_npts - 1;
                int dx = points->pts[i].x - filtered_pts[j].x;
                int dy = points->pts[i].y - filtered_pts[j].y;
                int magsq = dx * dx + dy * dy;

                if (magsq >= DIST_SQ_THRESHOLD) {
                        filtered_pts[filtered_npts] = points->pts[i];
                        filtered_npts++;
                }
        }

        points->npts = filtered_npts;
        delete_pen_point_array(points->pts);
        points->pts = filtered_pts;
        return(0);
}


/* scalex and scaley are x 100. */
/* Note that this does NOT update points->xrange and points->yrange! */
static int lialg_translate_points(point_list *points,
                                                                   int minx, int miny,
                                                                   int scalex, int scaley) {
        int i;

        for (i = 0; i < points->npts; i++) {
                                points->pts[i].x = ((points->pts[i].x - minx) * scalex + 50) / 100;
                                points->pts[i].y = ((points->pts[i].y - miny) * scaley + 50) / 100;
        }

        return(0);
}


static void lialg_get_bounding_box(point_list *points,
                                                                        int *pminx, int *pminy,
                                                                        int *pmaxx, int *pmaxy) {
        int minx, miny, maxx, maxy;
        int i;

        minx = maxx = points->pts[0].x;
        miny = maxy = points->pts[0].y;
        for (i = 1; i < points->npts; i++) {
                                pen_point *pt = &(points->pts[i]);
                                if (pt->x < minx) minx = pt->x;
                                else if (pt->x > maxx) maxx = pt->x;
                                if (pt->y < miny) miny = pt->y;
                                else if (pt->y > maxy) maxy = pt->y;
        }

        *pminx = minx;
        *pminy = miny;
        *pmaxx = maxx;
        *pmaxy = maxy;
}


int wtvals[] = {100, 104, 117, 143, 189, 271, 422};

static void lialg_compute_lpf_parameters(void) {
        int i;

                for (i = LP_FILTER_WIDTH; i >= 0; i--) {
//              double x = 0.04 * (i * i);
//              double tmp = 100.0 * exp(x);
//              int wt = floor((double)tmp);
                                int wt = wtvals[i];
                                lialg_lpfwts[LP_FILTER_WIDTH - i] = wt;
                                lialg_lpfwts[LP_FILTER_WIDTH + i] = wt;
                }
                lialg_lpfconst = 0;
                for (i = 0; i < (2 * LP_FILTER_WIDTH + 1); i++) {
                                lialg_lpfconst += lialg_lpfwts[i];
                }
}


/* Code from Joseph Hall (jnhall@sat.mot.com). */
static int isqrt(int n) {
                register int i;
                register long k0, k1, nn;

                for (nn = i = n, k0 = 2; i > 0; i >>= 2, k0 <<= 1)
                                ;
                nn <<= 2;
        for (;;) {
                                k1 = (nn / k0 + k0) >> 1;
                                if (((k0 ^ k1) & ~1) == 0)
                                break;
                                k0 = k1;
                }
                return (int) ((k1 + 1) >> 1);
}


/* Helper routines from Mark Hayter. */
static int likeatan(int tantop, int tanbot) { 
                int t;
                /* Use tan(theta)=top/bot --> order for t */
                /* t in range 0..0x40000 */

                if ((tantop == 0) && (tanbot == 0)) 
                                t = 0;
        else
                {
                                t = (tantop << 16) / (abs(tantop) + abs(tanbot));
                                if (tanbot < 0) 
                                                t = 0x20000 - t;
                                else 
                                                if (tantop < 0) t = 0x40000 + t;
                }
                return t;
}


static int quadr(int t) {
        return (8 - (((t + 0x4000) >> 15) & 7)) & 7;
}