ref: 75c92428225428c8fde2d015f010e608a0b12f1d
dir: purgatorio/appl/lib/strokes/strokes.b
implement Strokes;
#
# this Limbo code is derived from C code that had the following
# copyright notice, which i reproduce as requested
#
# 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.
#
#
# the copyright notice really did end in the middle of the sentence
#
#
# Limbo version for Inferno by forsyth@vitanuova.com, Vita Nuova, September 2001
#
#
# the code is reasonably close to the algorithms described in
# ``On-line Handwritten Alphanumeric Character Recognition Using Dominant Stroke in Strokes'',
# Xiaolin Li and Dit-Yan Yueng, Department of Computer Science,
# Hong Kong University of Science and Technology, Hong Kong (23 August 1996)
#
include "sys.m";
sys: Sys;
include "strokes.m";
MAXINT: con 16r7FFFFFFF;
# Dynamic programming parameters
DP_BAND: con 3;
SIM_THLD: con 60; # x100
#DIST_THLD: con 3200; # x100
DIST_THLD: con 3300; # x100
# Low-pass filter parameters -- empirically derived
LP_FILTER_WIDTH: con 6;
LP_FILTER_ITERS: con 8;
LP_FILTER_THLD: con 250; # x100
LP_FILTER_MIN: con 5;
# Pseudo-extrema parameters -- empirically derived
PE_AL_THLD: con 1500; # x100
PE_ATCR_THLD: con 135; # x100
# Pre-processing and canonicalization parameters
CANONICAL_X: con 108;
CANONICAL_Y: con 128;
DIST_SQ_THRESHOLD: con 3*3;
# direction-code table; indexed by dx, dy
dctbl := array[] of {array[] of {1, 0, 7}, array[] of {2, MAXINT, 6}, array[] of {3, 4, 5}};
# low-pass filter weights
lpfwts := array[2 * LP_FILTER_WIDTH + 1] of int;
lpfconst := -1;
lidebug: con 0;
stderr: ref Sys->FD;
# x := 0.04 * (i * i);
# wtvals[i] = floor(100.0 * exp(x));
wtvals := array[] of {100, 104, 117, 143, 189, 271, 422};
init()
{
sys = load Sys Sys->PATH;
if(lidebug)
stderr = sys->fildes(2);
for(i := LP_FILTER_WIDTH; i >= 0; i--){
wt := wtvals[i];
lpfwts[LP_FILTER_WIDTH - i] = wt;
lpfwts[LP_FILTER_WIDTH + i] = wt;
}
lpfconst = 0;
for(i = 0; i < (2 * LP_FILTER_WIDTH + 1); i++)
lpfconst += lpfwts[i];
}
Stroke.new(n: int): ref Stroke
{
return ref Stroke(n, array[n] of Penpoint, 0, 0);
}
Stroke.trim(ps: self ref Stroke, n: int)
{
ps.npts = n;
ps.pts = ps.pts[0:n];
}
Stroke.copy(ps: self ref Stroke): ref Stroke
{
n := ps.npts;
a := array[n] of Penpoint;
a[0:] = ps.pts[0:n];
return ref Stroke(n, a, ps.xrange, ps.yrange);
}
#
# return the bounding box of a set of points
# (note: unlike Draw->Rectangle, the region is closed)
#
Stroke.bbox(ps: self ref Stroke): (int, int, int, int)
{
minx := maxx := ps.pts[0].x;
miny := maxy := ps.pts[0].y;
for(i := 1; i < ps.npts; i++){
pt := ps.pts[i];
if(pt.x < minx)
minx = pt.x;
if(pt.x > maxx)
maxx = pt.x;
if(pt.y < miny)
miny = pt.y;
if(pt.y > maxy)
maxy = pt.y;
}
return (minx, miny, maxx, maxy); # warning: closed interval
}
Stroke.center(ps: self ref Stroke)
{
(minx, miny, maxx, maxy) := ps.bbox();
ps.xrange = maxx-minx;
ps.yrange = maxy-miny;
avgxoff := -((CANONICAL_X - ps.xrange + 1) / 2);
avgyoff := -((CANONICAL_Y - ps.yrange + 1) / 2);
ps.translate(avgxoff, avgyoff, 100, 100);
}
Stroke.scaleup(ps: self ref Stroke): int
{
(minx, miny, maxx, maxy) := ps.bbox();
xrange := maxx - minx;
yrange := maxy - miny;
scale: int;
if(((100 * xrange + CANONICAL_X / 2) / CANONICAL_X) >
((100 * yrange + CANONICAL_Y / 2) / CANONICAL_Y))
scale = (100 * CANONICAL_X + xrange / 2) / xrange;
else
scale = (100 * CANONICAL_Y + yrange / 2) / yrange;
ps.translate(minx, miny, scale, scale);
return scale;
}
# scalex and scaley are x 100.
# Note that this does NOT update points.xrange and points.yrange!
Stroke.translate(ps: self ref Stroke, minx: int, miny: int, scalex: int, scaley: int)
{
for(i := 0; i < ps.npts; i++){
ps.pts[i].x = ((ps.pts[i].x - minx) * scalex + 50) / 100;
ps.pts[i].y = ((ps.pts[i].y - miny) * scaley + 50) / 100;
}
}
TAP_PATHLEN: con 10*100; # x100
Classifier.match(rec: self ref Classifier, stroke: ref Stroke): (int, string)
{
if(stroke.npts < 1)
return (-1, nil);
# Check for tap.
# First thing is to filter out ``close points.''
stroke = stroke.filter();
# 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 || stroke.length() < TAP_PATHLEN)
return (-1, "tap");
preprocess_stroke(stroke);
# Compute its dominant points.
dompts := stroke.interpolate().dominant();
best_dist := MAXDIST;
best_i := -1;
best_name: string;
# Score input stroke against every class in classifier.
for(i := 0; i < len rec.cnames; i++){
name := rec.cnames[i];
(sim, dist) := score_stroke(dompts, rec.dompts[i]);
if(dist < MAXDIST)
sys->fprint(stderr, " (%s, %d, %d)", name, sim, dist);
if(dist < DIST_THLD){
if(lidebug)
sys->fprint(stderr, " (%s, %d, %d)", name, sim, dist);
# Is it the best so far?
if(dist < best_dist){
best_dist = dist;
best_i = i;
best_name = name;
}
}
}
if(lidebug)
sys->fprint(stderr, "\n");
return (best_i, best_name);
}
preprocess_stroke(s: ref Stroke)
{
# Filter out points that are too close.
# We did this earlier, when we checked for a tap.
# s = s.filter();
# assert(s.npts > 0);
# Scale up to avoid conversion errors.
s.scaleup();
# Center the stroke.
s.center();
if(lidebug){
(minx, miny, maxx, maxy) := s.bbox();
sys->fprint(stderr, "After pre-processing: [ %d %d %d %d]\n",
minx, miny, maxx, maxy);
printpoints(stderr, s, "\n");
}
}
#
# return the dominant points of Stroke s, assuming s has been through interpolation
#
Stroke.dominant(s: self ref Stroke): ref Stroke
{
regions := s.regions();
# Dominant points are: start, end, extrema of non plain regions, and midpoints of the preceding.
nonplain := 0;
for(r := regions; r != nil; r = r.next)
if(r.rtype != Rplain)
nonplain++;
dom := Stroke.new(1 + 2*nonplain + 2);
# Pick out dominant points.
# start point
dp := 0;
previx := 0;
dom.pts[dp++] = s.pts[previx];
currix: int;
cas := s.contourangles(regions);
for(r = regions; r != nil; r = r.next)
if(r.rtype != Rplain){
max_v := 0;
min_v := MAXINT;
max_ix := -1;
min_ix := -1;
for(i := r.start; i <= r.end; i++){
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)
sys->fprint(stderr, " %d\n", v);
}
if(r.rtype == Rconvex)
currix = max_ix;
else
currix = min_ix;
dom.pts[dp++] = s.pts[(previx+currix)/2]; # midpoint
dom.pts[dp++] = s.pts[currix]; # extreme
previx = currix;
}
# last midpoint, and end point
lastp := s.npts - 1;
dom.pts[dp++] = s.pts[(previx+lastp)/2];
dom.pts[dp++] = s.pts[lastp];
dom.trim(dp);
# Compute chain-code.
compute_chain_code(dom);
return dom;
}
Stroke.contourangles(s: self ref Stroke, regions: ref Region): array of int
{
V := array[s.npts] of int;
V[0] = 18000;
for(r := regions; r != nil; r = r.next){
for(i := r.start; i <= r.end; i++){
if(r.rtype == Rplain){
V[i] = 18000;
}else{
# For now, simply choose the mid-point.
ismidpt := i == (r.start + r.end)/2;
if(ismidpt ^ (r.rtype!=Rconvex))
V[i] = 18000;
else
V[i] = 0;
}
}
}
V[s.npts - 1] = 18000;
return V;
}
Stroke.interpolate(s: self ref Stroke): ref Stroke
{
# Compute an upper-bound on the number of interpolated points
maxpts := s.npts;
for(i := 0; i < s.npts - 1; i++){
a := s.pts[i];
b := s.pts[i+1];
maxpts += abs(a.x - b.x) + abs(a.y - b.y);
}
# Allocate an array of the maximum size
newpts := Stroke.new(maxpts);
# Interpolate each of the segments.
j := 0;
for(i = 0; i < s.npts - 1; i++){
j = bresline(s.pts[i], s.pts[i+1], newpts, j);
j--; # end point gets recorded as start point of next segment!
}
# Add-in last point and trim
newpts.pts[j++] = s.pts[s.npts - 1];
newpts.trim(j);
if(lidebug){
sys->fprint(stderr, "After interpolation:\n");
printpoints(stderr, newpts, "\n");
}
# Compute the chain code for P (the list of points).
compute_unit_chain_code(newpts);
return newpts;
}
# This implementation is due to an anonymous page on the net
bresline(startpt: Penpoint, endpt: Penpoint, newpts: ref Stroke, j: int): int
{
x0 := startpt.x;
x1 := endpt.x;
y0 := startpt.y;
y1 := endpt.y;
stepx := 1;
dx := x1-x0;
if(dx < 0){
dx = -dx;
stepx = -1;
}
dx <<= 1;
stepy := 1;
dy := y1-y0;
if(dy < 0){
dy = -dy;
stepy = -1;
}
dy <<= 1;
newpts.pts[j++] = (x0, y0, 0);
if(dx >= dy){
e := dy - (dx>>1);
while(x0 != x1){
if(e >= 0){
y0 += stepy;
e -= dx;
}
x0 += stepx;
e += dy;
newpts.pts[j++] = (x0, y0, 0);
}
}else{
e := dx - (dy>>1);
while(y0 != y1){
if(e >= 0){
x0 += stepx;
e -= dy;
}
y0 += stepy;
e += dx;
newpts.pts[j++] = (x0, y0, 0);
}
}
return j;
}
compute_chain_code(pts: ref Stroke)
{
for(i := 0; i < pts.npts - 1; i++){
dx := pts.pts[i+1].x - pts.pts[i].x;
dy := pts.pts[i+1].y - pts.pts[i].y;
pts.pts[i].chaincode = (12 - quadr(likeatan(dy, dx))) % 8;
}
}
compute_unit_chain_code(pts: ref Stroke)
{
for(i := 0; i < pts.npts - 1; i++){
dx := pts.pts[i+1].x - pts.pts[i].x;
dy := pts.pts[i+1].y - pts.pts[i].y;
pts.pts[i].chaincode = dctbl[dx+1][dy+1];
}
}
Stroke.regions(pts: self ref Stroke): ref Region
{
# Allocate a 2 x pts.npts array for use in computing the (filtered) Angle set, A_n.
R := array[] of {0 to LP_FILTER_ITERS+1 => array[pts.npts] of int};
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++){
d_i := pts.pts[i].chaincode;
d_iminusone := pts.pts[i-1].chaincode;
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; ){
for(i = 0; i < pts.npts; i++){
next[i] = 0;
for(k := i - LP_FILTER_WIDTH; k <= i + LP_FILTER_WIDTH; k++){
oldval: int;
if(k < 0 || k >= pts.npts)
oldval = 18000;
else
oldval = curr[k];
next[i] += oldval * lpfwts[k - (i - LP_FILTER_WIDTH)]; # overflow?
}
next[i] /= lpfconst;
}
j++;
curr = R[j];
next = R[j+1];
}
# Do final thresholding around PI.
# curr and next are set-up correctly at end of previous loop!
for(i = 0; i < pts.npts; i++)
if(abs(curr[i] - 18000) < LP_FILTER_THLD)
next[i] = 18000;
else
next[i] = curr[i];
curr = next;
# Debugging.
if(lidebug > 1){
for(i = 0; i < pts.npts; i++){
p := pts.pts[i];
sys->fprint(stderr, "%3d: (%d %d) %ud ",
i, p.x, p.y, p.chaincode);
for(j = 0; j < 2 + LP_FILTER_ITERS; j++)
sys->fprint(stderr, "%d ", R[j][i]);
sys->fprint(stderr, "\n");
}
}
# Do the region segmentation.
r := regions := ref Region(regiontype(curr[0]), 0, 0, nil);
for(i = 1; i < pts.npts; i++){
t := regiontype(curr[i]);
if(t != r.rtype){
r.end = i-1;
if(lidebug > 1)
sys->fprint(stderr, " (%d, %d) %d\n", r.start, r.end, r.rtype);
r.next = ref Region(t, i, 0, nil);
r = r.next;
}
}
r.end = i-1;
if(lidebug > 1)
sys->fprint(stderr, " (%d, %d), %d\n", r.start, r.end, r.rtype);
# Filter out convex/concave regions that are too short.
for(r = regions; r != nil; r = r.next)
if(r.rtype == Rplain){
while((nr := r.next) != nil && (nr.end - nr.start) < LP_FILTER_MIN){
# nr must not be plain, and it must be followed by a plain
# assert(nr.rtype != Rplain);
# assert(nr.next != nil && (nr.next).rtype == Rplain);
if(nr.next == nil){
sys->fprint(stderr, "recog: nr.next==nil\n"); # can't happen
break;
}
r.next = nr.next.next;
r.end = nr.next.end;
}
}
# Add-in pseudo-extremes.
for(r = regions; r != nil; r = r.next)
if(r.rtype == Rplain){
arclen := pts.pathlen(r.start, r.end);
dx := pts.pts[r.end].x - pts.pts[r.start].x;
dy := pts.pts[r.end].y - pts.pts[r.start].y;
chordlen := sqrt(100*100 * (dx*dx + dy*dy));
atcr := 0;
if(chordlen)
atcr = (100*arclen + chordlen/2) / chordlen;
if(lidebug)
sys->fprint(stderr, "%d, %d, %d\n", arclen, chordlen, atcr);
# Split region if necessary.
if(arclen >= PE_AL_THLD && atcr >= PE_ATCR_THLD){
mid := (r.start + r.end)/2;
end := r.end;
r.end = mid - 1;
r = r.next = ref Region(Rpseudo, mid, mid,
ref Region(Rplain, mid+1, end, r.next));
}
}
return regions;
}
regiontype(val: int): int
{
if(val == 18000)
return Rplain;
if(val < 18000)
return Rconcave;
return Rconvex;
}
#
# return the similarity of two strokes and,
# if similar, the distance between them;
# if dissimilar, the distance is MAXDIST)
#
score_stroke(a: ref Stroke, b: ref Stroke): (int, int)
{
sim := compute_similarity(a, b);
if(sim < SIM_THLD)
return (sim, MAXDIST);
return (sim, compute_distance(a, b));
}
compute_similarity(A: ref Stroke, B: ref Stroke): int
{
# A is the longer sequence, length N.
# B is the shorter sequence, length M.
if(A.npts < B.npts){
t := A; A = B; B = t;
}
N := A.npts;
M := B.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 := array[M] of array of int;
for(i := 1; i < M; i++){
G[i] = array[N+1] of int;
bcode := B.pts[i-1].chaincode;
G[i][0] = 0; # source column
for(j := 1; j < N; j++){
diff := abs(bcode - A.pts[j-1].chaincode);
if(diff > 4)
diff = 8 - diff; # symmetry
v := 0;
if(diff == 0)
v = 10;
else if(diff == 1)
v = 6;
G[i][j] = v;
}
G[i][N] = 0; # sink column
}
# 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--){
max := G[i][j + 1];
if(i < M-1){
t := G[i + 1][j + 1];
if(t > max)
max = t;
}
G[i][j] += max;
}
return (10*G[1][0] + (N-1)/2) / (N-1);
}
compute_distance(A: ref Stroke, B: ref Stroke): int
{
# A is the longer sequence, length N.
# B is the shorter sequence, length M.
if(A.npts < B.npts){
t := A; A = B; B = t;
}
N := A.npts;
M := B.npts;
# Construct the helper vectors, BE and TE, which say for each column
# what are the ``bottom'' and ``top'' rows of interest.
BE := array[N+1] of int;
TE := array[N+1] of int;
for(j := 1; j <= N; j++){
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 := array[M+1] of array of int;
for(i := 1; i <= M; i++){
C[i] = array[N+1] of int;
bx := B.pts[i-1].x;
by := B.pts[i-1].y;
for(j = 1; j <= N; j++){
ax := A.pts[j-1].x;
ay := A.pts[j-1].y;
dx := bx - ax;
dy := by - ay;
dist := sqrt(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--){
min := MAXDIST;
if(i > BE[j] || i < TE[j] || (j == N && i == M))
continue;
if(j < N){
if(i >= TE[j+1]){
tmp := C[i][j+1];
if(tmp < min)
min = tmp;
}
if(i < M){
tmp := C[i+1][j+1];
if(tmp < min)
min = tmp;
}
}
if(i < BE[j]){
tmp := C[i+1][j];
if(tmp < min)
min = tmp;
}
C[i][j] += min;
}
return (C[1][1] + N / 2) / N; # dist
}
# Result is x 100.
Stroke.length(s: self ref Stroke): int
{
return s.pathlen(0, s.npts-1);
}
# Result is x 100.
Stroke.pathlen(s: self ref Stroke, first: int, last: int): int
{
l := 0;
for(i := first + 1; i <= last; i++){
dx := s.pts[i].x - s.pts[i-1].x;
dy := s.pts[i].y - s.pts[i-1].y;
l += sqrt(100*100 * (dx*dx + dy*dy));
}
return l;
}
# Note that this does NOT update points.xrange and points.yrange!
Stroke.filter(s: self ref Stroke): ref Stroke
{
pts := array[s.npts] of Penpoint;
pts[0] = s.pts[0];
npts := 1;
for(i := 1; i < s.npts; i++){
j := npts - 1;
dx := s.pts[i].x - pts[j].x;
dy := s.pts[i].y - pts[j].y;
magsq := dx * dx + dy * dy;
if(magsq >= DIST_SQ_THRESHOLD){
pts[npts] = s.pts[i];
npts++;
}
}
return ref Stroke(npts, pts[0:npts], 0, 0);
}
abs(a: int): int
{
if(a < 0)
return -a;
return a;
}
# Code from Joseph Hall (jnhall@sat.mot.com).
sqrt(n: int): int
{
nn := n;
k0 := 2;
for(i := n; i > 0; i >>= 2)
k0 <<= 1;
nn <<= 2;
k1: int;
for(;;){
k1 = (nn / k0 + k0) >> 1;
if(((k0 ^ k1) & ~1) == 0)
break;
k0 = k1;
}
return (k1 + 1) >> 1;
}
# Helper routines from Mark Hayter.
likeatan(tantop: int, tanbot: int): int
{
# Use tan(theta)=top/bot -. order for t
# t in range 0..16r40000
if(tantop == 0 && tantop == 0)
return 0;
t := (tantop << 16) / (abs(tantop) + abs(tanbot));
if(tanbot < 0)
t = 16r20000 - t;
else if(tantop < 0)
t = 16r40000 + t;
return t;
}
quadr(t: int): int
{
return (8 - (((t + 16r4000) >> 15) & 7)) & 7;
}
printpoints(fd: ref Sys->FD, pts: ref Stroke, sep: string)
{
for(j := 0; j < pts.npts; j++){
p := pts.pts[j];
sys->fprint(fd, " (%d %d) %ud%s", p.x, p.y, pts.pts[j].chaincode, sep);
}
}