harvey/sys/src/libscribble/li_recognizer.c

/*
 * This file is part of the UCB release of Plan 9. It is subject to the license
 * terms in the LICENSE file found in the top-level directory of this
 * distribution and at http://akaros.cs.berkeley.edu/files/Plan9License. No
 * part of the UCB release of Plan 9, including this file, may be copied,
 * modified, propagated, or distributed except according to the terms contained
 * in the LICENSE file.
 */

/*
 *  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* i)
{
	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 r)
{
  /*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 r, rc* rc, uint u, Stroke* s, rec_element* re, bool b)
{
  /*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 r, char* c, char* d)
{ 
		/*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 r, char* c, char* d)
{
		/*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 r, char* c, char* d, wordset w)
{
		/*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 r, wordset w)
{
  /*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 r, letterset* l, wordset w)
{
		/*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 r, letterset* l, wordset w)
{
		/*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 b)
{
		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 r, rc* 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 r)
{
		/*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 r, uint* u, Stroke** s)
{
		/*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 r, uint u, Stroke* s)
{
		/*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 b, 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 r)
{
		/*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 r, char* c, xgesture x, void* v)
{
		/*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), %lu\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 (%lu)  ", 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)  %lu  ",
								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 int32_t 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;
}