scale-h.c

Hausdorff Distance for Image Recognition

     * Destroy model->trans as you go, so that there are never two
     * ListNodes which hold the same stransval. This makes cleanup
     * in the case of failure (bailout) possible, since otherwise we
     * wouldn't be able to tell which bits of model->trans had been read
     * out and which hadn't.
     */
    for (ln = liFirst(model->trans);
	 ln != NULLLISTNODE;
	 ln = liRemAndNext(ln)) {
	t = (stransval *)liGet(ln);
	assert(t != (stransval *)NULL);

	/* Compute the scale space grid coordinates */
	xs = (int)(t->scalex / model->scaleStepX + 0.5);
	ys = (int)(t->scaley / model->scaleStepY + 0.5);
	assert(xs >= 0);
	assert(ys >= 0);
	assert(xs <= maxxscale);
	assert(ys <= maxyscale);
	    
	if (imRef(scaleim, xs, ys) == (void *)NULL) {
	    /* First one at this scale */
	    slist = liNew();
	    if (slist == NULLLIST) {
		goto bailout;
		}
	    imRef(scaleim, xs, ys) = (void *)slist;
	    }
	/* Hook it into the bin */
	newln = liAdd((List)imRef(scaleim, xs, ys), (void *)t);
	if (newln == NULLLISTNODE) {
	    goto bailout;
	    }
	}

    /*
     * We don't need model->trans any more - free it.
     * Note that we *do not* use free_stranslist, since that would free
     * all the stransvals, which we just hooked into scaleim. Actually,
     * at this point model->trans should be empty.
     */
    assert(liLen(model->trans) == 0);
    liFree(model->trans);
    model->trans = NULLLIST;

    /* Now work on the bins in scaleim. */
    for (ys = 0; ys <= maxyscale; ys++) {
	for (xs = 0; xs <= maxxscale; xs++) {
	    if (imRef(scaleim, xs, ys) != (void *)NULL) {
		slist = (List)imRef(scaleim, xs, ys);

		/* Build the dtrans for scalings of xs, ys */
		if (revstyle == REVERSE_BOX) {
		    model_dtrans = get_smodel_dtrans(model, xs, ys);
		    }
		else {
		    assert(revstyle == REVERSE_ALLIMAGE);
		    /*
		     * Put borders of model->image_thresh on the dtrans.
		     * This means that if an image pixel's value is below
		     * thresh, then we will get its real value; any image
		     * pixel outside the model trans would have been over
		     * thresh if we'd seen it; we just need to keep track
		     * of how many such pixels there are.
		     */
		    model_dtrans =
			get_smodel_dtrans_padded(model, xs, ys,
						 (int)ceil(model->image_thresh /
						      (double)DSCALE));
		    }
		if (model_dtrans == (LongImage)NULL) {
		    goto bailout;
		    }
		
		/* Find out the area we need to scan */
		curw = imGetWidth(model_dtrans);
		curh = imGetHeight(model_dtrans);

		/* Now work on the list */
		for (ln = liFirst(slist); ln != NULLLISTNODE; ) {
		    t = (stransval *)liGet(ln);
		    assert(t != (stransval *)NULL);
		    x = t->transpos.x;
		    y = t->transpos.y;
		    
		    /* Scan over the place in the image where this lands */
		    nvals = 0;
		    if (revstyle == REVERSE_BOX) {
			assert(imGetXBase(model_dtrans) == 0);
			assert(imGetYBase(model_dtrans) == 0);
			/* We know that this lies within borders */
			for (y2 = y; y2 < y + curh; y2++) {
			    for (x2 = x; x2 < x + curw; x2++) {
				if (imRef(image_dtrans, x2, y2) == 0) {
				    vals[nvals++] =
					imRef(model_dtrans, x2 - x, y2 - y);
				    }
				}
			    }
			}
		    else {
			/* We have to do border checks */
			minx = MAX(x + imGetXBase(model_dtrans),
				   imGetXBase(image_dtrans));
			miny = MAX(y + imGetYBase(model_dtrans),
				   imGetYBase(image_dtrans));
			maxx = MIN(x + imGetXMax(model_dtrans),
				   imGetXMax(image_dtrans));
			maxy = MIN(y + imGetYMax(model_dtrans),
				   imGetYMax(image_dtrans));
			
			for (y2 = miny; y2 <= maxy; y2++) {
			    for (x2 = minx; x2 <= maxx; x2++) {
				if (imRef(image_dtrans, x2, y2) == 0) {
				    vals[nvals++] =
					imRef(model_dtrans, x2 - x, y2 - y);
				    }
				}
			    }
			}

		    if (revstyle == REVERSE_BOX) {
			whichval =
			    (int)(ceil(model->image_frac * (nvals - EPSILON)));
			pointval = pickNth(vals, nvals, whichval - 1);
			}
		    else {
			/*
			 * 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.;
				t->reverse_frac2 = 1.;
				}
			    else {
				t->reverse_frac =
				    frac_lower(vals, nvals,
					       model->image_thresh);
				if (pointval == model->image_thresh) {
				    t->reverse_frac2 = t->reverse_frac;
				    }
				else {
				    t->reverse_frac2 =
					frac_lower(vals, nvals, pointval);
				    }
				}
			    t->reverse_num = nvals;
			    }
			else {
			    if (nvals == 0) {
				if (nimage == 0) {
				    t->reverse_frac = 1.;
				    t->reverse_frac2 = 1.;
				    }
				else {
				    t->reverse_frac = 0.;
				    t->reverse_frac2 = 0.;
				    }
				}
			    else {
				t->reverse_frac =
				    frac_lower(vals, nvals,
					       model->image_thresh) * nvals /
						   nimage;
				if (pointval == model->image_thresh) {
				    t->reverse_frac2 = t->reverse_frac;
				    }
				else {
				    t->reverse_frac2 =
					frac_lower(vals, nvals,
						   pointval) * nvals /
						       nimage;
				    }
				}
			    
			    t->reverse_num = nimage;
			    }

			assert(t->reverse_frac >= model->image_frac - EPSILON);
			
			/* and advance */
			ln = liNext(ln);
			}
		    else {
			/* Free the hook from slist */
			free_strans(t);
			ln = liRemAndNext(ln);
			}
		    
		    }
		}
	    }
	}

    /*
     * The stransvals which remain in scaleim have been validated
     * and filled out. Read them out into model->trans.
     */
    model->trans = liNew();
    if (model->trans == NULLLIST) {
	goto bailout;
	}

    for (ys = 0; ys <= maxyscale; ys++) {
	for (xs = 0; xs <= maxxscale; xs++) {
	    if (imRef(scaleim, xs, ys) != (void *)NULL) {
		slist = (List)imRef(scaleim, xs, ys);
		liApp(model->trans, slist);
		imRef(scaleim, xs, ys) = (void *)NULL;
		}
	    }
	}

    imFree(scaleim);
    return;

  bailout:
    /* Don't free model_dtrans - we don't own it */

    if (model->trans != NULLLIST) {
	free_stranslist(model->trans);
	model->trans = NULLLIST;
	}

    if (scaleim != (PtrImage)NULL) {
	for (ys = 0; ys <= maxyscale; ys++) {
	    for (xs = 0; xs <= maxxscale; xs++) {
		if (imRef(scaleim, xs, ys) != (void *)NULL) {
		    slist = (List)imRef(scaleim, xs, ys);
		    free_stranslist(slist);
		    }
		}
	    }
	imFree(scaleim);
	}

    return;
    }

/*
 * Evaluate the transformation (x,y,xs,ys) for model wrt image, using
 * the supplied matching parameters and revstyle.
 * Return TRUE in *result if it checks out, FALSE if not.
 * Return TRUE on success, FALSE on failure.
 */
static boolean
checkOneSTransItoM(simage_info *image, smodel_info *model,
		   long reverse_thresh, double reverse_frac, int revstyle,
		   int x, int y, int xs, int ys, boolean *result)
{
    int curw, curh;
    int minx, maxx;
    int miny, maxy;
    int nvals;
    int nunder;
    LongImage image_dtrans;
    LongImage model_dtrans;
    int x2, y2;

    assert(image != (simage_info *)NULL);
    assert(model != (smodel_info *)NULL);
    assert(image->dtrans != (LongImage)NULL);
    assert(result != (boolean *)NULL);

    assert(xs >= 0);
    assert(ys >= 0);

    if (revstyle == REVERSE_NONE) {
	/* All succeed */
	*result = TRUE;
	return(TRUE);
	}

    /* At the moment, it has to be this */
    assert((revstyle == REVERSE_BOX) || (revstyle == REVERSE_ALLIMAGE));

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

    /* Get the model dtrans */
    if (revstyle == REVERSE_BOX) {
	model_dtrans = get_smodel_dtrans(model, xs, ys);
	}
    else {
	model_dtrans = get_smodel_dtrans_padded(model, xs, ys,
						(int)ceil(reverse_thresh /
							  (double)DSCALE));
	}

    if (model_dtrans == (LongImage)NULL) {
	goto bailout;
	}

    curw = imGetWidth(model_dtrans);
    curh = imGetHeight(model_dtrans);

    /* Scan over the place in the image where this lands */
    nvals = nunder = 0;
    if (revstyle == REVERSE_BOX) {
	/* We know that this lies within borders */
	for (y2 = y; y2 < y + curh; y2++) {
	    for (x2 = x; x2 < x + curw; x2++) {
		if (imRef(image_dtrans, x2, y2) == 0) {
		    nvals++;
		    if (imRef(model_dtrans, x2 - x, y2 - y) <= reverse_thresh) {
			nunder++;
			}
		    }
		}
	    }
	}
    else {
	/* We have to do border checks */
	minx = MAX(x + imGetXBase(model_dtrans),
		   imGetXBase(image_dtrans));
	miny = MAX(y + imGetYBase(model_dtrans),
		   imGetYBase(image_dtrans));
	maxx = MIN(x + imGetXMax(model_dtrans),
		   imGetXMax(image_dtrans));
	maxy = MIN(y + imGetYMax(model_dtrans),
		   imGetYMax(image_dtrans));
	
	for (y2 = miny; y2 <= maxy; y2++) {
	    for (x2 = minx; x2 <= maxx; x2++) {
		if (imRef(image_dtrans, x2, y2) == 0) {
		    if (imRef(model_dtrans, x2 - x, y2 - y) <= reverse_thresh) {
			nunder++;
			}
		    }
		}
	    }

	nvals = image->npts;
	}

    if (nunder >= nvals * reverse_frac) {
	*result = TRUE;
	}
    else {
	*result = FALSE;
	}

    return(TRUE);

  bailout:

    /* Don't free model_dtrans - we don't own it */

    return(FALSE);

    }

/* Do the forwstyle == FORWARD_ONE matching */
static void
findSTransOne(simage_info *image, smodel_info *model, int revstyle)
{
    BFS(image, model, revstyle, findOneCompare, findOneMatch);
    /* We only wanted one match, and it's been thoroughly checked out. */
    }

static void
findSTransBestDist(simage_info *image, smodel_info *model, int revstyle)
{
    ListNode ln;
    stransval *t;

    assert(model != (smodel_info *)NULL);
    assert(image != (simage_info *)NULL);

    /* Set the model distance to a known large value, if it is < 0 */
    if (model->model_thresh < 0) {
	model->model_thresh =
	    DSCALE * (image->xsize + image->ysize +
		      image->leftborder + image->topborder +
		      image->rightborder + image->bottomborder);
	}
    /* Otherwise, it's set to a good starting value */

    BFS(image, model, revstyle, findDistCompare, findDistMatch);

    if (model->trans == NULLLIST) {
	/* Failed */
	return;
	}

    /*
     * model->trans now has a bunch of matches, some of which may not
     * meet the specs, due to the search method: whenever we find a match
     * which is better than the previous best, we throw it on to model->trans,
     * and rely on this cleanup phase to remove the previous best.
     *
     * Also, the forward_frac values won't be right: they will be based on
     * whatever the threshold was set to when this match was generated; it
     * may have decreased since then. However, the frac2 will be exactly
     * what we want; they'll all be based on the right model_thresh, since
     * all the matches we want to keep have the same forward_val.
     */

    for (ln = liFirst(model->trans); ln != NULLLISTNODE; ) {
	t = (stransval *)liGet(ln);
	assert(t != (stransval *)NULL);

	assert(t->forward_val >= model->model_thresh);
	if (t->forward_val > model->model_thresh) {
	    /* This one didn't work */
	    free_strans(t);
	    ln = liRemAndNext(ln);
	    }
	else {
	    /* t->forward_frac may be out of date - clobber it */
	    t->forward_frac = t->forward_frac2;
	    ln = liNext(ln);
	    }
	}

    }

static void
findSTransBestFrac(simage_info *image, smodel_info *model, int revstyle)
{
    ListNode ln;
    stransval *t;

    assert(model != (smodel_info *)NULL);
    assert(image != (simage_info *)NULL);

    /* Set the model fraction to a known small value, if it is < 0 */
    if (model->model_frac < 0) {
	model->model_frac = 0.;
	}
    /* Otherwise, it's set to a good starting value */