scale-h.c

Hausdorff Distance for Image Recognition

	     * 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;
	    }

	/* Any cached box dtranses are out of date */
	clear_sboxdtcache(image);

	image->leftborder = max_leftborder;
	image->topborder = max_topborder;
	image->rightborder = max_rightborder;
	image->bottomborder = max_bottomborder;

	image->dtrans =
	    dtrans_pts(image->leftborder + image->rightborder + image->xsize,
		       image->topborder + image->bottomborder + image->ysize,
		       -image->leftborder, -image->topborder,
		       image->npts, image->pts);
	if (image->dtrans == (LongImage)NULL) {
	    return;
	    }
	}

    for (i = 0; i < nmodels; i++) {
	/* We'll need the max values */
	/* and the scaled X and Y coordinates of all the points */
	if (!compute_scaled_pts(&models[i])) {
	    return;
	    }

	/* Do the match */
	findSTrans(image, &models[i], forwstyle, revstyle);

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

	}
#ifdef INSTRUMENT
    DEBUG("Total: %d probes\n", nprobes);
#endif
    }

/*
 * Perform the matching between this model and this image, using the
 * specified matching styles.
 *
 * Should handle trans!=NULL in cases other than FORWARD_ALL, even if only
 * to assertfail. FIX.
 */
static void
findSTrans(simage_info *image, smodel_info *model, int forwstyle, int revstyle)
{
    ListNode ln;
    stransval *t;

#ifdef	MP
    /* We're running in multiprocess mode. Start up the MP stuff */
    if (!mp_startup(image, model, forwstyle)) {
	return;
	}
#endif

    switch(forwstyle) {
      case FORWARD_ALL:
      case FORWARD_ALL | FORWARD_HILLCLIMB: {
	/*
	 * If we're searching for all matches, can do the two (forward and
	 * reverse) passes separately
	 */
	  
	  
	if (model->trans != NULLLIST) {
	    /* We've been handed a list of transformations to look at */
	    checkSTransMtoI(image, model);
	    }
	else {
	    /* Search for matches */
	    assert(forwstyle == FORWARD_ALL);	/* Can't hillclimb */
	    findSTransMtoI(image, model);
	    }
	
	/* Failed? */
	if (model->trans == NULLLIST) {
	    return;
	    }
	
	break;
	}

      case FORWARD_ONE:
      case FORWARD_ONE | FORWARD_HILLCLIMB: {
	/* Find a single match */
	findSTransOne(image, model, revstyle);
	break;
	}

      case FORWARD_BESTDIST: {
	/* Find the best distance */
	findSTransBestDist(image, model, revstyle);
	break;
	}

      case FORWARD_BESTFRAC: {
	/* Find the best distance */
	findSTransBestFrac(image, model, revstyle);
	break;
	}

      default: {
	panic("bad forwstyle");
	break;
	}

      }

#ifdef	MP
    mp_shutdown();
#endif

    if (forwstyle & FORWARD_HILLCLIMB) {
	/*
	 * We got some matches. Improve them.
	 */
	hillclimb(image, model, revstyle);
	}

    /*
     * OK - now we need to fill in the reverse values for the matches.
     * Some of these may already have been done, as matches used to generate
     * the parameter must pass the reverse test, but there's no way to
     * tell which, at the moment, so we do them all.
     */
    if (revstyle != REVERSE_NONE) {
	checkSTransItoM(image, model, revstyle);
	}
    else {
	/* Clear out all the reverse vals, so we don't report garbage */
	foreach (ln, model->trans) {
	    t = (stransval *)liGet(ln);
	    assert(t != (stransval *)NULL);
	    t->reverse_val = 0;
	    t->reverse_frac = t->reverse_frac2 = 0.;
	    t->reverse_num = 0;
	    }
	}
    
    }
    

/*
 * Find all transformations 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, crossed with the
 * space of scales (indpendent scales in x and y) such that
 * maxtransx >= x >= mintransx
 * maxtransy >= y >= mintransy
 * MIN(1, maxscalex) >= scalex >= model->minscalex
 * MIN(1, maxscaley) >= scaley >= model->minscaley
 * scalex / scaley <= maxskew
 * scaley / scalex <= maxskew.
 * and the scaled model box fits inside the image borders.
 *
 * We discretise each space. We discretise the translation space at a
 * resolution of one pixel in x and y, and discretise the scale space in such
 * a way that one "tick" in the x-scale direction causes each model point to
 * move at most one pixel in x, and similarly for the y dimension.
 *
 * 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.
 *
 */
void
findSTransMtoI(simage_info *image, smodel_info *model)
{
    PriQ lowpq = NULLPRIQ;
    PriQ highpq = NULLPRIQ;
    List trans = NULLLIST;
    PriQNode pqn;
    spqHook *hook;
    ListNode newln;
    stransval *newt;
    boolean needdata;

    int lowX, highX, stepX, oldStepX;
    int lowY, highY, stepY, oldStepY;
    int lowScaleX, highScaleX, stepScaleX, oldStepScaleX;
    int lowScaleY, highScaleY, stepScaleY, oldStepScaleY;
    unsigned box_width, box_height;
    
    assert(image != (simage_info *)NULL);
    assert(model != (smodel_info *)NULL);
    assert(model->pts != (point *)NULL);
    assert(image->dtrans != (LongImage)NULL);
    assert(model->stepx > 0);
    assert(model->stepy > 0);

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

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

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

    calcLowCells(image, model,
		 &lowX, &highX, &stepX, &lowY, &highY, &stepY,
		 &lowScaleX, &highScaleX, &stepScaleX,
		 &lowScaleY, &highScaleY, &stepScaleY);

    /* Allocate the lowest resolution priority queue */
    lowpq = pqNew(falseCompare);
    if (lowpq == NULLPRIQ) {
	goto bailout;
	}

    /* Get the top-level box sizes */
    calcBoxes(model, stepX, stepY, stepScaleX, stepScaleY,
	      &box_width, &box_height);

    /* and make the boxdtrans */
    if (!ensure_sdtrans(image, box_width, box_height, model->model_thresh)) {
	goto bailout;
	}

    DEBUG("box_width = %d, box_height = %d\n", box_width, box_height);

    /*
     * Determine if we're in highest res yet; store in needdata.
     */
    if ((stepX == model->stepx) && (stepY == model->stepy) &&
	(stepScaleX == model->stepx) && (stepScaleY == model->stepy)) {
	needdata = TRUE;
	}
    else {
	needdata = FALSE;
	}
    
#ifdef	INSTRUMENT
    level = 0;
    nprobes1 = nabort1 = npick1 = nearly1 = nwrong1 = nfull1 = 0;
    nprobes = nabort = npick = nearly = nwrong = nfull = 0;
#ifdef	INST2
    {
    int i;
    for (i = 0; i < MAXPT; i++) {
	ntt1[i] = nft1[i] = ntf1[i] = nff1[i] = 0;
	ntt[i] = nft[i] = ntf[i] = nff[i] = 0;
	nbelow1[i] = 0;
	nbelow[i] = 0;
	}
    }
#endif
#endif

    /* Get the coarsest-resolution scan */
    if (!findSTransMtoISampled(image, model,
			       lowX, highX, stepX,
			       lowY, highY, stepY,
			       lowScaleX, highScaleX, stepScaleX,
			       lowScaleY, highScaleY, stepScaleY,
			       lowpq, model->model_thresh, needdata)) {
	/* Failed. Bah. */
	goto bailout;
	}

#ifdef	INSTRUMENT

    DEBUG("after level %d:\n", level);
    DEBUG("%d probes, %d aborted, %d picked (%d early, of which %d wrong), %d full (this level)\n",
	  nprobes1, nabort1, npick1, nearly1, nwrong1, nfull1);
    DEBUG("%d probes, %d aborted, %d picked (%d early, of which %d wrong), %d full (total)\n",
	  nprobes, nabort, npick, nearly, nwrong, nfull);
#ifdef	INST2
    {
    int i;
    fprintf(stderr, "\"ntt%d\n", level);
    for (i = 0; i < model->npts; i++) {
	fprintf(stderr, "%d %d\n", i, ntt1[i]);
	}
    fprintf(stderr, "\n\"ntf%d\n", level);
    for (i = 0; i < model->npts; i++) {
	fprintf(stderr, "%d %d\n", i, ntf1[i]);
	}
    fprintf(stderr, "\n\"nft%d\n", level);
    for (i = 0; i < model->npts; i++) {
	fprintf(stderr, "%d %d\n", i, nft1[i]);
	}
    fprintf(stderr, "\n\"nff%d\n", level);
    for (i = 0; i < model->npts; i++) {
	fprintf(stderr, "%d %d\n", i, nff1[i]);
	}
    }
    {
    GrayImage m;
    int i;
    char name[100];
    m = (GrayImage)imNew(IMAGE_GRAY, model->xsize, model->ysize);
    for (i = 0; i < model->npts; i++) {
	imRef(m, model->pts[i].x, model->pts[i].y) =
	    nbelow1[i] * (COLRNG - 1) / nprobes1;
	}
    (void)sprintf(name, "/tmp/mod%d.pgm", level);
    (void)imSave(m, name);
    imFree(m);
    }
#endif
#endif

    /* We now have the lowest-resolution priority queue. Refine... */
    while ((stepX > model->stepx) || (stepY > model->stepy) ||
	   (stepScaleX > model->stepx) || (stepScaleY > model->stepy)) {
	/* update stepX etc. keep old ones in oldStepX etc */
	oldStepX = stepX;
	oldStepY = stepY;
	oldStepScaleX = stepScaleX;
	oldStepScaleY = stepScaleY;

	updateCellSteps(model,
			oldStepX, oldStepY, oldStepScaleX, oldStepScaleY,
			&stepX, &stepY, &stepScaleX, &stepScaleY);
	/* Update the threshold */
	calcBoxes(model, stepX, stepY, stepScaleX, stepScaleY,
		  &box_width, &box_height);

	if (!ensure_sdtrans(image, box_width, box_height, model->model_thresh)) {
	    goto bailout;
	    }
	
	DEBUG("box_width = %d, box_height = %d\n", box_width, box_height);
	
	/*
	 * Determine if we're in highest res yet; store in needdata.
	 */
	if ((stepX == model->stepx) && (stepY == model->stepy) &&
	    (stepScaleX == model->stepx) && (stepScaleY == model->stepy)) {
	    needdata = TRUE;
	    }
	else {
	    needdata = FALSE;
	    }

	/* Allocate the new priority queue */
	highpq = pqNew(falseCompare);
	if (highpq == NULLPRIQ) {
	    goto bailout;
	    }
	
	shuffle_model(model);
#ifdef	INSTRUMENT
	level++;
	nprobes1 = nabort1 = npick1 = nearly1 = nwrong1 = nfull1 = 0;
#ifdef	INST2
	{
	int i;
	for (i = 0; i < MAXPT; i++) {
	    ntt1[i] = nft1[i] = ntf1[i] = nff1[i] = 0;
	    nbelow1[i] = 0;
	    }
	}
#endif
#endif
	/* Pick out regions */
	if (!probeRegions(image, model, lowpq, highpq,
			  lowX, highX, stepX, oldStepX,
			  lowY, highY, stepY, oldStepY,
			  lowScaleX, highScaleX, stepScaleX, oldStepScaleX,
			  lowScaleY, highScaleY, stepScaleY, oldStepScaleY,
			  model->model_thresh, needdata, -1)) {
	    goto bailout;
	    }

#ifdef	INSTRUMENT
    DEBUG("after level %d:\n", level);
    DEBUG("%d probes, %d aborted, %d picked (%d early, of which %d wrong), %d full (this level)\n",
	  nprobes1, nabort1, npick1, nearly1, nwrong1, nfull1);
    DEBUG("%d probes, %d aborted, %d picked (%d early, of which %d wrong), %d full (total)\n",
	  nprobes, nabort, npick, nearly, nwrong, nfull);
#ifdef	INST2
	{
	int i;
	fprintf(stderr, "\"ntt%d\n", level);
	for (i = 0; i < model->npts; i++) {
	    fprintf(stderr, "%d %d\n", i, ntt1[i]);
	    }
	fprintf(stderr, "\n\"ntf%d\n", level);
	for (i = 0; i < model->npts; i++) {
	    fprintf(stderr, "%d %d\n", i, ntf1[i]);
	    }
	fprintf(stderr, "\n\"nft%d\n", level);
	for (i = 0; i < model->npts; i++) {
	    fprintf(stderr, "%d %d\n", i, nft1[i]);
	    }
	fprintf(stderr, "\n\"nff%d\n", level);
	for (i = 0; i < model->npts; i++) {
	    fprintf(stderr, "%d %d\n", i, nff1[i]);
	    }
	}
	{
	GrayImage m;
	int i;
	char name[100];
	m = (GrayImage)imNew(IMAGE_GRAY, model->xsize, model->ysize);
	for (i = 0; i < model->npts; i++) {
	    imRef(m, model->pts[i].x, model->pts[i].y) =
		nbelow1[i] * (COLRNG - 1) / nprobes1;
	    }
	(void)sprintf(name, "/tmp/mod%d.pgm", level);
	(void)imSave(m, name);
	imFree(m);
	}
#endif
#endif

	/* That generated highpq for us - make it the new lowpq */
	assert(pqFirst(lowpq) == NULLPRIQNODE);
	pqFree(lowpq);
	lowpq = highpq;
	highpq = NULLPRIQ;
	}

    /*
     * We now have lowpq - the priority queue at highest resolution.
     * Read it out into trans.
     */
    for (pqn = pqFirst(lowpq); pqn != NULLPRIQNODE; pqn = pqFirst(lowpq)) {

	/* Get the hook from the first thing in the pq */
	hook = (spqHook *)pqGet(pqn);
	assert(hook != (spqHook *)NULL);