dtrans.c

Hausdorff Distance for Image Recognition

/*
 *
 * $Header: /usr/u/wjr/src/support/RCS/dtrans.c,v 1.11 1993/07/26 22:16:16 wjr Exp $
 *
 * Copyright (c) 1990, 1991, 1992, 1993 Cornell University.  All Rights
 * Reserved.
 *
 * Copyright (c) 1991, 1992 Xerox Corporation.  All Rights Reserved.
 *  
 * Use, reproduction, preparation of derivative works, and distribution
 * of this software is permitted.  Any copy of this software or of any
 * derivative work must include both the above copyright notices of
 * Cornell University and Xerox Corporation and this paragraph.  Any
 * distribution of this software or derivative works must comply with all
 * applicable United States export control laws.  This software is made
 * available AS IS, and XEROX CORPORATION DISCLAIMS ALL WARRANTIES,
 * EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE,
 * AND NOTWITHSTANDING ANY OTHER PROVISION CONTAINED HEREIN, ANY
 * LIABILITY FOR DAMAGES RESULTING FROM THE SOFTWARE OR ITS USE IS
 * EXPRESSLY DISCLAIMED, WHETHER ARISING IN CONTRACT, TORT (INCLUDING
 * NEGLIGENCE) OR STRICT LIABILITY, EVEN IF XEROX CORPORATION IS ADVISED
 * OF THE POSSIBILITY OF SUCH DAMAGES.
 */

static char rcsid[] = "@(#)$Header: /usr/u/wjr/src/support/RCS/dtrans.c,v 1.11 1993/07/26 22:16:16 wjr Exp $";

#include "misc.h"
#include "image.h"
#include "readpts.h"
#include "panic.h"

#include <prof.h>

/* Some values, scaled by DSCALE */
#define	ROOT_2	141
#define	ROOT_5	224
#define	ONE	100

/* A number larger than any possible distance transform value */
#define	INFINITY	(1 << 30)
#define	INFINITY_SHORT	(1 << 14)
/* N.B. It must be possible to increment INFINITY_SHORT by at least width */

static void do_dtrans(LongImage im);
static LongImage pts_to_lim(unsigned width, unsigned height,
			    int xbase, int ybase,
			    unsigned npts, point *points);
static LongImage bim_to_lim(BinaryImage bim);

/*
 * Perform a distance transform of a point set.
 *
 * This routine generates the distance transform of a point set. It uses the
 * 5x5 approximation to the L_2 norm, scaled by DSCALE.
 */
LongImage
dtrans_pts(unsigned width, unsigned height,
	   int xbase, int ybase,
	   unsigned npts, point *points)
{
    LongImage im;

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

    im = pts_to_lim(width, height, xbase, ybase, npts, points);
    if (im == (LongImage)NULL) {
	return((LongImage)NULL);
	}

    do_dtrans(im);

    return(im);
    }

/*
 * This distance computes the distance transform of the
 * image: the distance to the nearest black pixel.
 */
LongImage
dtrans_image(BinaryImage bim)
{
    LongImage im;

    assert(bim != (BinaryImage)NULL);

    im = bim_to_lim(bim);
    if (im == (LongImage)NULL) {
	return((LongImage)NULL);
	}

    do_dtrans(im);
    return(im);
    }

/*
 * This routine actually performs the required distance transform.
 * It is passed a LongImage which has been initialised to INFINITY
 * where there is no set pixel and 0 where there is. It modifies the
 * LongImage so that the value at each point is the distance to the nearest
 * set pixel (pixel where there was a 0 initially). The distance is an
 * approximation to the L_2 norm, scaled by 100.
 *
 * The transform is done using a mask which is
 *			ROOT_5		ROOT_5
 *		ROOT_5	ROOT_2	  ONE	ROOT_2	ROOT_5
 *			ONE	   0	ONE
 *		ROOT_5	ROOT_2	  ONE	ROOT_2	ROOT_5
 *			ROOT_5		ROOT_5
 *
 * It is done in 2 passes: first forward using the first half of the mask
 * (the parts preceding the centre in standard reading order) and second
 * backward using the second half of the mask.
 * At each point, it computes MIN(im(x+dx,y+dy)+mask(dx,dy)) for values in
 * the appropriate half of the mask.
 */
static void
do_dtrans(LongImage im)
{
    int xbase, ybase;
    unsigned width, height;
    long minval;
    int x, y;

    assert(im != (LongImage)NULL);

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

    /* Now, do the forward pass. Handcode boundary conditions for speed. */

    /* First line */
    if (height >= 1) {
	y = ybase;
	for (x = xbase + 1; x < xbase + (int)width; x++) {
	    imRef(im, x, y) = MIN(imRef(im, x, y),
				  imRef(im, x - 1, y) + ONE);
	    }
	}

    /* Second line */
    if (height >= 2) {
	y = ybase + 1;
	for (x = xbase; x < xbase + (int)width; x++) {
	    minval = imRef(im, x, y);
	    minval = MIN(minval, imRef(im, x, y - 1) + ONE);

	    /* Make sure we don't look over the edge */
	    if (x >= xbase + 1) {
		minval = MIN(minval, imRef(im, x - 1, y) + ONE);
		minval = MIN(minval, imRef(im, x - 1, y - 1) + ROOT_2);
		}
	    if (x >= xbase + 2) {
		minval = MIN(minval, imRef(im, x - 2, y - 1) + ROOT_5);
		}
	    if (x < xbase + (int)width - 1) {
		minval = MIN(minval, imRef(im, x + 1, y - 1) + ROOT_2);
		}
	    if (x < xbase + (int)width - 2) {
		minval = MIN(minval, imRef(im, x + 2, y - 1) + ROOT_5);
		}
	    imRef(im, x, y) = minval;
	    }
	}

    /* Other lines */
    for (y = ybase + 2; y < ybase + (int)height; y++) {
	/* First entry in the line */
	x = xbase;
	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 entry */
	if (width >= 2) {
	    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);
	    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;
	    }
	    
	/* All the rest of the entries */
	/*
	 * Don't try to understand this code unless you really have to.
	 * It loads a bunch of values out of the dtrans, then does its own
	 * shifting between them. It also biases each value depending on its
	 * position relative to the current centre. These biases are adjusted
	 * when the values are shifted. The rationale for doing it this way
	 * is because on many current processors, adds cost the same as
	 * moves.
	 * One exception: xym2 (representing the value at (x, y - 2) is
	 * not biased by 2 (or ONE+ONE), but by ROOT_5, because it isn't
	 * actually used in the computation, and the two values adjacent to
	 * it are, and are both biased by ROOT_5.
	 */
	if (width >= 5) {
	    long xm1y, xy;
	    long xm2ym1, xm1ym1, xym1, xp1ym1, xp2ym1;
	    long xm1ym2, xym2, xp1ym2;
	    LongPixPtr xyp, xp2ym1p, xp1ym2p;
	    LongPixPtr xlim;
	    
	    /* Start it out */
	    x = xbase + 2;

	    xm1y = imRef(im, x - 1, y) + ONE;
	    xyp = imGetPixPtr(im, x, y);
	    xy = imPtrRef(im, xyp);
	    
	    xm2ym1 = imRef(im, x - 2, y - 1) + ROOT_5;
	    xm1ym1 = imRef(im, x - 1, y - 1) + ROOT_2;
	    xym1 = imRef(im, x, y - 1) + ONE;
	    xp1ym1 = imRef(im, x + 1, y - 1) + ROOT_2;
	    xp2ym1p = imGetPixPtr(im, x + 2, y - 1);
	    xp2ym1 = imPtrRef(im, xp2ym1p) + ROOT_5;
	    
	    xm1ym2 = imRef(im, x - 1, y - 2) + ROOT_5;
	    xym2 = imRef(im, x, y - 2) + ROOT_5;
	    xp1ym2p = imGetPixPtr(im, x + 1, y - 2);
	    xp1ym2 = imPtrRef(im, xp1ym2p) + ROOT_5;
	    
	    xlim = imGetPixPtr(im, xbase + (int)width - 2, y);
	    for (; !imPtrEq(im, xyp, xlim); ) {
		xy = MIN(xy, xm1y);
		xy = MIN(xy, xm2ym1);
		xy = MIN(xy, xm1ym1);
		xy = MIN(xy, xym1);
		xy = MIN(xy, xp1ym1);
		xy = MIN(xy, xp2ym1);
		xy = MIN(xy, xm1ym2);
		xy = MIN(xy, xp1ym2);

		imPtrRef(im, xyp) = xy;

		/* Shift */
		xm1y = xy + ONE;
		imPtrRight(im, xyp);
		xy = imPtrRef(im, xyp);

		xm2ym1 = xm1ym1 + ROOT_5 - ROOT_2;
		xm1ym1 = xym1 + ROOT_2 - ONE;
		xym1 = xp1ym1 + ONE - ROOT_2;
		xp1ym1 = xp2ym1 + ROOT_2 - ROOT_5;
		imPtrRight(im, xp2ym1p);
		xp2ym1 = imPtrRef(im, xp2ym1p) + ROOT_5;

		xm1ym2 = xym2;
		xym2 = xp1ym2;
		imPtrRight(im, xp1ym2p);
		xp1ym2 = imPtrRef(im, xp1ym2p) + 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 last entry */
	if (width >= 4) {
	    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);
	    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;
	    }

	/* Last entry */
	if (width >= 3) {
	    x = xbase + (int)width - 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);
	    imRef(im, x, y) = minval;
	    }
	}

    /* Now the reverse pass. The same mess, but with signs changed. Yay. */

    /* Last line */
    if (height >= 1) {
	y = ybase + (int)height - 1;
	for (x = xbase + (int)width - 2; x >= xbase; x--) {
	    imRef(im, x, y) = MIN(imRef(im, x, y),
				  imRef(im, x + 1, y) + ONE);
	    }
	}

    /* Second last line */
    if (height >= 2) {
	y = ybase + (int)height - 2;