r-h.c

Hausdorff Distance for Image Recognition

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)
{
    int i;
    int x, y;
    long *vals = (long *)NULL;
    transval *newt;
    int model_required;
    int model_early;
    int model_early_required;
    int novermax;
    int valsunder;
    long minover;
    long pointval;
    boolean keep;
    
    /* Stuff cached from image and model */
    point *model_pts;
    unsigned nmodel_pts;
    unsigned image_plusx_width;
    LongImage image_dtrans;
    ShortImage image_xdtrans;
    long model_thresh;
    int maxx, maxy;
    
    assert(image != (image_info *)NULL);
    assert(model != (model_info *)NULL);
    assert(model->pts != (point *)NULL);
    assert(image->dtrans != (LongImage)NULL);
    assert(th != (TransHash *)NULL);

    DEBUG("sampling(%d..%d by %d, %d..%d by %d)\n", lowX, highX, stepX, lowY, highY, stepY);

    model_pts = model->pts;
    nmodel_pts = model->npts;
    model_thresh = model->model_thresh;
    maxx = model->xsize - 1;
    maxy = model->ysize - 1;

    image_dtrans = image->dtrans;

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

    model_early = (int)(ceil(nmodel_pts * EARLY_FRAC));
    model_early_required = (int)(ceil(model_early * model->model_frac));

    vals = (long *)malloc(nmodel_pts * sizeof(long));
    if (vals == (long *)NULL) {
	goto bailout;
	}
    
    /* Generate the skip-in-x transform for this threshold */
    if (!ensure_dtrans(image, thresh)) {
	goto bailout;
	}

    assert(image->plusx_dtrans != (ShortImage)NULL);
    image_xdtrans = image->plusx_dtrans;
    image_plusx_width = imGetWidth(image_xdtrans);

    /* Perform the scan. */
    for (y = lowY + stepY / 2; y < highY + stepY / 2; y += stepY) {
	/* Is there a place in this cell where the model fits in the borders? */
	if (y - stepY / 2 + maxy >= (int)(image->ysize + image->bottomborder)) {
	    /* No - even at the top of the cell it goes too far */
	    break;
	    }

	for (x = lowX + stepX / 2;
	     x < highX + stepX / 2;
	     x += stepX) {
	    
	    /* Can this fit in? */
	    if (x - stepX / 2 + maxx >= (int)(image->xsize + image->rightborder)) {
		/* No - even at the left of the cell it goes too far */
		break;
		}
	    
	    valsunder = 0;
	    minover = -1;
	    
	    for (i = 0; i < model_early; i++) {
		pointval = imRef(image_xdtrans,
				 model_pts[i].x + x, model_pts[i].y + y);
		if (pointval == 0) {
		    valsunder++;
		    }
		else {
		    if ((minover < 0) || (pointval < minover)) {
			minover = pointval;
			}
		    if (valsunder + (nmodel_pts - i - 1) <
			model_required) {
			/* DIE YOU GRAVY-SUCKING PIG! */
			break;
			}
		    }

		}
	    /*
	     * At this point, check valsunder for early decision
	     * Note that we aren't allowed to take early decisions at
	     * highest resolution.
	     */
	    if ((needdata) || (valsunder < model_early_required)) {
		/* No early decision - keep probing */
		for (i = model_early; i < nmodel_pts; i++) {
		    pointval = imRef(image_xdtrans,
				     model_pts[i].x + x, model_pts[i].y + y);
		    if (pointval == 0) {
			valsunder++;
			}
		    else {
			if ((minover < 0) || (pointval < minover)) {
			    minover = pointval;
			    }
			if (valsunder + (nmodel_pts - i - 1) <
			    model_required) {
			    /* DIE YOU GRAVY-SUCKING PIG! */
			    break;
			    }
			}
		    }
		if (i == nmodel_pts) {
		    keep = TRUE;
		    }
		else {
		    keep = FALSE;
		    }
		}
	    else {
		keep = TRUE;
		}

	    if (keep) {
		VDEBUG("Found one at (%d ", x, 10);
		VDEBUG("%d)\n", y, 10);
		
		/* Make a mark in th */
		if (needdata) {
		    /* Create a transval to hold the info */
		    newt = malloc_trans();
		    if (newt == (transval *)NULL) {
			goto bailout;
			}
		    newt->transpos.x = x;
		    newt->transpos.y = y;
		    
		    /* Work out the actual pointval */
		    for (i = 0; i < nmodel_pts; i++) {
			vals[i] = imRef(image_dtrans,
					model_pts[i].x + x, model_pts[i].y + y);
			}
		    pointval = pickNth(vals, (int)nmodel_pts,
				       model_required - 1);
		    newt->forward_val = pointval;
		    newt->forward_frac =
			frac_lower(vals, (int)nmodel_pts,
				   model_thresh);
		    }
		else {
		    newt = (transval *)NULL;
		    }
		
		/* The mark goes at the top left corner of the box */
		if (!thIsIn(th, x - stepX / 2, y - stepY / 2)) {
		    if (!thAdd(th, x - stepX / 2, y - stepY / 2,
			       (void *)newt)) {
			if (newt != (transval *)NULL) {
			    free((void *)newt);
			    }
			goto bailout;
			}
		    }
		}
	    else {
		/* Do some simple pruning */
		if (minover > image_plusx_width) {
		    /* Big skip */
		    break;
		    }
		
		/* Round up to the next point */
		if (stepX != 1) {
		    minover += stepX - 1;
		    minover -= minover % stepX;
		    }
		/* Precompensate for the loop increment */
		x += minover - stepX;
		
		}
	    }
	}
    
    /* Free now-useless stuff */
    free((void *)vals);

    /* and done. */
    return(TRUE);

  bailout:
    if (vals != (long *)NULL) {
	free((void *)vals);
	}

    return(FALSE);
    }

/*
 * Given lowth, a hash table containing interesting regions at one
 * resolution (lowX..highX by oldStepX etc), we want to compute
 * highth, a hash table containing interesting regions at another resolution
 * (lowX..highX by stepX etc), and store this in highth.
 *
 * We do this by finding "rectangular" regions in lowth which are
 * interesting, and calling findTransMtoISampled() on each, telling
 * it to store the results in highth. Once every region flagged as interesting
 * in lowth has been processed, we're done.
 * Each region is represented by its top left corner.
 * thresh defines what the threshold of "interesting" is.
 */
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)
{
    int botX, botY;
    int topX, topY;
    int x, y;
    boolean blockedX, blockedY;
    int nLow;
    int nLowTot = 0;
    int nProbes = 0;

    assert(image != (image_info *)NULL);
    assert(model != (model_info *)NULL);
    assert(lowth != (TransHash *)NULL);
    assert(highth != (TransHash *)NULL);

    while (thGetLowest(lowth, &x, &y)) {
	/* Expand this in increasing x, y as far as possible */
	topX = botX = x;
	topY = botY = y;

	/* We're not blocked in any dimension */
	blockedX = blockedY = FALSE;

	while ((!blockedX) || (!blockedY)) {
	    /* Expand in X as far as we can go */
	    while (!blockedX) {
		/* Try to expand one tick in X */
		if (topX + oldStepX >= highX) {
		    /* We'd go out of bounds */
		    blockedX = TRUE;
		    }
		else {
		    for (y = botY; y <= topY; y += oldStepY) {
			if (!thIsIn(lowth, topX + oldStepX, y)) {
			    /* Can't expand any more in X */
			    blockedX = TRUE;
			    break;
			    }
			}
		    }

		if (!blockedX) {
		    topX += oldStepX;
		    }
		}
	    
	    if (!blockedY) {
		/* Try to expand one tick in Y */
		if (topY + oldStepY >= highY) {
		    /* We'd go out of bounds */
		    blockedY = TRUE;
		    }
		else {
		    for (x = botX; x <= topX; x += oldStepX) {
			if (!thIsIn(lowth, x, topY + oldStepY)) {
			    /* Can't expand any more in Y */
			    blockedY = TRUE;
			    break;
			    }
			}
		    }

		if (!blockedY) {
		    topY += oldStepY;
		    }
		}
	    
	    }

	/*
	 * We now have a rectangular region, all of which was interesting,
	 * in (botX..topX)(botY..topY). Now,
	 * these are actually top-left-corner of cell coordinates. We
	 * want to end up with botX being the lowest X coordinate in the
	 * region that we want scanned, and topX being one more than the
	 * highest one: we want to check transformations which have X values
	 * in the botX..topX-1 range; find...Sampled will actually check
	 * ones spaced every stepX apart, starting at botX + stepX / 2.
	 *
	 * However, we must ensure that each resolution level has consistent
	 * coordinates: no matter where the current region of interest is,
	 * we want to have the higher-resolution cells we generate be
	 * on spacings which are consistent with all other high-resolution
	 * cells we may generate for other regions of interest.
	 */

	/* First, remove this region from lowth */
	nLow = 0;
	for (x = botX; x <= topX; x += oldStepX) {
	    for (y = botY; y <= topY; y += oldStepY) {
		nLow++;
		thRem(lowth, x, y);
		}
	    }

	/* Now compensate for the cell coordinate frame */
	/*
	 * Round down botX to be a multiple of stepX above lowX, and
	 * round topX to be a multiple of stepX above botX. This ensures
	 * that everything will line up with cells generated for other regions
	 * of interest.
	 * We want topX to be the first translation *not* considered at higher
	 * resolution, so we make it one cell higher than we would if
	 * it were the inclusive limit.
	 * We also must increase topX by oldStepX - 1, since that
	 * is the highest actual translation considered in this region of
	 * interest; we then round it down so that it is the left corner
	 * of the cell containing that translation.
	 */
    DEBUG("interest(%d..%d by %d, %d..%d by %d) = %d\n", botX, topX, stepX, botY, topY, stepY, nLow);
	botX -= (botX - lowX) % stepX;
	topX += oldStepX + stepX - 1;
	topX -= (topX - lowX) % stepX;

	botY -= (botY - lowY) % stepY;
	topY += oldStepY + stepY - 1;
	topY -= (topY - lowY) % stepY;

	if (!findTransMtoISampled(image, model,
				  botX, topX, stepX,
				  botY, topY, stepY,
				  highth, thresh, needdata)) {
	    /* Something failed down the line */
	    return(FALSE);
	    }
	nLowTot += nLow;
	nProbes++;

	}
    DEBUG("level covered %d cells in %d probes (%f avg)\n",
	    nLowTot, nProbes, nLowTot / (double)nProbes);

    return(TRUE);
    }

/*
 * Make sure that the image's +x distance transform exists and is based
 * on the correct threshold.
 */
static boolean
ensure_dtrans(image_info *image, long thresh)
{
    DEBUG("ensure_dtrans(%d)\n", thresh);

    if ((image->plusx_dtrans == (ShortImage)NULL) ||
	(thresh != image->dtrans_thresh)) {
	DEBUG("rebuild from %d\n", image->dtrans_thresh);
	/* Need to rebuild the +x dtrans */
	if (image->plusx_dtrans != (ShortImage)NULL) {
	    imFree(image->plusx_dtrans);
	    }

	image->plusx_dtrans =
	    dtrans_plusx(image->dtrans, thresh);

	if (image->plusx_dtrans == (ShortImage)NULL) {
	    return(FALSE);
	    }

	image->dtrans_thresh = thresh;

	}
    return(TRUE);
    }