dtrans.c

Hausdorff Distance for Image Recognition

	for (x = xbase + (int)width - 1; x >= xbase; x--) {
	    minval = imRef(im, x, y);
	    minval = MIN(minval, imRef(im, x, y + 1) + ONE);
	    if (x <= xbase + (int)width - 2) {
		minval = MIN(minval, imRef(im, x + 1, y) + ONE);
		minval = MIN(minval, imRef(im, x + 1, y + 1) + ROOT_2);
		}
	    if (x <= xbase + (int)width - 3) {
		minval = MIN(minval, imRef(im, x + 2, y + 1) + ROOT_5);
		}
	    if (x > xbase) {
		minval = MIN(minval, imRef(im, x - 1, y + 1) + ROOT_2);
		}
	    if (x > xbase + 1) {
		minval = MIN(minval, imRef(im, x - 2, y + 1) + ROOT_5);
		}
	    imRef(im, x, y) = minval;
	    }
	}

    /* Other lines */
    for (y = ybase + (int)height - 3; y >= ybase; y--) {
	/* Last entry */
	x = xbase + (int)width - 1;
	if (width >= 1) {
	    minval = imRef(im, x, y);
	    minval = MIN(minval, imRef(im, x, y + 1) + ONE);
	    if (width >= 2) {
		minval = MIN(minval, imRef(im, x - 1, y + 1) + ROOT_2);
		minval = MIN(minval, imRef(im, x - 1, y + 2) + ROOT_5);
		}
	    if (width >= 3) {
		minval = MIN(minval, imRef(im, x - 2, y + 1) + ROOT_5);
		}
	    imRef(im, x, y) = minval;
	    }

	/* Second last entry */
	if (width >= 2) {
	    x = xbase + (int)width - 2;
	    minval = imRef(im, x, y);
	    minval = MIN(minval, imRef(im, x, y + 1) + ONE);
	    minval = MIN(minval, imRef(im, x + 1, y) + ONE);
	    minval = MIN(minval, imRef(im, x + 1, y + 1) + ROOT_2);
	    minval = MIN(minval, imRef(im, x + 1, y + 2) + ROOT_5);
	    if (width >= 3) {
		minval = MIN(minval, imRef(im, x - 1, y + 1) + ROOT_2);
		minval = MIN(minval, imRef(im, x - 1, y + 2) + ROOT_5);
		}
	    if (width >= 4) {
		minval = MIN(minval, imRef(im, x - 2, y + 1) + ROOT_5);
		}
	    imRef(im, x, y) = minval;
	    }
	    
	if (width >= 5) {
	    long xp1y, xy;
	    long xp2yp1, xp1yp1, xyp1, xm1yp1, xm2yp1;
	    long xp1yp2, xyp2, xm1yp2;
	    LongPixPtr xyp, xm2yp1p, xm1yp2p;
	    LongPixPtr xlim;
	    
	    /* Start it out */
	    x = xbase + (int)width - 3;

	    xp1y = imRef(im, x + 1, y) + ONE;
	    xyp = imGetPixPtr(im, x, y);
	    xy = imPtrRef(im, xyp);
	    
	    xp2yp1 = imRef(im, x + 2, y + 1) + ROOT_5;
	    xp1yp1 = imRef(im, x + 1, y + 1) + ROOT_2;
	    xyp1 = imRef(im, x, y + 1) + ONE;
	    xm1yp1 = imRef(im, x - 1, y + 1) + ROOT_2;
	    xm2yp1p = imGetPixPtr(im, x - 2, y + 1);
	    xm2yp1 = imPtrRef(im, xm2yp1p) + ROOT_5;
	    
	    xp1yp2 = imRef(im, x + 1, y + 2) + ROOT_5;
	    xyp2 = imRef(im, x, y + 2) + ROOT_5;
	    xm1yp2p = imGetPixPtr(im, x - 1, y + 2);
	    xm1yp2 = imPtrRef(im, xm1yp2p) + ROOT_5;
	    
	    xlim = imGetPixPtr(im, xbase + 1, y);
	    for (; !imPtrEq(im, xyp, xlim); ) {
		xy = MIN(xy, xp1y);
		xy = MIN(xy, xp2yp1);
		xy = MIN(xy, xp1yp1);
		xy = MIN(xy, xyp1);
		xy = MIN(xy, xm1yp1);
		xy = MIN(xy, xm2yp1);
		xy = MIN(xy, xp1yp2);
		xy = MIN(xy, xm1yp2);

		/* Shift */
		imPtrRef(im, xyp) = xy;
		xp1y = xy + ONE;
		imPtrLeft(im, xyp);
		xy = imPtrRef(im, xyp);

		xp2yp1 = xp1yp1 + ROOT_5 - ROOT_2;
		xp1yp1 = xyp1 + ROOT_2 - ONE;
		xyp1 = xm1yp1 + ONE - ROOT_2;
		xm1yp1 = xm2yp1 + ROOT_2 - ROOT_5;
		imPtrLeft(im, xm2yp1p);
		xm2yp1 = imPtrRef(im, xm2yp1p) + ROOT_5;

		xp1yp2 = xyp2;
		xyp2 = xm1yp2;
		imPtrLeft(im, xm1yp2p);
		xm1yp2 = imPtrRef(im, xm1yp2p) + ROOT_5;

		/* Here is the original code:
		   minval = imRef(im, x, y);
		   minval = MIN(minval, imRef(im, x, y + 1) + ONE);
		   minval = MIN(minval, imRef(im, x + 1, y) + ONE);
		   minval = MIN(minval, imRef(im, x + 1, y + 1) + ROOT_2);
		   minval = MIN(minval, imRef(im, x + 1, y + 2) + ROOT_5);
		   minval = MIN(minval, imRef(im, x + 2, y + 1) + ROOT_5);
		   minval = MIN(minval, imRef(im, x - 1, y + 1) + ROOT_2);
		   minval = MIN(minval, imRef(im, x - 1, y + 2) + ROOT_5);
		   minval = MIN(minval, imRef(im, x - 2, y + 1) + ROOT_5);
		   imRef(im, x, y) = minval;
		   */
		}
	    }

	/* Second entry */
	if (width >= 4) {
	    x = xbase + 1;
	    minval = imRef(im, x, y);
	    minval = MIN(minval, imRef(im, x, y + 1) + ONE);
	    minval = MIN(minval, imRef(im, x + 1, y) + ONE);
	    minval = MIN(minval, imRef(im, x + 1, y + 1) + ROOT_2);
	    minval = MIN(minval, imRef(im, x + 1, y + 2) + ROOT_5);
	    minval = MIN(minval, imRef(im, x + 2, y + 1) + ROOT_5);
	    minval = MIN(minval, imRef(im, x - 1, y + 1) + ROOT_2);
	    minval = MIN(minval, imRef(im, x - 1, y + 2) + ROOT_5);
	    imRef(im, x, y) = minval;
	    }

	/* First entry */
	if (width >= 3) {
	    x = xbase;
	    minval = imRef(im, x, y);
	    minval = MIN(minval, imRef(im, x, y + 1) + ONE);
	    minval = MIN(minval, imRef(im, x + 1, y) + ONE);
	    minval = MIN(minval, imRef(im, x + 1, y + 1) + ROOT_2);
	    minval = MIN(minval, imRef(im, x + 1, y + 2) + ROOT_5);
	    minval = MIN(minval, imRef(im, x + 2, y + 1) + ROOT_5);
	    imRef(im, x, y) = minval;
	    }
	}

    /* And that's done it. What a mess... */
    }

/*
 * This routine actually performs the increasing-x distance transform.
 * It is passed a LongImage which contains the distance transform of
 * some image. It generates the increasing-x distance transform of
 * the binary image which consists of all points in the original image
 * which are under thresh.
 */
ShortImage
dtrans_plusx(LongImage dt, long thresh)
{
    int xbase, ybase;
    unsigned width, height;
    short minval;
    int x, y;
    ShortImage im;
    ShortPixPtr xdtp;
    LongPixPtr dtp;

    assert(dt != (LongImage)NULL);

    xbase = imGetXBase(dt);
    ybase = imGetYBase(dt);
    width = imGetWidth(dt);
    height = imGetHeight(dt);

    im = (ShortImage)imNewOffset(IMAGE_SHORT, width, height, xbase, ybase);
    if (im == (ShortImage)NULL) {
	return((ShortImage)NULL);
	}

    /* Now propagate things to the left */
    for (y = ybase; y < ybase + (int)height; y++) {
	/* All the entries - scan backwards */
	minval = INFINITY_SHORT;
	xdtp = imGetPixPtr(im, xbase + (int)width - 1, y);
	dtp = imGetPixPtr(dt, xbase + (int)width - 1, y);
	for (x = xbase + (int)width - 1; x >= xbase; x--) {
	    if (imPtrRef(dt, dtp) <= thresh) {
		imPtrRef(im, xdtp) = minval = 0;
		}
	    else {
		imPtrRef(im, xdtp) = ++minval;
		}
	    imPtrLeft(im, xdtp);
	    imPtrLeft(dt, dtp);
	    }
	}

    return(im);
    }

/*
 * This routine performs the increasing-y distance transform.
 * It is passed a LongImage which contains the distance transform of
 * some image. It generates the increasing-y distance transform of
 * the binary image which consists of all points in the original image
 * which are under thresh.
 */
ShortImage
dtrans_plusy(LongImage dt, long thresh)
{
    int xbase, ybase;
    unsigned width, height;
    short minval;
    int x, y;
    ShortImage im;

    assert(dt != (LongImage)NULL);

    xbase = imGetXBase(dt);
    ybase = imGetYBase(dt);
    width = imGetWidth(dt);
    height = imGetHeight(dt);

    im = (ShortImage)imNewOffset(IMAGE_SHORT, width, height, xbase, ybase);
    if (im == (ShortImage)NULL) {
	return((ShortImage)NULL);
	}

    /* Now propagate things upward */
    for (x = xbase; x < xbase + (int)width; x++) {
	/* All the entries - scan backwards */
	minval = INFINITY_SHORT;
	for (y = ybase + (int)height - 1; y >= ybase; y--) {
	    if (imRef(dt, x, y) <= thresh) {
		imRef(im, x, y) = minval = 0;
		}
	    else {
		imRef(im, x, y) = ++minval;
		}
	    }
	}

    return(im);
    }

/*
 * Transform a list of points to a LongImage, 0 where there is a point,
 * INFINITY where there isn't.
 */
static LongImage
pts_to_lim(unsigned width, unsigned height,
	   int xbase, int ybase,
	   unsigned npts, point *points)
{
    int x, y;
    int i;
    LongImage im;

    assert(points != (point *)NULL);

    im = (LongImage)imNewOffset(IMAGE_LONG, width, height, xbase, ybase);
    if (im == (LongImage)NULL) {
	return((LongImage)NULL);
	}

    /* Initialise the distance transform */
    for (y = ybase; y < ybase + (int)height; y++) {
	for (x = xbase; x < xbase + (int)width; x++) {
	    imRef(im, x, y) = INFINITY;
	    }
	}
    for (i = 0; i < npts; i++) {
	if (((points[i].x >= xbase) && (points[i].x < xbase + (int)width)) &&
	    ((points[i].y >= ybase) && (points[i].y < ybase + (int)height))) {
	    imRef(im, points[i].x, points[i].y) = 0;
	    }
	}

    return(im);
    }

/*
 * Transform a BinaryImage to a LongImage, 0 where the BinaryImage is
 * black, INFINITY where it is white.
 */
static LongImage
bim_to_lim(BinaryImage bim)
{
    int xbase, ybase;
    unsigned width, height;
    int x, y;
    LongImage im;

    assert(bim != (BinaryImage)NULL);
    xbase = imGetXBase(bim);
    ybase = imGetYBase(bim);
    width = imGetWidth(bim);
    height = imGetHeight(bim);

    im = (LongImage)imNewOffset(IMAGE_LONG, width, height, xbase, ybase);
    if (im == (LongImage)NULL) {
	return((LongImage)NULL);
	}

    /* Initialise the distance transform */
    for (y = ybase; y < ybase + (int)height; y++) {
	for (x = xbase; x < xbase + (int)width; x++) {
	    if (imRef(bim, x, y) == 0) {
		/* 0 means white - no point */
		imRef(im, x, y) = INFINITY;
		}
	    else {
		/* It's black - a point */
		imRef(im, x, y) = 0;
		}
	    }
	}

    return(im);
    }