scale-h.c

Hausdorff Distance for Image Recognition

	/* and unhook it */
	pqRem(pqn);

	/* The values must be valid */
	assert(hook->model_thresh == model->model_thresh);
	assert(hook->model_frac == model->model_frac);

	newt = malloc_strans();
	if (newt == (stransval *)NULL) {
	    free_spqHook(hook);		/* Free it - we own it now */
	    goto bailout;
	    }
	newt->transpos.x = hook->x;
	newt->transpos.y = hook->y;
	newt->scalex = hook->scaleX * model->scaleStepX;
	newt->scaley = hook->scaleY * model->scaleStepY;
	newt->forward_val = hook->forward_val;
	newt->forward_frac = hook->forward_frac;
	newt->forward_frac2 = hook->forward_frac2;
	/* The hook is no longer needed */
	free_spqHook(hook);

	newln = liAdd(trans, (void *)newt);
	if (newln == NULLLISTNODE) {
	    free_strans(newt);
	    goto bailout;
	    }

	}

    /* Free useless stuff */
    assert(pqFirst(lowpq) == NULLPRIQNODE);
    pqFree(lowpq);

    model->trans = trans;

    
#ifdef INST2
    {
    int i;
    fprintf(stderr, "\"ntt\n");
    for (i = 0; i < model->npts; i++) {
	fprintf(stderr, "%d %d\n", i, ntt[i]);
	}
    fprintf(stderr, "\n\"ntf\n");
    for (i = 0; i < model->npts; i++) {
	fprintf(stderr, "%d %d\n", i, ntf[i]);
	}
    fprintf(stderr, "\n\"nft\n");
    for (i = 0; i < model->npts; i++) {
	fprintf(stderr, "%d %d\n", i, nft[i]);
	}
    fprintf(stderr, "\n\"nff\n");
    for (i = 0; i < model->npts; i++) {
	fprintf(stderr, "%d %d\n", i, nff[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) =
	    nbelow[i] * (COLRNG - 1) / nprobes;
	}
    (void)sprintf(name, "/tmp/modall.pgm");
    (void)imSave(m, name);
    imFree(m);
    }
#endif

    return;

  bailout:

    if (lowpq != NULLPRIQ) {
	free_spqHooklist(lowpq);
	}

    if (highpq != NULLPRIQ) {
	free_spqHooklist(highpq);
	}

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

    model->trans = NULLLIST;
    }

/*
 * Check the transformations that we were handed (in model->trans); remove
 * those where the thresholds aren't satisfied.
 */
static void
checkSTransMtoI(simage_info *image, smodel_info *model)
{
    ListNode ln;
    long pointval = 0;
    long *vals = (long *)NULL;
    int *curxvals;
    int *curyvals;
    stransval *t;
    int x, y;
    int xs, ys;
    int i;
    int model_required;

    int nover;
    int novermax;

    /* Stuff cached from image and model */
    unsigned nmodel_pts;
    LongImage image_dtrans;
    long model_thresh;
    double model_frac;
    double scalestepx, scalestepy;
    double minscalex, minscaley;
    double maxskew;
    
    assert(image != (simage_info *)NULL);
    assert(model != (smodel_info *)NULL);
    assert(image->dtrans != (LongImage)NULL);
    assert(model->trans != NULLLIST);

    nmodel_pts = model->npts;
    model_thresh = model->model_thresh;
    model_frac = model->model_frac;
    image_dtrans = image->dtrans;
    minscalex = model->minscalex;
    minscaley = model->minscaley;
    maxskew = model->maxskew;
    scalestepx = model->scaleStepX;
    scalestepy = model->scaleStepY;

    model_required = (int)(ceil(nmodel_pts * (model_frac - EPSILON)));
    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 = (stransval *)liGet(ln);
	assert(t != (stransval *)NULL);
	x = t->transpos.x;
	y = t->transpos.y;
	/* Convert the X and Y scale values to integers... */
	xs = (int)(t->scalex / scalestepx + 0.5);
	ys = (int)(t->scaley / scalestepy + 0.5);

	if (!ensure_xscales(model, xs)) {
	    goto bailout;
	    }
	if (!ensure_yscales(model, ys)) {
	    goto bailout;
	    }
	curxvals = model->scaled_x[xs];
	curyvals = model->scaled_y[ys];

	if (model_frac == 1.) {
	    pointval = 0;
	    for (i = 0; i < nmodel_pts; i++) {
		pointval = MAX(pointval, imRef(image_dtrans,
					       curxvals[i] + x,
					       curyvals[i] + y));
		if (pointval > model_thresh) {
		    break;
		    }
		}
	    }
	else {
	    nover = 0;
	    for (i = 0; i < nmodel_pts; i++) {
		/* Probe at this model point's transformed location */
		pointval = vals[i] = imRef(image_dtrans,
					   curxvals[i] + x,
					   curyvals[i] + 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.;

		if (pointval == model_thresh) {
		    t->forward_frac2 = 1.;
		    }
		else {
		    /* Great - we have to rescan */
		    nover = 0;
		    for (i = 0; i < nmodel_pts; i++) {
			if (imRef(image_dtrans,
				  curxvals[i] + x,
				  curyvals[i] + y) > pointval) {
			    nover++;
			    }
			}
		    t->forward_frac2 =
			(nmodel_pts - nover) / ((double)nmodel_pts);
		    }

		}
	    else {
		t->forward_frac = frac_lower(vals, (int)nmodel_pts,
					     model_thresh);
		assert(t->forward_frac >= model_frac - EPSILON);
		if (pointval == model_thresh) {
		    t->forward_frac2 = t->forward_frac;
		    }
		else {
		    t->forward_frac2 =
			frac_lower(vals, (int)nmodel_pts, pointval);
		    }
		
		}

	    ln = liNext(ln);
	    }
	else {
	    free_strans(t);
	    ln = liRemAndNext(ln);
	    }
	}

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

    return;

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

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

/*
 * This routine evaluates the forward distance for a single stransval.
 * It returns TRUE on success, FALSE on failure. It does no checking
 * of ranges etc - just evaluates the stransval wrt the forward match
 * parameters. It fills in the forward_val and forward_frac* fields.
 */
static boolean
checkOneSTransMtoI(simage_info *image, smodel_info *model, stransval *t)
{
    int x, y;
    int xs, ys;

    assert(image != (simage_info *)NULL);
    assert(model != (smodel_info *)NULL);
    assert(t != (stransval *)NULL);
    
    x = t->transpos.x;
    y = t->transpos.y;
    /* Convert the X and Y scale values to integers... */
    xs = (int)(t->scalex / model->scaleStepX + 0.5);
    ys = (int)(t->scaley / model->scaleStepY + 0.5);

    /* and evaluate, storing into t */
    return(evalOneSTransMtoI(image, model, x, y, xs, ys,
			     &(t->forward_val), &(t->forward_frac),
			     &(t->forward_frac2)));
    }

/*
 * This routine evaluates the forward distance for a single location.
 * It returns TRUE on success, FALSE on failure. It does no checking
 * of ranges etc - just evaluates the location wrt the forward match
 * parameters.
 */
static boolean
evalOneSTransMtoI(simage_info *image, smodel_info *model,
		  int x, int y, int xs, int ys,
		  long *dist, double *frac, double *frac2)
{
    static long *vals = (long *)NULL;
    static int nvals = 0;

    int *curxvals;
    int *curyvals;
    int i;
    int model_required;

    /* Stuff cached from image and model */
    LongImage image_dtrans;
    unsigned nmodel_pts;

    assert(image != (simage_info *)NULL);
    assert(model != (smodel_info *)NULL);
    assert(dist != (long *)NULL);
    assert(frac != (double *)NULL);
    assert(frac2 != (double *)NULL);
    
    assert(image->dtrans != (LongImage)NULL);
    image_dtrans = image->dtrans;
    nmodel_pts = model->npts;

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

    /* Make sure our work area is valid */
    if ((vals == (long *)NULL) || (nvals < nmodel_pts)) {
	if (vals != (long *)NULL) {
	    free((void *)vals);
	    }
	vals = (long *)malloc(nmodel_pts * sizeof(*vals));
	if (vals == (long *)NULL) {
	    goto bailout;
	    }
	nvals = nmodel_pts;
	}

    /* Get the scaled pts */
    if (!ensure_xscales(model, xs)) {
	goto bailout;
	}
    if (!ensure_yscales(model, ys)) {
	goto bailout;
	}
    curxvals = model->scaled_x[xs];
    curyvals = model->scaled_y[ys];

    /* Probe at this model point's transformed location */
    for (i = 0; i < nmodel_pts; i++) {
	vals[i] = imRef(image_dtrans, curxvals[i] + x, curyvals[i] + y);
	}

    /* and evaluate */
    *dist = pickNth(vals, (int)nmodel_pts, model_required - 1);
    *frac = frac_lower(vals, (int)nmodel_pts, model->model_thresh);
    if (*dist == model->model_thresh) {
	*frac2 = *frac;
	}
    else {
	*frac2 = frac_lower(vals, (int)nmodel_pts, *dist);
	}

    return(TRUE);

  bailout:
    return(FALSE);

    }

/*
 * This routine scans the transformations listed in modeltrans, looking for
 * those where at least image_frac of the image points under the model lie
 * within image_thresh of model points.
 *
 * It sorts those transformations by their scale parameters, since each
 * distinct set of scale values requires a new model to be constructed and
 * dtransed; this is expensive if there are lots of transformations to be
 * checked. We therefore group them and reuse the dtranses wherever possible.
 *
 * Note that this means that modeltrans is actually rebuilt as a new list
 * in a possibly different order.
 */
static void
checkSTransItoM(simage_info *image, smodel_info *model, int revstyle)
{
    static long *vals = (long *)NULL;
    static int nstaticvals = 0;

    ListNode ln;
    ListNode newln;
    stransval *t;
    int x, y;
    int x2, y2;
    int xs, ys;
    int curw, curh;
    int minx, maxx;
    int miny, maxy;
    long pointval;
    int nvals;
    int whichval;
    LongImage image_dtrans;
    LongImage model_dtrans;
    PtrImage scaleim = (PtrImage)NULL;
    List slist;
    int nimmax;
    int nimage;
    int maxxscale;
    int maxyscale;

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

    assert((revstyle == REVERSE_BOX) || (revstyle == REVERSE_ALLIMAGE));

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

    /* Figure out how many image pixels might be visible */
    if (revstyle == REVERSE_BOX) {
	nimmax = MIN(nimage, model->xsize * model->ysize);
	}
    else {
	nimmax = nimage;
	}

    /* Figure out how big an array we need */
    maxxscale = floor(model->maxscalex / model->scaleStepX + EPSILON) + 1;
    maxyscale = floor(model->maxscaley / model->scaleStepY + EPSILON) + 1;

    if ((vals == (long *)NULL) || (nstaticvals < nimmax)) {
	if (vals != (long *)NULL) {
	    free((void *)vals);
	    }
	vals = (long *)malloc(nimmax * sizeof(long));
	if (vals == (long *)NULL) {
	    goto bailout;
	    }
	nstaticvals = nimmax;
	}

    /* Build scaleim - bucket sorter for the scales */
    scaleim = (PtrImage)imNew(IMAGE_PTR,
			      (unsigned)maxxscale + 1,
			      (unsigned)maxyscale + 1);
    if (scaleim  == (PtrImage)NULL) {
	goto bailout;
	}

    /*
     * Read the list of translations out into scaleim.