r-h.c

Hausdorff Distance for Image Recognition

/*
 *
 * $Header: /usr/u/wjr/vision/r-h/RCS/r-h.c,v 2.9 1994/09/01 21:14:09 wjr Exp $
 *
 * Copyright (c) 1990, 1991, 1992, 1993 Cornell University.  All Rights
 * Reserved.  
 *  
 * Copyright (c) 1991, 1992 Xerox Corporation.  All Rights Reserved.  
 *  
 * Use, reproduction, preparation of derivative works, and distribution
 * of this software is permitted.  Any copy of this software or of any
 * derivative work must include both the above copyright notices of
 * Cornell University and Xerox Corporation and this paragraph.  Any
 * distribution of this software or derivative works must comply with all
 * applicable United States export control laws.  This software is made
 * available AS IS, and XEROX CORPORATION DISCLAIMS ALL WARRANTIES,
 * EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE,
 * AND NOTWITHSTANDING ANY OTHER PROVISION CONTAINED HEREIN, ANY
 * LIABILITY FOR DAMAGES RESULTING FROM THE SOFTWARE OR ITS USE IS
 * EXPRESSLY DISCLAIMED, WHETHER ARISING IN CONTRACT, TORT (INCLUDING
 * NEGLIGENCE) OR STRICT LIABILITY, EVEN IF XEROX CORPORATION IS ADVISED
 * OF THE POSSIBILITY OF SUCH DAMAGES.
 */

static char rcsid[] = "@(#)$Header: /usr/u/wjr/vision/r-h/RCS/r-h.c,v 2.9 1994/09/01 21:14:09 wjr Exp $";

#include <stdio.h>
#include "misc.h"
#include <math.h>
#include "image.h"
#include "r-h.h"
#include "list.h"
#include "chunk.h"
#include "transhash.h"
#include "dbug.h"

#include "prof.h"

static void findTransMtoI(image_info *image, model_info *model);
static void checkTransMtoI(image_info *image, model_info *model);
static void findTransItoM(image_info *image, model_info *model,
			  int revstyle);
static boolean findTransMtoISampled(image_info *image, model_info *model,
				    int lowX, int highX, int stepX,
				    int lowY, int highY, int stepY,
				    TransHash *th,
				    long thresh, boolean needdata);
static boolean probeRegions(image_info *image, model_info *model,
			    TransHash *lowth, TransHash *highth,
			    int lowX, int highX, int baseX, int stepX,
			    int oldBaseX, int oldStepX,
			    int lowY, int highY, int baseY, int stepY,
			    int oldBaseY, int oldStepY,
			    long thresh, boolean needdata);

static boolean ensure_dtrans(image_info *image, long thresh);

/* About how fine to slice each dimension at the lowest resolution */
#define	SLICE		4

/* What the resolution refinement between steps is */
#define	REFINE		2

/* The usual FP epsilon */
#define	EPSILON	(1e-5)

/* How much of the model to scan before deciding on early decision */
#define	EARLY_FRAC	0.2

/* How big a border to put on the model's distance transform (in pixels) */
#define	MODEL_BORD	4

void
findTransAllModels(image_info *image,
		   model_info *models, unsigned nmodels,
		   int revstyle)
{
    int i;
    int max_leftborder;
    int max_topborder;
    int max_rightborder;
    int max_bottomborder;

    assert(models != (model_info *)NULL);
    assert(image != (image_info *)NULL);
    assert((revstyle == REVERSE_BOX) || (revstyle == REVERSE_ALLIMAGE) ||
	   (revstyle == REVERSE_NONE));

    if (nmodels == 0) {
	return;		/* What's to do? */
	}

    max_leftborder = models[0].leftborder;
    max_topborder = models[0].topborder;
    max_rightborder = models[0].rightborder;
    max_bottomborder = models[0].bottomborder;
    for (i = 1; i < nmodels; i++) {
	max_leftborder = MAX(max_leftborder, models[i].leftborder);
	max_topborder = MAX(max_topborder, models[i].topborder);
	max_rightborder = MAX(max_rightborder, models[i].rightborder);
	max_bottomborder = MAX(max_bottomborder, models[i].bottomborder);
	}

    if ((image->dtrans == (LongImage)NULL) ||
	(image->leftborder < max_leftborder) ||
	(image->topborder < max_topborder) ||
	(image->rightborder < max_rightborder) ||
	(image->bottomborder < max_bottomborder)) {

	/* Were they being very silly? */
	if ((max_leftborder + max_rightborder + (int)image->xsize <= 0) ||
	    (max_topborder + max_bottomborder + (int)image->ysize <= 0)) {
	    /*
	     * They asked for borders which leave no part of the image.
	     * I refuse to deal with this case properly (i.e. deal with the
	     * model->trans fields appropriately). FIX someday?
	     */
	    return;
	    }
	/* The image's distance transform is not present, or is too small */
	if (image->dtrans != (LongImage)NULL) {
	    imFree(image->dtrans);
	    image->dtrans = (LongImage)NULL;
	    }

	if (image->plusx_dtrans != (ShortImage)NULL) {
	    imFree(image->plusx_dtrans);
	    image->plusx_dtrans = (ShortImage)NULL;
	    }

	/*
	 * Generate the image's distance transform. Allow space for the
	 * borders. Also, we need to make sure that every feasible translation
	 * can be reached by some cell whose centre may be outside the actual
	 * range of feasible translations: in order to cover everything, we
	 * may end up probing at a non-feasible translation.
	 */
	image->leftborder = max_leftborder;
	image->topborder = max_topborder;
	image->rightborder = max_rightborder;
	image->bottomborder = max_bottomborder;

	image->dtrans =
	    dtrans_pts((unsigned)ceil((image->leftborder + image->rightborder + image->xsize) *
				      (1. + .5 / SLICE)),
		       (unsigned)ceil((image->topborder + image->bottomborder + image->ysize) *
				      (1. + .5 / SLICE)),
		       -(int)image->leftborder, -(int)image->topborder,
		       image->npts, image->pts);
	if (image->dtrans == (LongImage)NULL) {
	    return;
	    }
	}

    for (i = 0; i < nmodels; i++) {
	/* We'll need the model's dtrans */
	if (models[i].dtrans == (LongImage)NULL) {
	    models[i].dtrans = dtrans_pts(2 * MODEL_BORD + models[i].xsize,
					  2 * MODEL_BORD + models[i].ysize,
					  -MODEL_BORD, -MODEL_BORD,
					  models[i].npts, models[i].pts);
	    if (models[i].dtrans == (LongImage)NULL) {
		return;
		}
	    }

	if (models[i].trans != NULLLIST) {
	    /* We've been handed a list of translations to look at */
	    checkTransMtoI(image, &models[i]);
	    }
	else {
	    findTransMtoI(image, &models[i]);
	    }

	if (models[i].trans == NULLLIST) {
	    return;
	    }

	if (revstyle != REVERSE_NONE) {
	    findTransItoM(image, &models[i], revstyle);
	    }

	if (models[i].trans == NULLLIST) {
	    return;
	    }
	}

    }

/*
 * Find all translations where the model->image distance is <= thresh.
 * Use pruning tricks. The tricks are complicated a little by the fact
 * that we're using a partial match. We therefore have to consider all the
 * model points.
 * Also, we can't simply abort out of a scan when we find a value that is
 * too large, though we can do something similar by counting. The rest of
 * the pruning is still possible, though.
 *
 * The space we search is the space of all translations. We discretise this
 * space to a resolution of one pixel in x and y.
 *
 * We do the whole mess in a multi-resolution manner:
 * the first pass is very coarse, with very high threshold. It tells us
 * what sections look interesting. Further passes refine these sections
 * at increasingly higher resolution. Since there is a definite speed
 * advantage to looking at a large area in transformation space as a
 * single chunk, as opposed to several smaller chunks, we try to merge
 * adjacent interesting sections before going on to the next pass.
 */
static void
findTransMtoI(image_info *image, model_info *model)
{
    TransHash *lowth = (TransHash *)NULL;
    TransHash *highth = (TransHash *)NULL;
    List trans = NULLLIST;
    ListNode newln;
    transval *newt;
    boolean needdata;
    int x, y;
    int lowX, highX, baseX, stepX, oldBaseX, oldStepX;
    int lowY, highY, baseY, stepY, oldBaseY, oldStepY;
    long curthresh;

    /* Stuff cached from image and model */
    int nmodel_pts;
    
    assert(image != (image_info *)NULL);
    assert(model != (model_info *)NULL);
    assert(model->pts != (point *)NULL);
    assert(image->dtrans != (LongImage)NULL);
    assert(model->stepx > 0);
    assert(model->stepy > 0);

    nmodel_pts = model->npts;

    trans = liNew();
    if (trans == NULLLIST) {
	goto bailout;
	}

    /* Get rid of an annoying case */
    if (nmodel_pts == 0) {
	model->trans = trans;
	return;		/* No transformations */
	}

    /* Now figure out what step size we're using at the lowest resolution */

    /*
     * Find out how big the box we're working in is (borders are inclusive
     * on the low side, exclusive on the high side)
     * and calculate lowX, highX, stepX, ...Y for
     * the lowest resolution.
     *
     * What we are determining here:
     * The absolute bounds on X are lowX..highX-1 inclusive. Ditto Y.
     * The lowest resolution cell is stepX by stepY.
     * The first cell is centred at (baseX, baseY).
     */
    /*
     * Possible bug here - we may not hit all the way up to the high end
     * due to truncation in the division by SLICE. Check and fix.
     *
     * No, I think it's right the way it is now, because there's no
     * number of iterations computed; ...Sampled() does however many it
     * needs to for the current stepX. We need to make stepX a multiple
     * of model->stepx though.
     */
    lowX = -image->leftborder;
    highX = (int)(image->xsize + image->rightborder);
    stepX = MAX(model->stepx, (highX - lowX) / SLICE);
    baseX = lowX;
    /* Round stepX down to a multiple of model->stepx */
    stepX = (stepX / model->stepx) * model->stepx;
    
    lowY = -image->topborder;
    highY = (int)(image->ysize + image->bottomborder);
    stepY = MAX(model->stepy, (highY - lowY) / SLICE);
    baseY = lowY;
    stepY = (stepY / model->stepy) * model->stepy;

    /* Allocate the lowest resolution hash table */
    lowth = thNew(lowX, highX, stepX,
		  lowY, highY, stepY);
    if (lowth == (TransHash *)NULL) {
	goto bailout;
	}

    /*
     * Work out the current threshold.
     *
     * We always probe at the centre of the box, so
     * the actual coordinates that we're probing for are:
     * (x - stepX / 2 .. x + stepX - stepX / 2 - 1) x (....Y)
     * so we need to find the distance to the farthest corner. We will be
     * probing at (x, y).
     * If stepX is odd, then stepX / 2 = (stepX - 1) / 2., and so
     * stepX - stepX / 2 - 1 = (stepX - 1) / 2.
     * If stepX is even, then stepX / 2 = stepX / 2., and so
     * stepX - stepX / 2 - 1 = stepX / 2. - 1. In either case, the
     * farther corner is stepX / 2 away.
     * The farthest away we can get in X is therefore stepX / 2, and
     * similarly in Y. Note that all of these are
     * integer divisions.
     * Assuming L_2, this gives us:
     * model_thresh +
     *	ceil(DSCALE * hypot((double)(stepX / 2), (double)(stepY / 2)))
     * However, we have to correct for the fact that stepX and stepY
     * are constrained to be at least model->stepx and model->stepy, and
     * that model->model_thresh *is the threshold to be used at
     * that finest resolution*.
     */
    curthresh = model->model_thresh +
	ceil(DSCALE * hypot((double)((stepX - model->stepx + 1) / 2),
			    (double)((stepY - model->stepy + 1) / 2)));

    VDEBUG("curthresh = %ld\n", curthresh, 5);

    if (!findTransMtoISampled(image, model,
			      lowX, highX, stepX,
			      lowY, highY, stepY,
			      lowth, curthresh, FALSE)) {
	/* Failed. Bah. */
	goto bailout;
	}

    /* We now have the lowest-resolution hash table. Refine... */
    while ((stepX > model->stepx) || (stepY > model->stepy)) {
	/* update stepX etc. keep old ones in oldStepX etc */
	oldStepX = stepX;
	oldBaseX = baseX;
	stepX = MAX(model->stepx, stepX / REFINE);
	baseX = lowX;
	stepX = (stepX / model->stepx) * model->stepx;
	
	oldStepY = stepY;
	oldBaseY = baseY;
	stepY = MAX(model->stepy, stepY / REFINE);
	baseY = lowY;
	stepY = (stepX / model->stepy) * model->stepy;

	/* Update the threshold */
	curthresh = model->model_thresh +
	    ceil(DSCALE * hypot((double)((stepX - model->stepx + 1) / 2),
				(double)((stepY - model->stepy + 1) / 2)));
	
	VDEBUG("curthresh = %ld\n", curthresh, 5);
	
	/*
	 * Determine if we're in highest res yet; store in needdata.
	 */
	if ((stepX == model->stepx) && (stepY == model->stepy)) {
	    needdata = TRUE;
	    }
	else {
	    needdata = FALSE;
	    }

	/* Allocate the new hash table */
	highth = thNew(lowX, highX, stepX,
		       lowY, highY, stepY);
	if (highth == (TransHash *)NULL) {
	    goto bailout;
	    }
	
	/* Pick out regions */
	if (!probeRegions(image, model, lowth, highth,
			  lowX, highX, baseX, stepX, oldBaseX, oldStepX,
			  lowY, highY, baseY, stepY, oldBaseY, oldStepY,
			  curthresh, needdata)) {
	    goto bailout;
	    }

	/* That generated highth for us - make it the new lowth */
	thFree(lowth);
	lowth = highth;
	highth = (TransHash *)NULL;
	}

    /*
     * We now have lowth - the hash table at highest resolution.
     * Read it out into trans.
     */
    while (thGetLowest(lowth, &x, &y)) {
	newt = (transval *)thGet(lowth, x, y);
	thRem(lowth, x, y);
	assert(newt != (transval *)NULL);

	newln = liAdd(trans, newt);
	if (newln == NULLLISTNODE) {
	    /* In this case, we really need to free all hooks... FIX. */
	    goto bailout;
	    }

	}

    /* Free useless stuff */