utility.c

Hausdorff Distance for Image Recognition

 * stick them together to get a scaled model, which you call dtrans on.
 * Pad the dtrans by nbord pixels on each edge.
 * Return the dtrans if you can, NULL if not.
 *
 * Also, maintain the dtrans cache. This is basically the same as the
 * image dtrans cache, except that there are no fields in the smodel_info
 * structure to hold it - just return the dtrans. However, we want to
 * make sure that the cache holds only one entry for each (xs, ys) pair
 * (i.e. don't keep around entries that differ only by their nbord values).
 * There are things that rely on getting back a dtrans padded with exactly
 * nbord pixels, so don't return something with an nbord which is too large,
 * even though it has the right values. Sigh.
 *
 * N.B. The called *must not* try to free the value they get back from this.
 *
 * We want to avoid nasty states which could occur if something failed to
 * get allocated, but they keep on calling us.
 *
 * It tries to keep around a temp area so it doesn't keep allocating it.
 */
LongImage
get_smodel_dtrans_padded(smodel_info *model, int xs, int ys, int nbord)
{
    static point *scaledpts = (point *)NULL;
    int nscaledpts = 0;

    int i;
    int slot;
    int curw, curh;

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

    if (model->modeldtcache == (smodeldt_info *)NULL) {
	/* Allocate it */
	assert(model->ncached == 0);

	model->modeldtcache =
	    (smodeldt_info *)malloc(MAXCACHE * sizeof(smodeldt_info));
	if (model->modeldtcache == (smodeldt_info *)NULL) {
	    return((LongImage)NULL);
	    }
	}
				
    /* Search the cache for it. */
    slot = -1;
    for (i = 0; i < model->ncached; i++) {
	if ((model->modeldtcache[i].xs == xs) &&
	    (model->modeldtcache[i].ys == ys) &&
	    (model->modeldtcache[i].dtrans != (LongImage)NULL)) {
	    /* Found it, if the borders match. */
	    if (model->modeldtcache[i].nbord == nbord) {
		return(model->modeldtcache[i].dtrans);
		}
	    else {
		slot = i;
		break;
		}
	    }
	}

    /* Not present in the cache. */
    if (slot < 0) {
	/* Find a slot, if we didn't have an all-but-borders match */
	if (model->ncached < MAXCACHE) {
	    /* Allocate a new slot. */
	    slot = model->ncached++;
	    }
	else {
	    /* Kill an old slot */
	    slot = random() % MAXCACHE;
	    if (model->modeldtcache[slot].dtrans != (LongImage)NULL) {
		imFree(model->modeldtcache[slot].dtrans);
		}
	    }
	}
    else {
	/*
	 * Toast the old dtrans.
	 * Since the only mismatch is with the borders, if the current
	 * nbord is smaller than the previous one, then it should be possible
	 * to just trim off the extra bits, but I don't feel like implementing
	 * this now. FIX.
	 */
	if (model->modeldtcache[slot].dtrans != (LongImage)NULL) {
	    imFree(model->modeldtcache[slot].dtrans);
	    }
	}

    /* Clear the slot we're about to write into, in case of failure. */
    model->modeldtcache[slot].dtrans = (LongImage)NULL;

    /* Make sure we have the appropriate scaled points */
    if (!ensure_xscales(model, xs)) {
	return((LongImage)NULL);
	}
    if (!ensure_yscales(model, ys)) {
	return((LongImage)NULL);
	}

    /* Build somewhere we can put the scaled model */
    if ((scaledpts == (point *)NULL) || (nscaledpts < model->npts)) {
	if (scaledpts != (point *)NULL) {
	    free((void *)scaledpts);
	    }
	scaledpts = (point *)malloc(sizeof(*scaledpts) * model->npts);
	if (scaledpts == (point *)NULL) {
	    return((LongImage)NULL);
	    }
	nscaledpts = model->npts;
	}

    assert(model->scaled_x[xs] != (int *)NULL);
    assert(model->scaled_y[ys] != (int *)NULL);
    
    /* Build the model pointlist */
    for (i = 0; i < model->npts; i++) {
	scaledpts[i].x = model->scaled_x[xs][i];
	scaledpts[i].y = model->scaled_y[ys][i];
	assert(scaledpts[i].x <= xs);
	assert(scaledpts[i].y <= ys);
	}

    /*
     * Work out how big a box in the image we might need to scan
     * (i.e. how big to make the model dtrans)
     */
    curw = (int)((model->xsize - 1) * xs * model->scaleStepX + 0.5) + 1;
    curh = (int)((model->ysize - 1) * ys * model->scaleStepY + 0.5) + 1;
    
    /* Now make up the model distance transform, adding in the borders */
    model->modeldtcache[slot].dtrans =
	dtrans_pts((unsigned)curw + 2 * nbord, (unsigned)curh + 2 * nbord,
		   -nbord, -nbord, model->npts, scaledpts);

    if (model->modeldtcache[slot].dtrans == (LongImage)NULL) {
	return((LongImage)NULL);
	}

    /* Make sure the cache tags are right */
    model->modeldtcache[slot].xs = xs;
    model->modeldtcache[slot].ys = ys;
    model->modeldtcache[slot].nbord = nbord;

    return(model->modeldtcache[slot].dtrans);
    }

/*
 * Make sure that the model's scaled points array is filled in for
 * this scale. Generate or extend scaled_x if necessary.
 */
boolean
ensure_xscales(smodel_info *model, int xscale)
{
    int i;
    int *p;
    double scale;
    int **oscaled_x;

    assert(model != (smodel_info *)NULL);
    assert(model->pts != (point *)NULL);
    /* We need maxx and maxy; compute_scaled_pts must have been called */
    assert(xscale >= 0);

    if (model->scaled_x == (int **)NULL) {
	/* It hasn't been allocated yet. Do so. */
	/* Allocate the pointer blocks for the scaled points */
	model->scaled_x = (int **)malloc((unsigned)(xscale + 1) *
					 sizeof(int *));
	if (model->scaled_x == (int **)NULL) {
	    goto bailout;
	    }
	
	for (i = 0; i <= xscale; i++) {
	    model->scaled_x[i] = (int *)NULL;
	    }
	model->nscaled_x = xscale + 1;
	}
    else if (model->nscaled_x <= xscale) {
	/* Extend it. */

	oscaled_x = model->scaled_x;

	/* Generate the new block */
	model->scaled_x = (int **)malloc((unsigned)(xscale + 1) *
					 sizeof(int *));
	if (model->scaled_x == (int **)NULL) {
	    model->scaled_x = oscaled_x;
	    goto bailout;
	    }
	
	/* Clear the new block */
	for (i = 0; i <= xscale; i++) {
	    model->scaled_x[i] = (int *)NULL;
	    }
	/* and copy the old block */
	for (i = 0; i < model->nscaled_x; i++) {
	    model->scaled_x[i] = oscaled_x[i];
	    }
	/* Free the old block */
	free((void *)oscaled_x);
	model->nscaled_x = xscale + 1;
	}

    if (model->scaled_x[xscale] != (int *)NULL) {
	/* Already been calculated */
	return(TRUE);
	}
    
    model->scaled_x[xscale] = (int *)malloc(model->npts * sizeof(int));
    if (model->scaled_x[xscale] == (int *)NULL) {
	goto bailout;
	}

    /* Compute the things scaled by X */
    scale = xscale * model->scaleStepX;
    for (i = 0, p = model->scaled_x[xscale]; i < model->npts; i++, p++) {
	*p = (int)(model->pts[i].x * scale + 0.5);
	}

    return(TRUE);

  bailout:

    return(FALSE);
    }

/*
 * Make sure that the model's scaled points array is filled in for
 * this scale.
 */
boolean
ensure_yscales(smodel_info *model, int yscale)
{
    int i;
    int *p;
    double scale;
    int **oscaled_y;

    assert(model != (smodel_info *)NULL);
    assert(model->pts != (point *)NULL);
    /* We need maxx and maxy; compute_scaled_pts must have been called */
    assert(yscale >= 0);

    if (model->scaled_y == (int **)NULL) {
	/* It hasn't been allocated yet. Do so. */
	/* Allocate the pointer blocks for the scaled points */
	model->scaled_y = (int **)malloc((unsigned)(yscale + 1) *
					 sizeof(int *));
	if (model->scaled_y == (int **)NULL) {
	    goto bailout;
	    }
	
	for (i = 0; i <= yscale; i++) {
	    model->scaled_y[i] = (int *)NULL;
	    }
	model->nscaled_y = yscale + 1;
	}
    else if (model->nscaled_y <= yscale) {
	/* Extend it. */

	oscaled_y = model->scaled_y;

	/* Generate the new block */
	model->scaled_y = (int **)malloc((unsigned)(yscale + 1) *
					 sizeof(int *));
	if (model->scaled_y == (int **)NULL) {
	    model->scaled_y = oscaled_y;
	    goto bailout;
	    }
	
	/* Clear the new block */
	for (i = 0; i <= yscale; i++) {
	    model->scaled_y[i] = (int *)NULL;
	    }
	/* and copy the old block */
	for (i = 0; i < model->nscaled_y; i++) {
	    model->scaled_y[i] = oscaled_y[i];
	    }
	/* Free the old block */
	free((void *)oscaled_y);
	model->nscaled_y = yscale + 1;
	}

    if (model->scaled_y[yscale] != (int *)NULL) {
	/* Already been calculated */
	return(TRUE);
	}
    
    model->scaled_y[yscale] = (int *)malloc(model->npts * sizeof(int));
    if (model->scaled_y[yscale] == (int *)NULL) {
	goto bailout;
	}

    /* Compute the things scaled by Y */
    scale = yscale * model->scaleStepY;
    for (i = 0, p = model->scaled_y[yscale]; i < model->npts; i++, p++) {
	*p = (int)(model->pts[i].y * scale + 0.5);
	}

    return(TRUE);

  bailout:

    return(FALSE);
    }