immax.c

Hausdorff Distance for Image Recognition

/*
 *
 * $Header: /usr/u/wjr/src/support/RCS/immax.c,v 1.6 1994/09/01 23:22:42 wjr Exp $
 *
 * Copyright (c) 1993, 1994 Cornell University.  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 notice of
 * Cornell University 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 CORNELL UNIVERSITY 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
 * CORNELL UNIVERSITY IS ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
 */

static char rcsid[] = "@(#)$Header: /usr/u/wjr/src/support/RCS/immax.c,v 1.6 1994/09/01 23:22:42 wjr Exp $";

/*
 * This file contains routines which "convolve" a rectangular max or min
 * operator with a LongImage. For example, the image
 *	3   1   4   1   5
 *	9   2   6   5   3
 *	5   8   9   7   9
 *	3   2   3   8   4
 *	6   2   6   4   3
 * if convolved with a 2x2 max operator gives:
 *	9   6   6   5   5
 *	9   8   9   9   9
 *	8   9   9   9   9
 *	6   6   8   8   4
 *	6   6   6   4   3
 * and with a 2x2 min operator,
 *	1   1   1   1   3
 *	2   2   5   3   3
 *	2   2   3   4   4
 *	2   2   3   3   3
 *	2   2   4   3   3
 * i.e. the value at (x,y) is MAX(A(x,y),A(x+1,y),A(x,y+1),A(x+1,y+1)) (or
 * MIN(...)).
 * In general, a rectangle of size (w, h) is centred at the current pixel
 * (x,y) and the pixels from (x,y)...(x+w-1,y+h-1) are examined.
 * The max (or min) value is stored at (x,y). Pixels beyond the boundary are
 * ignored.
 *
 */

#include "image.h"
#include "misc.h"
#include "prof.h"

static long opmax(long a, long b);
static long opmin(long a, long b);

static LongImage imOp(LongImage im, unsigned width, unsigned height,
		      long (*op)(long a, long b));

static boolean ppinit(unsigned oplen, unsigned rowlen);
static void ppfinal(void);
static void ppmap(long *row, unsigned rowlen, unsigned oplen,
		  long (*op)(long a, long b));

static int floorlog2(unsigned n);

static int nrows = 0;
static long **rows = (long **)NULL;

/*
 * Compute the max as described above, and return the resulting
 * LongImage, which will be the same size as the original. Return
 * (LongImage)NULL in case of failure.
 */
LongImage
imMax(LongImage im, unsigned width, unsigned height)
{
    assert(im != (LongImage)NULL);

    /* width or height of zero doesn't make enough sense to be worth it */
    assert((width != 0) && (height != 0));

    return(imOp(im, width, height, opmax));
    }

/*
 * Compute the min as described above, and return the resulting
 * LongImage, which will be the same size as the original. Return
 * (LongImage)NULL in case of failure.
 */
LongImage
imMin(LongImage im, unsigned width, unsigned height)
{
    assert(im != (LongImage)NULL);

    /* width or height of zero doesn't make enough sense to be worth it */
    assert((width != 0) && (height != 0));

    return(imOp(im, width, height, opmin));
    }

/*
 * A general routine: map op over rectangular regions in the image.
 * op must be associative and commutative.
 * Return (LongImage)NULL in case of failure.
 *
 * It must be possible to separate the rectangle: map op over the rows
 * then over the columns. We read out each row into an array, and pass
 * it off to the parallel-prefix mapping routing, then do the same with
 * each column of the array this generated.
 * If the box is 1 wide (or 1 high) we skip the mapping phase for that
 * direction.
 */
static LongImage
imOp(LongImage im, unsigned width, unsigned height,
     long (*op)(long a, long b))
{
    int x;
    int y;
    LongImage newim;
    LongPixPtr p;
    long *q;
    long *row;

    assert(im != (LongImage)NULL);
    assert(op != (long (*)(long a, long b))NULL);

    /* width or height of zero doesn't make enough sense to be worth it */
    assert((width != 0) && (height != 0));

    /* Copy the original array */
    newim = (LongImage)imDup(im);
    if (newim == (LongImage)NULL) {
	return((LongImage)NULL);
	}

    row = (long *)malloc(sizeof(long) * MAX(imGetWidth(im), imGetHeight(im)));
    if (row == (long *)NULL) {
	imFree(newim);
	return((LongImage)NULL);
	}

    if (!ppinit(MAX(width, height),
		MAX(imGetWidth(im), imGetHeight(im)))) {
	free((void *)row);
	imFree(newim);
	return((LongImage)NULL);
	}

    if (width > 1) {
	/* Map across the rows */
	for (y = imGetYBase(newim); y <= imGetYMax(newim); y++) {
	    /* Read out the row */
	    for (q = row, p = imGetPixPtr(newim, imGetXBase(newim), y),
		 x = imGetXBase(newim);
		 x <= imGetXMax(newim);
		 q++, imPtrRight(newim, p), x++) {
		*q = imPtrRef(newim, p);
		}
	    
	    /* Process the row */
	    ppmap(row, imGetWidth(im), width, op);
	    
	    /* and write it back */
	    for (q = row, p = imGetPixPtr(newim, imGetXBase(newim), y),
		 x = imGetXBase(newim);
		 x <= imGetXMax(newim);
		 q++, imPtrRight(newim, p), x++) {
		imPtrRef(newim, p) = *q;
		}
	    }
	}

    if (height > 1) {
	/* Map down the columns */
	for (x = imGetXBase(newim); x <= imGetXMax(newim); x++) {
	    /* Read out the column */
	    for (q = row, p = imGetPixPtr(newim, x, imGetYBase(newim)),
		 y = imGetYBase(newim);
		 y <= imGetYMax(newim);
		 q++, imPtrDown(newim, p), y++) {
		*q = imPtrRef(newim, p);
		}
	    
	    /* Process the column */
	    ppmap(row, imGetHeight(im), height, op);
	    
	    /* and write it back */
	    for (q = row, p = imGetPixPtr(newim, x, imGetYBase(newim)),
		 y = imGetYBase(newim);
		 y <= imGetYMax(newim);
		 q++, imPtrDown(newim, p), y++) {
		imPtrRef(newim, p) = *q;
		}
	    }
	}

    free((void *)row);
    ppfinal();
    return(newim);
    }

/*
 * For efficiency's sake, we don't want to allocate and free all the
 * temporary rows that ppmap uses every time that it's called. Therefore,
 * you have to call ppinit() with the maximum oplen and rowlen
 * that ppmap() will be called on. It returns TRUE if successful,
 * FALSE if failed.
 */
static boolean
ppinit(unsigned oplen, unsigned rowlen)
{
    int i, j;

    assert(oplen > 0);

    /* Make sure we're currently not inited */
    assert(nrows == 0);

    /* Find floor(log_2(oplen)) */
    nrows = floorlog2(oplen) + 1;
    rows = (long **)malloc(sizeof(long *) * (unsigned)nrows);
    if (rows == (long **)NULL) {
	nrows = 0;
	return(FALSE);
	}
    for (i = 0; i < nrows; i++) {
	rows[i] = (long *)malloc(sizeof(long) * rowlen);
	if (rows[i] == (long *)NULL) {
	    for (j = 0; j < i; j++) {
		free((void *)rows[j]);
		}
	    free((void *)rows);
	    nrows = 0;
	    return(FALSE);
	    }
	}
    
    return(TRUE);
    }

static void
ppfinal(void)
{
    int i;

    /* We must have been inited */
    assert(nrows > 0);
    assert(rows != (long **)NULL);

    for (i = 0; i < nrows; i++) {
	assert(rows[i] != (long *)NULL);
	free((void *)rows[i]);
	}
    free((void *)rows);
    nrows = 0;
    }
    
/*
 * Perform a parallel prefix mapping of op over row. Use the temp storage
 * in rows[0..nrows - 1] for the mapping. Only go as far as is required for
 * the given oplen.
 */
static void
ppmap(long *row, unsigned rowlen, unsigned oplen,
      long (*op)(long a, long b))
{
    int i, j;
    unsigned val;
    int offset;
    long *p, *q, *r;
    int edge;

    assert(row != (long *)NULL);
    assert(nrows > 0);
    assert(rows != (long **)NULL);
    assert(rowlen > 0);
    assert(oplen > 0);
    assert(nrows >= floorlog2(oplen) + 1);
    assert(op != (long (*)(long a, long b))NULL);

    /* Copy the first row */
    assert(rows[0] != (long *)NULL);
    for (p = rows[0], q = row, j = 0; j < rowlen; p++, q++, j++) {
	*p = *q;
	}

    /*
     * We want:
     * rows[i][j] = op(row[j], row[j+1], ..., row[j+2^i-1])
     * (ignoring any values >= rowlen).
     * This is now true for rows[0], since it's a copy of row.
     * We make it true by computing
     * rows[i][j] = op(rows[i-1][j], rows[i-1][j + 2^(i - 1)])
     * or just
     * rows[i][j] = op(rows[i-1][j]) if j + 2^(i - 1) >= rowlen.
     */
    /*
     * Hack in special support for opmin and opmax, to reduce funcall overhead.
     */
    edge = rowlen - (1 << (i - 1));
    if (op == opmin) {
	for (i = 1; i <= floorlog2(oplen); i++) {
	    assert(rows[i] != (long *)NULL);

	    edge = rowlen - (1 << (i - 1));
	    /* Do the bit before the edge effect sets in */
	    for (p = rows[i], q = rows[i-1],
		 r = rows[i-1] + (1 << (i - 1)), j = 0;
		 j < edge;
		 p++, q++, r++, j++) {
		*p = MIN(*q, *r);
		}

	    /* and now do the edge effect bit */
	    for (; j < rowlen; p++, q++, j++) {
		*p = *q;
		}
	    }
	}
    else if (op == opmax) {
	for (i = 1; i <= floorlog2(oplen); i++) {
	    assert(rows[i] != (long *)NULL);

	    edge = rowlen - (1 << (i - 1));
	    /* Do the bit before the edge effect sets in */
	    for (p = rows[i], q = rows[i-1],
		 r = rows[i-1] + (1 << (i - 1)), j = 0;
		 j < edge;
		 p++, q++, r++, j++) {
		*p = MAX(*q, *r);
		}

	    /* and now do the edge effect bit */
	    for (; j < rowlen; p++, q++, j++) {
		*p = *q;
		}
	    }
	}
    else {
	for (i = 1; i <= floorlog2(oplen); i++) {
	    assert(rows[i] != (long *)NULL);

	    edge = rowlen - (1 << (i - 1));
	    /* Do the bit before the edge effect sets in */
	    for (p = rows[i], q = rows[i-1],
		 r = rows[i-1] + (1 << (i - 1)), j = 0;
		 j < edge;
		 p++, q++, r++, j++) {
		*p = op(*q, *r);
		}

	    /* and now do the edge effect bit */
	    for (; j < rowlen; p++, q++, j++) {
		*p = *q;
		}
	    }
	}
	
    /*
     * OK - now look at the bit decomposition of oplen: if oplen is
     * 2^i1 + 2^i2 + ... + 2^in, to calculate the final value of row[j],
     * look at
     * op(rows[i1][j] , rows[i2][j + 2^i1], ...,
     *	  rows[i3][j + 2^i1 + ... + 2^(in-1)]).
     */
    offset = 0;
    for (i = 0, val = oplen; val; i++, val >>= 1) {
	assert(i < nrows);

	if (val & 1) {
	    /* oplen had bit i set */
	    if ((val << i) == oplen) {
		/* This is the first set bit we've found - copy */
		for (p = row, q = rows[i], j = 0; j < rowlen; p++, q++, j++) {
		    *p = *q;
		    }
		offset = (1 << i);
		}
	    else {
		/* op into the current values in row */
		if (op == opmin) {
		    for (p = row, q = rows[i] + offset, j = offset;
			 j < rowlen;
			 p++, q++, j++) {
			*p = MIN(*p, *q);
			}
		    }
		else if (op == opmax) {
		    for (p = row, q = rows[i] + offset, j = offset;
			 j < rowlen;
			 p++, q++, j++) {
			*p = MAX(*p, *q);
			}
		    }
		else {
		    for (p = row, q = rows[i] + offset, j = offset;
			 j < rowlen;
			 p++, q++, j++) {
			*p = op(*p, *q);
			}
		    }
		offset += (1 << i);
		}
	    }
	}

    /* That did it... */
    }

static long
opmax(long a, long b)
{
    return(MAX(a, b));
    }

static long
opmin(long a, long b)
{
    return(MIN(a, b));
    }

/*
 * Compute floor(log_2(n)).
 */
static int
floorlog2(unsigned n)
{
    int i;

    assert(n > 0);

    for (i = 0; n; i++, n >>= 1) {
	}
    /*
     * We know that at every loop iteration but the exiting one,
     * 2^(floor(log_2(n))) <= 2^i * n <= n.
     * When the loop exits, n has been shifted right by i bits, and is zero;
     * n shifted right by i-1 bits must therefore have been 1, and so
     * i-1 = floor(log_2(n)).
     */
    return(i - 1);
    }