r-h.c

Hausdorff Distance for Image Recognition

    thFree(lowth);

    model->trans = trans;

    return;

  bailout:

    if (lowth != (TransHash *)NULL) {
	thFree(lowth);
	}

    if (highth != (TransHash *)NULL) {
	thFree(highth);
	}

    if (trans != NULLLIST) {
	free_translist(trans);
	}

    model->trans = NULLLIST;
    }
    
/*
 * Check the translations that we were handed (in model->trans); remove
 * those where the thresholds aren't satisfied.
 */
static void
checkTransMtoI(image_info *image, model_info *model)
{
    ListNode ln;
    long pointval = 0;
    long *vals = (long *)NULL;
    transval *t;
    int x, y;
    int i;
    int model_required;

    int nover;
    int novermax;

    /* Stuff cached from image and model */
    point *model_pts;
    unsigned nmodel_pts;
    LongImage image_dtrans;
    long model_thresh;
    double model_frac;
    
    assert(image != (image_info *)NULL);
    assert(model != (model_info *)NULL);
    assert(image->dtrans != (LongImage)NULL);
    assert(model->trans != NULLLIST);

    model_pts = model->pts;
    nmodel_pts = model->npts;
    model_thresh = model->model_thresh;
    model_frac = model->model_frac;
    image_dtrans = image->dtrans;

    assert(model_pts != (point *)NULL);

    model_required = (int)(ceil(nmodel_pts * model_frac));
    novermax = nmodel_pts - model_required;

    if (model_frac != 1.) {
	vals = (long *)malloc(nmodel_pts * sizeof(long));
	if (vals == (long *)NULL) {
	    goto bailout;
	    }
	}

    for (ln = liFirst(model->trans); ln != NULLLISTNODE; ) {
	t = (transval *)liGet(ln);
	assert(t != (transval *)NULL);
	x = t->transpos.x;
	y = t->transpos.y;

	if (model_frac == 1.) {
	    pointval = 0;
	    for (i = 0; i < nmodel_pts; i++) {
		pointval = MAX(pointval, imRef(image_dtrans,
					       model_pts[i].x + x,
					       model_pts[i].y + y));
		if (pointval > model_thresh) {
		    break;
		    }
		}
	    }
	else {
	    nover = 0;
	    for (i = 0; i < nmodel_pts; i++) {
		/* Probe at this model point's translated location */
		pointval = vals[i] = imRef(image_dtrans,
					   model_pts[i].x + x,
					   model_pts[i].y + y);
		if (pointval > model_thresh) {
		    nover++;
		    if (nover > novermax) {
			break;
			}
		    }

		}
	    
	    if (i == nmodel_pts) {
		pointval = pickNth(vals, (int)nmodel_pts, model_required - 1);
		assert(pointval <= model_thresh);
		}
	    }

	if (pointval <= model_thresh) {
	    t->forward_val = pointval;
	    if (model_frac == 1.) {
		t->forward_frac = 1.;
		}
	    else {
		t->forward_frac = frac_lower(vals, (int)nmodel_pts,
					     model_thresh);
		assert(t->forward_frac >= model_frac - EPSILON);
		}
		
	    ln = liNext(ln);
	    }
	else {
	    free_trans(t);
	    ln = liRemAndNext(ln);
	    }
	}

    if (model_frac != 1) {
	assert(vals != (long *)NULL);
	free((void *)vals);
	}

    return;

  bailout:
    if (vals) {
	free((void *)vals);
	}

    free_translist(model->trans);
    model->trans = NULLLIST;
    return;
    }

/*
 * This routine scans the points listed in modeltrans, looking for points
 * where at least image_frac of the image points under the model lie within
 * image_thresh of model points.
 */
static void
findTransItoM(image_info *image, model_info *model, int revstyle)
{
    ListNode ln;
    transval *t;
    int x, y;
    int x2, y2;
    long pointval;
    int nvals;
    int nimage = 0;
    int whichval;
    long *vals;
    LongImage image_dtrans;
    LongImage model_dtrans;

    assert(image != (image_info *)NULL);
    assert(model != (model_info *)NULL);
    assert(model->trans != NULLLIST);
    assert(model->dtrans != (LongImage)NULL);
    assert(image->dtrans != (LongImage)NULL);

    /* Cache stuff */
    image_dtrans = image->dtrans;
    model_dtrans = model->dtrans;

    /* Figure out how many image pixels might be visible */
    vals = (long *)malloc(sizeof(*vals) *
			  MIN(image->npts,
			      (model->xsize + 2 * MODEL_BORD) *
			      (model->ysize + 2 * MODEL_BORD)));
    if (vals == (long *)NULL) {
	free_translist(model->trans);
	model->trans = NULLLIST;
	return;
	}

    for (ln = liFirst(model->trans); ln != NULLLISTNODE; ) {
	t = (transval *)liGet(ln);
	assert(t != (transval *)NULL);
	x = t->transpos.x;
	y = t->transpos.y;

	if (revstyle == REVERSE_BOX) {
	    
	    /*
	     * For this, it's almost always wrong to iterate over the
	     * point list, since most of them will be outside the domain
	     * of the model (not underneath it at the moment). We therefore
	     * scan the actual image looking for set points. Actually,
	     * we cheat: we scan the distance transform of the image,
	     * looking for 0s.
	     */
	    nvals = 0;
	    for (y2 = y; y2 < y + (int)model->ysize; y2++) {
		for (x2 = x; x2 < x + (int)model->xsize; x2++) {
		    if (imRef(image_dtrans, x2, y2) == 0) {
			vals[nvals++] = imRef(model_dtrans, x2 - x, y2 - y);
			}
		    }
		}
	    
	    whichval = (int)(ceil(model->image_frac * nvals));
	    pointval = pickNth(vals, nvals, whichval - 1);
	    }
	else {
	    assert(revstyle == REVERSE_ALLIMAGE);
	    /* We need to scan over the entire image. */
	    /* Oh, does my brain hurt... */

	    /*
	     * OK - the deal is that we've got the model's distance transform,
	     * which is the size of the model plus MODEL_BORD pixels border.
	     * If image_thresh is <= MODEL_BORD * DSCALE then we can use this
	     * distance transform: for each set pixel in the image, if it
	     * falls within the distance transform, probe and use the value.
	     * If it doesn't, then its value must be > image_thresh, and so
	     * we can use a large value instead of bothering to probe.
	     *
	     * If image_thresh is > MODEL_BORD, then you're out of luck -
	     * it's going to be hard to figure out some of these values.
	     */
	    if (model->image_thresh < DSCALE * MODEL_BORD) {
		nvals = 0;
		for (y2 = MAX(0, y - MODEL_BORD);
		     y2 < MIN(y + (int)model->ysize + MODEL_BORD,
			      image->ysize);
		     y2++) {
		    for (x2 = MAX(0, x - MODEL_BORD);
			 x2 < MIN(x + (int)model->xsize + MODEL_BORD,
				  image->xsize);
			 x2++) {
			if (imRef(image_dtrans, x2, y2) == 0) {
			    vals[nvals++] =
				imRef(model_dtrans, x2 - x, y2 - y);
			    }
			}
		    }
		nimage = image->npts;
		}
	    else {
		/*
		 * Do it the hard way.
		 * Scan the entire image. If something lies inside the
		 * model's dtrans, probe and use. If it lies far enough
		 * outside it that we know its value is large, ignore.
		 * If it lies somewhere in between, then we get to do
		 * a full find-nearest-neighbour with all the model points...
		 * yay.
		 */
		int allowed = (int)ceil(model->image_thresh / (double)DSCALE);
		int i;

		panic("untested code");	/* I want to know if someone uses this */

		nvals = 0;
		for (y2 = MAX(0, y - allowed);
		     y2 < MIN(y + (int)model->ysize + allowed,
			      image->ysize);
		     y2++) {
		    for (x2 = MAX(0, x - allowed);
			 x2 < MIN(x + (int)model->xsize + allowed,
				  image->xsize);
			 x2++) {
			if (imRef(image_dtrans, x2, y2) == 0) {
			    if ((x2 >= x - MODEL_BORD) &&
				(x2 < x + (int)model->xsize + MODEL_BORD) &&
				(y2 >= y - MODEL_BORD) &&
				(y2 < y + (int)model->ysize + MODEL_BORD)) {
				/* It falls inside - yay */
				vals[nvals++] = imRef(model_dtrans,
						      x2 - x, y2 - y);
				}
			    else {
				/* Go search. */
				/*
				 * N.B. This assumes that we're using L2
				 * as our underlying norm.
				 */
				assert(model->npts > 0);
				assert(model->pts != (point *)NULL);
				
				vals[nvals] =
				    (int)(DSCALE *
					  sqrt((double)
					       ((x2 - x - model->pts[0].x) *
						(x2 - x - model->pts[0].x) +
						(y2 - y - model->pts[0].y) *
						(y2 - y - model->pts[0].y))));
				for (i = 1; i < model->npts; i++) {
				    vals[nvals] =
					MIN(vals[nvals],
					    (int)(DSCALE *
						  sqrt((double)
						       ((x2 - x -
							 model->pts[0].x) *
							(x2 - x -
							 model->pts[0].x) +
							(y2 - y -
							 model->pts[0].y) *
							(y2 - y -
							 model->pts[0].y)))));
				    }
				nvals++;
				}
			    }
			}
		    }
		nimage = image->npts;
		}
	    /*
	     * OK - we have nvals values which have been probed
	     * and nimage - nvals which haven't been, but which would
	     * be large if they were.
	     */
	    whichval = (int)(ceil(model->image_frac * nimage));
	    if (whichval > nvals) {
		/* Reject */
		pointval = model->image_thresh * 2 + 1;
		}
	    else {
		/* Give it a chance */
		pointval = pickNth(vals, nvals, whichval - 1);
		}
	    }

	if (pointval <= model->image_thresh) {
	    /* A good point. Keep it. */
	    t->reverse_val = pointval;
	    if (revstyle == REVERSE_BOX) {
		if (nvals == 0) {
		    t->reverse_frac = 1.;
		    }
		else {
		    t->reverse_frac =
			frac_lower(vals, nvals, model->image_thresh);
		    }
		t->reverse_num = nvals;
		}
	    else {
		if (nvals == 0) {
		    if (nimage == 0) {
			t->reverse_frac = 1.;
			}
		    else {
			t->reverse_frac = 0.;
			}
		    }
		else {
		    t->reverse_frac =
			frac_lower(vals, nvals, model->image_thresh) * nvals /
			    nimage;
		    }
		t->reverse_num = nimage;
		}
	    assert(t->reverse_frac >= model->image_frac - EPSILON);

	    /* and advance */
	    ln = liNext(ln);
	    }
	else {
	    /* Free the hook from modeltrans */
	    free_trans(t);
	    ln = liRemAndNext(ln);
	    }
	}

    free((void *)vals);

    return;
    }

/*
 * Scan a region in transformation space from lowX to highX (inclusive,
 * exclusive) in X, stepping by stepX in X, and similarly for Y.
 * Probe at locations which are cell centres:
 * lowX + stepX / 2, lowX + 3 * stepX / 2, etc.
 *
 * Look for places where there may be something interesting within a cell
 * of size (stepX,stepY) centred around the current position. For
 * each interesting cell you find, add a flag to th. If needdata is TRUE,
 * attach a stransval structure to that flag. Cells are represented
 * by their top left corner position, not the centre position, even though
 * the centre position is what is actually probed.
 *
 * A cell is interesting if the centre position's value is <= thresh.
 */
static boolean