scale-h.c

Hausdorff Distance for Image Recognition

/*
 *
 * $Header: /usr/u/wjr/vision/scale-h/RCS/scale-h.c,v 1.49 1994/04/07 21:09:55 wjr Exp $
 *
 * Copyright (c) 1992, 1993 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/vision/scale-h/RCS/scale-h.c,v 1.49 1994/04/07 21:09:55 wjr Exp $";

#include <stdio.h>
#include "misc.h"
#include <math.h>
#include "image.h"
#include "scale-h.h"
#include "list.h"
#include "chunk.h"
#include "dbug.h"

#include "prof.h"

#ifdef	MP
#include <signal.h>
#include <sys/wait.h>
#endif

/* About how fine to slice each dimension at the lowest resolution */
#define	SLICE		4 /* 63 */

/* What the resolution refinement between steps is */
#define	REFINE		2

/* What the minimum step size at top level we're willing to accept is */
#define	MINTOPSTEP	4

/* The usual FP epsilon */
#define	EPSILON		1e-7

/* How much of the model to scan before deciding on early decision */
#define	EARLY_FRAC	0.2

/* How wide a beam in the beam search, initially */
#define	BEAMWIDTH	16

/* and how much bigger it gets on a backtrack */
#define	BEAMSPREAD	2

/* Instrumentation */
#ifdef	INSTRUMENT
static long nprobes = 0;
static long nabort = 0;
static long npick = 0;
static long nearly = 0;
static long nwrong = 0;
static long nfull = 0;

static long nprobes1 = 0;
static long nabort1 = 0;
static long npick1 = 0;
static long nearly1 = 0;
static long nwrong1 = 0;
static long nfull1 = 0;

static int level;

#ifdef INST2
#define MAXPT 10000

static boolean decision[MAXPT];

static int ntt[MAXPT];
static int ntf[MAXPT];
static int nft[MAXPT];
static int nff[MAXPT];

static int ntt1[MAXPT];
static int ntf1[MAXPT];
static int nft1[MAXPT];
static int nff1[MAXPT];

static int nbelow[MAXPT];
static int nbelow1[MAXPT];

#endif
#endif

/* Search and verification functions */
static void findSTrans(simage_info *image, smodel_info *model,
		       int forwstyle, int revstyle);
static void checkSTransMtoI(simage_info *image, smodel_info *model);
static void checkSTransItoM(simage_info *image, smodel_info *model,
			    int revstyle);
static boolean checkOneSTransMtoI(simage_info *image, smodel_info *model,
				  stransval *t);
static boolean evalOneSTransMtoI(simage_info *image, smodel_info *model,
				 int x, int y, int xs, int ys,
				 long *dist, double *frac, double *frac2);
static boolean checkOneSTransItoM(simage_info *image, smodel_info *model,
				  long reverse_thresh, double reverse_frac,
				  int revstyle,
				  int x, int y, int xs, int ys,
				  boolean *result);
extern void findSTransMtoI(simage_info *image, smodel_info *model);
static void findSTransOne(simage_info *image, smodel_info *model,
			  int revstyle);
static void findSTransBestDist(simage_info *image, smodel_info *model,
			       int revstyle);
static void findSTransBestFrac(simage_info *image, smodel_info *model,
			       int revstyle);

extern boolean findSTransMtoISampled(simage_info *image, smodel_info *model,
				     int lowX, int highX, int stepX,
				     int lowY, int highY, int stepY,
				     int lowScaleX, int highScaleX,
				     int stepScaleX,
				     int lowScaleY, int highScaleY,
				     int stepScaleY,
				     PriQ pq,
				     long thresh, boolean needdata);

static void BFS(simage_info *image, smodel_info *model, int revstyle,
		boolean (*cmpfunc)(void *pqh1, void *pqh2),
		boolean (*matchfunc)(simage_info *image, smodel_info *model,
				     int revstyle, spqHook *hook));
static boolean bfsrecur(simage_info *image, smodel_info *model,
			PriQ lowpq,
			boolean (*cmpfunc)(void *pqh1, void *pqh2),
			boolean (*matchfunc)(simage_info *image,
					     smodel_info *model,
					     int revstyle, spqHook *hook),
			int revstyle,
			int lowX, int highX, int stepX,
			int lowY, int highY, int stepY,
			int lowScaleX, int highScaleX, int stepScaleX,
			int lowScaleY, int highScaleY, int stepScaleY,
			int curbeam);

static void hillclimb(simage_info *image, smodel_info *model, int revstyle);

extern boolean probeRegions(simage_info *image, smodel_info *model,
			    PriQ lowpq, PriQ highpq,
			    int lowX, int highX, int stepX, int oldStepX,
			    int lowY, int highY, int stepY, int oldStepY,
			    int lowScaleX, int highScaleX, int stepScaleX,
			    int oldStepScaleX,
			    int lowScaleY, int highScaleY, int stepScaleY,
			    int oldStepScaleY,
			    long thresh, boolean needdata, int nexpand);

#ifdef	MP
static void mp_probeRegions(void);
#endif

/* Support functions */
static boolean compute_scaled_pts(smodel_info *model);

static void shuffle_model(smodel_info *model);

static void calcLowCells(simage_info *image, smodel_info *model,
			 int *lowX, int *highX, int *stepX,
			 int *lowY, int *highY, int *stepY,
			 int *lowScaleX, int *highScaleX, int *stepScaleX,
			 int *lowScaleY, int *highScaleY, int *stepScaleY);
static void updateCellSteps(smodel_info *model,
			    int stepX, int stepY,
			    int stepScaleX, int stepScaleY,
			    int *newStepX, int *newStepY,
			    int *newStepScaleX, int *newStepScaleY);
static void calcBoxes(smodel_info *model, int stepX, int stepY,
		      int stepScaleX, int stepScaleY,
		      unsigned *box_width, unsigned *box_height);
static void convertPtToBox(int x, int y, int scaleX, int scaleY,
			   int lowX, int stepX, int oldStepX,
			   int lowY, int stepY, int oldStepY,
			   int lowScaleX, int stepScaleX, int oldStepScaleX,
			   int lowScaleY, int stepScaleY, int oldStepScaleY,
			   int *botX, int *topX, int *botY, int *topY,
			   int *botSX, int *topSX,
			   int *botSY, int *topSY);

static boolean inrange(simage_info *image, smodel_info *model,
		       int x, int y, int xs, int ys);

static int roundlog(int x, int n, int m);

/* Comparison functions for priq sorting */
static boolean falseCompare(void *a, void *b);
static boolean findDistCompare(void *a, void *b);
static boolean findOneCompare(void *a, void *b);
static boolean findFracCompare(void *a, void *b);

/* Match functions for search */
static boolean findOneMatch(simage_info *image, smodel_info *model,
			    int revstyle, spqHook *hook);
static boolean findDistMatch(simage_info *image, smodel_info *model,
			     int revstyle, spqHook *hook);
static boolean findFracMatch(simage_info *image, smodel_info *model,
			     int revstyle, spqHook *hook);

#ifdef	MP
/* Multiprocess stuff */
static boolean mp_startup(simage_info *image, smodel_info *model, int forwstyle);
static void mp_shutdown(void);

static void killall();

extern void (*signal(int sig, void (*func)()))();

static pid_t pid[NPROC];	/* Who is out there? */
static pid_t mypid;		/* who am I? */
static int mynum;		/* and what's my number? */

static void *shmem;		/* Shared mem stuff */
static int shmid;

static int semid;

/* Some semaphore numbers */
#define	MAIN_Q_LOCK	0		/* Access to main queue */
#define	SECOND_Q_LOCK	1		/* Access to second queue */
#define	MAIN_Q_SLEEP	2		/* Waiting for main queue non-empty */
#define	SECOND_Q_SLEEP	3		/* Waiting for second queue non-full */
#define	MAINPROC_WAIT	4		/* Main process waiting for everyone else */

#define	NSEMS		5

/* Some stuff to deal with the job queues */

/* Max length of each one */
#define	Q_SIZE		5000

/* How many things you should grab from a queue at once, max */
#define	JOBCHUNK	100

/* The sorts of things that can be stuck in a queue */
#define	Q_WORK		1
#define	Q_DIE		2

/* and an actual queue entry */
typedef struct {
    int whattodo;
    spqHook pqh;
    } qent;

/*
 * Queue pointers and the queues. The main process is the only writer on
 * the main queue and only reader on the secondary. Note that the queue
 * pointers are shared. Any access to this structure must be protected
 * by a lock.
 * The box sizes etc are basically the values passed to probeRegions. They
 * stay constant until the level of the search shifts. This will never
 * happen unless all the subsidiary processes are sitting around waiting
 * for jobs, and so the subsidiary processes can just look at them when they
 * pull a job chunk off the queue.
 *
 * Note that I'm passing pointers (specifically model and image) through
 * shared memory. This is fine as long as they point to the right thing
 * in every process's address space. This is true since they start out
 * as clones of each other, at some time *after* model and image are
 * generated.
 */
typedef struct {
    int main_first;
    int main_last;
    int main_waiters;
    int second_first;
    int second_last;
    int second_waiters;
    int jobsize;
    boolean (*cmpfunc)(void *pqh1, void *pqh2);
    unsigned box_width;
    unsigned box_height;
    int lowX, highX, stepX, oldStepX;
    int lowY, highY, stepY, oldStepY;
    int lowScaleX, highScaleX, stepScaleX, oldStepScaleX;
    int lowScaleY, highScaleY, stepScaleY, oldStepScaleY;
    boolean needdata;
    long thresh;
    smodel_info *model;
    simage_info *image;
    long model_thresh;
    double model_frac;
    qent main[Q_SIZE];
    qent second[Q_SIZE];
    } jobs;

#define	SHM_SIZE	(sizeof(jobs))

static volatile jobs *queues;

/* Some operations */

/* Lock and unlock the queues */
static struct sembuf lock_main_q_op = { MAIN_Q_LOCK, -1, 0 };
static struct sembuf unlock_main_q_op = { MAIN_Q_LOCK, 1, 0 };
static struct sembuf lock_second_q_op = { SECOND_Q_LOCK, -1, 0 };
static struct sembuf unlock_second_q_op = { SECOND_Q_LOCK, 1, 0 };

static union semun dummyarg = { 0 };

/*
 * This is using the semaphore-zero stuff of SunOS: you can wait
 * for a semaphore to become zero. These semaphores are normally 1,
 * and if the queue is empty (full), the processes start sleeping
 * on them (waiting for them to become zero). Once they are ready
 * again (non-empty, non-full resp), the main process decrements the
 * semaphore, everybody wakes up, and it increments it again.
 * Note that they have to unlock the main lock for the queue before they
 * go to sleep since otherwise nobody could get in to wake them up. I'm
 * assuming (since man semop isn't clear) that:
 * Given an array of sembufs, semop will process the operations in them
 * 1) in order and
 * 2) atomically.
 * Thus, we unlock the queue and go to sleep as one action - there is no
 * timing hole between unlocking and sleeping (i.e. nobody can sneak in
 * and remove the need for us to go to sleep, so we'll never get woken
 * up since that event won't happen again).
 *
 * Hum. These assumptions appear to be incorrect. Of course.
 * What it seems to do is more like wait until it can perform them all
 * at once, then do them. I guess this makes sense, since then it really
 * is an atomic action. I guess I get to do it as two operations and live
 * with the (tiny) timing hole until I decide to figure out how to fix it.
 * FIX.
 *
 * There's an added complication (of course) - the main process needs to
 * know when everyone has finished processing the jobs they took, so it
 * needs to wait for them. Of course, they might block partway through because
 * the secondary queue filled up, and so they'll need the main process to
 * come and clean it out for them. Thus, whenever a process goes to sleep
 * on either queue (waiting for main non-empty or secondary non-full), it
 * increments the MAINPROC_WAIT semaphore. The main process can therefore
 * safely sleep on this semaphore, and whenever anyone else increments
 * it, it will wake up and do what is necessary: either proceed because all
 * the work has been done, or clean out the secondary queue and go back to
 * sleep, or go back to sleep because someone is still working.
 */
static struct sembuf presleep_main_q_op[2] = { { MAINPROC_WAIT, 1, 0 },
					       { MAIN_Q_LOCK, 1, 0 } };
static struct sembuf sleep_main_q_op = { MAIN_Q_SLEEP, -1, 0 };
static struct sembuf wakeup_main_q_op = { MAIN_Q_SLEEP, 0, 0 };
static struct sembuf presleep_second_q_op[2] = { { MAINPROC_WAIT, 1, 0 },
						 { SECOND_Q_LOCK, 1, 0 } };
static struct sembuf sleep_second_q_op = { SECOND_Q_SLEEP, -1, 0 };
static struct sembuf wakeup_second_q_op = { SECOND_Q_SLEEP, 0, 0 };

static struct sembuf mainproc_do_wait_op = { MAINPROC_WAIT, -1, 0 };

/*
 * Initial values:
 * The lock semaphores should be 1, so one process can grab each lock.
 * The sleep semaphores should be 0, so processes initially sleep on them.
 * The mainproc semaphore should be 0, so the main process sleeps on it
 * until someone else sleeps on one of the others.
 */
static ushort semarr[NSEMS] = { 1, 1, 0, 0, 0 };

/* Some macros to deal with the queues */
#define	LOCK_MAIN_Q	semop(semid, &lock_main_q_op, 1)
#define	UNLOCK_MAIN_Q	semop(semid, &unlock_main_q_op, 1)
#define	LOCK_SECOND_Q	semop(semid, &lock_second_q_op, 1)
#define	UNLOCK_SECOND_Q	semop(semid, &unlock_second_q_op, 1)

#define	PRESLEEP_MAIN_Q	semop(semid, &presleep_main_q_op[0], 2)
#define	SLEEP_MAIN_Q	semop(semid, &sleep_main_q_op, 1)
#define	WAKEUP_MAIN_Q	semop(semid, &wakeup_main_q_op, 1)
#define	PRESLEEP_SECOND_Q	semop(semid, &presleep_second_q_op[0], 2)
#define	SLEEP_SECOND_Q	semop(semid, &sleep_second_q_op, 1)
#define	WAKEUP_SECOND_Q	semop(semid, &wakeup_second_q_op, 1)

#define	MAINPROC_DO_WAIT	semop(semid, &mainproc_do_wait_op, 1)

#endif

/*
 * forwstyle controls how the search is performed.
 * If it is FORWARD_ALL, we search for all transformations where
 * the specified model_frac of the model pixels are under model_thresh.
 * In this case, we can also pass in a list of transformations in model->trans;
 * these will be checked out and those that meet the criteria will be returned;
 * others will be deleted.
 * If it is FORWARD_ONE, we search for *any* transformation which
 * meets this requirement, and report it. It must also satisfy the reverse
 * criteria (unless revstyle is REVERSE_NONE).
 * If it is FORWARD_BESTDIST, we search for all transformations where
 * the model_frac'th distance is minimised: find a transformation (which must
 * satisfy the reverse criteria) and use its distance until a better one
 * is found etc. model->model_thresh is set to this best distance.
 * If it is FORWARD_BESTFRAC, we do a similar search, except that we want
 * to find transformations where the number of points less than model_thresh
 * are maximised. model->model_frac is set to this best fraction.
 *
 * FORWARD_ONE can be modified by FORWARD_HILLCLIMB: find a match, and
 * then hillclimb from there to a local minimum. Also, with FORWARD_ALL and
 * trans != NULLLIST (i.e. passing in a list of transformations to be
 * evaluated), if FORWARD_HILLCLIMB is set it will hillclimb out of all
 * transformations which meet the criteria.
 */
void
findSTransAllModels(simage_info *image, smodel_info *models, unsigned nmodels,
		    int forwstyle, int revstyle)
{
    int i;
    int max_leftborder;
    int max_topborder;
    int max_rightborder;
    int max_bottomborder;

    assert(models != (smodel_info *)NULL);
    assert(image != (simage_info *)NULL);

    if (nmodels == 0) {
	return;		/* What's to do? */
	}

    max_leftborder = models[0].leftborder;
    max_topborder = models[0].topborder;
    max_rightborder = models[0].rightborder;
    max_bottomborder = models[0].bottomborder;
    for (i = 1; i < nmodels; i++) {
	max_leftborder = MAX(max_leftborder, models[i].leftborder);
	max_topborder = MAX(max_topborder, models[i].topborder);
	max_rightborder = MAX(max_rightborder, models[i].rightborder);
	max_bottomborder = MAX(max_bottomborder, models[i].bottomborder);
	}

    if ((image->dtrans == (LongImage)NULL) ||
	(image->leftborder < max_leftborder) ||
	(image->topborder < max_topborder) ||
	(image->rightborder < max_rightborder) ||
	(image->bottomborder < max_bottomborder)) {

	/* Were they being very silly? */
	if ((max_leftborder + max_rightborder + (int)image->xsize <= 0) ||
	    (max_topborder + max_bottomborder + (int)image->ysize <= 0)) {
	    /*
	     * They asked for borders which leave no part of the image.