li_recognizer.c

这是一个同样来自贝尔实验室的和UNIX有着渊源的操作系统, 其简洁的设计和实现易于我们学习和理解

								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);