ghough.c

snake算法

/****************************************************************************
 File Name  : ghough.c
 Purpose    : to implement generalized Hough transform
 Release    : Version 1.0
 Date	    : Aug 31,1995

GSNAKE API is jointly developed by the Information Technology Institute (ITI), Singapore, and the School of Applied Science, Nanyang Technological
University (NTU), Singapore. 

These software programs are available to the user without any license or royalty fees. Permission is hereby granted to use, copy, modify, and distribute this software and its documentation for any purpose. ITI and NTU gives no warranty, express, implied, or statuary for the software and/or documentation provided, including, without limitation, waranty of merchantibility and warranty of fitness for a particular purpose. The software provided hereunder is on an "as is"  basis, and ITI and NTU has no obligation to provide maintenance, support, updates, enhancements, or modifications.

GSNAKE API is available for any UNIX compatible system. User feedback, bugs, or software and manual suggestions should be sent via electronic mail to one of the following, who may or may not act on them as he/she desires :

		[e-mail]
		[e-mail]
***************************************************************************/
		 
#include "gsnake.h"
		  
  /*
   * we separate linear transformations into 3 components :
   *      T = T(r) T(t) T(d)
   *      where
   *              T(R) : scale change, indexed by 2 elements
   *                     R | rxy+  0   |   centered at 1 | 1 0 |
   *                       | 0    rxy- |                 | 0 1 |
   *              T(t) : rotation matrix indexed by *theta* centered at *0*
   *              T(d) : dilation matrix, indexed by 2 elements
   *                      | 1   dx |   centered at | 1 0 |
   *                      | dy  1  |               | 0 1 |
   * The Q values below control the *resolution* of the transformation, while
   * the N values control the # of cells, or allowable range. The total planes
   * needed are NR * Nrxy * Nt * Ndx * Ndy.
   */

GHOUGH::GHOUGH(int _Qx, int _Qy, 
               int _NR, double _QR, 
               int _Nrxy, double _Qrxy,
               int _Nt,  double _Qt,
               int _Ndx, double _Qdx,
               int _Ndy, double _Qdy,
               double _R, double _THETA)
{
	QR = _QR ;		/* resolution for scale change */
	Qrxy = _Qrxy ;		/* resolution for diagonal stretching */
	
	Qt  = _Qt;		/* resolution for rotation (in radian) */
	Qdx  = _Qdx;		/* resolutions for off diagonal stretching */
        Qdy  = _Qdy;
	Qx = _Qx;
	Qy = _Qy;

	NR = _NR;		/* # of cells in scale change axis */
	Nrxy = _Nrxy;
	Nt = _Nt;          	/* # of cells in the rotation axis */
	Ndx = _Ndx;		/* # of cells in each of the T(d) axis */
	Ndy = _Ndy;

	R = _R ; 		/* scaling factor for referenced snake */
	THETA = _THETA ; 	/* initial angle for referenced line snake */

	plane_max = -1; 	/* store index for blanking */
	row_max = -1;
	col_max = -1 ;
}		

GHOUGH::reset()
{
        register short i;

        for (i = 0; i<NumPlane;i++) {
                delete Planes[i];
                delete Template[i];
 		_iMemFree(Angle[i]);
        }

	_iMemFree( Template );
	_iMemFree( Planes );
	_iMemFree( Angle );
	
	plane_max = row_max = col_max = -1 ;
	NumPlane = 0 ;

}

GHOUGH::~GHOUGH()
{
	reset() ;
}

CONTOUR *GHOUGH::localize(CONTOUR *reference, IMAGE *img)
{
	/* if memory allocation fail, do nothing to reference */
	if( init(img->getRow(), img->getCol()) != NOERROR )
		return reference->duplicate() ;
	
	transform(reference, 0); /* get various instances of ref. contour */
	
	correlate(img) ;	/* perform correlation */

	CONTOUR *result = getMax() ;
	
	reset() ;
	
	return result ;
}

CONTOUR *GHOUGH::localize(CONTOUR *reference, EDGE *edgeMap)
{
	
	/* if memory allocation fail, do nothing to reference */
	if( init(edgeMap->getRow(), edgeMap->getCol()) != NOERROR )
		return reference->duplicate() ;
	
	transform(reference); /* get various instances of ref. contour */


	correlate(edgeMap) ;	/* perform correlation */
	
	CONTOUR *result = getMax() ;

	
	reset() ;
	
	return result ;
}

CONTOUR **GHOUGH::localize(CONTOUR *reference, IMAGE *img, int numFind)
{
	CONTOUR **resultlist ;
	
	/* if memory allocation fail, return NULL */
	if( init(img->getRow(), img->getCol()) != NOERROR ||
	    !(resultlist = (CONTOUR **) _iMemAlloc(numFind*sizeof(CONTOUR *))) )
		return NULL ;
	
	transform(reference, 0); /* get various instances of ref. contour */
	
	correlate(img) ;	/* perform correlation */
	
	register short i ;
	for(i=0; i<numFind; i++)
		resultlist[i] = getMax() ;
	
	reset() ;
	
	return resultlist ;
}

CONTOUR **GHOUGH::localize(CONTOUR *reference, EDGE *edgeMap, int numFind)
{
	CONTOUR **resultlist ;
	
	/* if memory allocation fail, return NULL */
	if( init(edgeMap->getRow(), edgeMap->getCol()) != NOERROR ||
	    !(resultlist = (CONTOUR **) _iMemAlloc(numFind*sizeof(CONTOUR *))) )
		return NULL ;
	
	transform(reference); /* get various instances of ref. contour */
	
	correlate(edgeMap) ;	/* perform correlation */
	
	register short i ;
	for(i=0; i<numFind; i++)
		resultlist[i] = getMax() ;
	
	reset() ;
	
	return resultlist ;
}



int GHOUGH::init(int row, int col)
{
	register short i;
	short m,n;
	
	NumPlane = NR*Nrxy*Nt*Ndx*Ndy;
	
	
	Planes = (IMAGE **) _iMemAlloc(NumPlane * sizeof(IMAGE *));
	Template = (CONTOUR **) _iMemAlloc (NumPlane * sizeof(CONTOUR *));
	Angle = (double **) _iMemAlloc(NumPlane * sizeof(double*));

			
	if ( !Planes || !Template || !Angle ) {
                fprintf(stderr,
                        "<GHOUGH::init> : cannot allocate GHT pointers\n");
                return MEMORYERROR;
        }
	
	m = short (row / Qy) +1;
	n = short (col / Qx) +1;
	for (i=0; i<NumPlane; i++) {
		
		Planes[i] = new IMAGE;
		if( !Planes[i] || Planes[i]->init(m, n) != NOERROR ) {
                        fprintf(stderr,
                           "<GHOUGH::init> : cannot allocate GHT planes\n");
                        return MEMORYERROR;
                }
	}

	return NOERROR;
}

void GHOUGH::transform(CONTOUR *reference, char comAngleFlag)
{
        short iR, irxy, it, idx, idy;
        int index ;
    
        reference->computeCG();   	/* compute center of gravity */
        reference->contourCentered();     /* convert to U coordinate */

        /* generate various transformations */
         index = 0 ;
         for(iR=0; iR<NR; iR++)
          for(irxy=0; irxy<Nrxy; irxy++)
            for(it=0; it<Nt; it++)
              for(idx=0; idx<Ndx; idx++)
                for(idy=0; idy<Ndy; idy++, index++) {

        		double rx  = R * (1 + (int)(irxy-(int)(Nrxy/2))*Qrxy) *
        			     ( 1 + (int)(iR-(int)(NR/2))*QR) ;
        		double ry  = R * (1 + (int)((int)(Nrxy/2)-irxy)*Qrxy) *
        			     ( 1 + (int)(iR-(int)(NR/2))*QR) ;
        		double t   = THETA + (int)(it-(int)(Nt/2))*Qt;
        		double tdx = 0 + (int)(idx-(int)(Ndx/2))*Qdx;
        		double tdy = 0 + (int)(idy-(int)(Ndy/2))*Qdy;

			Template[index] = reference->duplicate() ;
			Template[index]->affineTransform( 
					rx, ry, t, 0, 0, tdx, tdy );
			
	           }

	if( !comAngleFlag ) {
		memset(Angle, 0, index*sizeof(double *) );
		return ;
	}
	
	/*
 	 * get desired angle (gradient of snaxel)
 	 */
	 int i, j ;
	 for(i=0; i<index; i++) {
	
	      Angle[i] = (double *) _iMemAlloc(sizeof(double)*
	 				reference->getNumSnaxel() ) ;
  
              SNAXEL *active ; 
              double *angle = Angle[i] ;
              for(active = Template[i]->getHead(), j=0 ; active ;
              		active = active->getNext(), j++) {


                /* compute tangent direction */
                angle[j] = active->getNormalAng(
			Template[i]->getMode(), 
			Template[i]->getHead(), 
			Template[i]->getTail() );

               }
         }
}

/*
 * perform GHT correlation
 */
void GHOUGH::correlate(EDGE *edgeMap)
{
	register short m, n, i, j;
	int sM, sN, eM, eN;

	/* portion of image to where object can be centered */
	sM = 0 ;
	eM = Planes[0]->getRow() - 1 ;
	sN = 0;
	eN = Planes[0]->getCol() - 1 ;

	for (m=0; m<edgeMap->getRow(); m++)
	  for (n=0; n<edgeMap->getCol(); n++) {

		double inc; /* increase plane when match */
		double mag = (double) edgeMap->getMag(m, n);

		if (mag < 0.1) /* skip insignificant values */
			continue;

		double angle = (double) edgeMap->getAng(m, n);

		/* increment for all templates */
                for (i=0; i<NumPlane; i++) {
		   int Cx, Cy;
		   SNAXEL *sxptr = Template[i]->getHead();

		   for ( j=0; sxptr; j++, sxptr = sxptr->getNext() ) {

			Cx = n - ROUNDOFF( sxptr->getCol() );
			Cy = m - ROUNDOFF( sxptr->getRow() );

			if (Qx > 1)
			   Cx = (int) ((double)Cx/(double)Qx) ;
			if (Qy > 1)
			   Cy = (int) ((double)Cy/(double)Qy);

			/* calculate projection angle and inc.  planes */
			if (Cx == CAP(sN, Cx, eN) && Cy == CAP(sM, Cy, eM)) {
			   inc = mag*cos(Angle[i][j] - angle);
			   Planes[i]->put(Cy,Cx,Planes[i]->get(Cy,Cx) + inc);
			}
		   }
		}
	} 
}

/*
 * perform GHT correlation
 */
void GHOUGH::correlate(IMAGE *img)
{
	register short m, n, i, j;
	int sM, sN, eM, eN;
	

	/* portion of image to where object can be centered */
	sM = 0 ;
	eM = Planes[0]->getRow() - 1 ;
	sN = 0;
	eN = Planes[0]->getCol() - 1 ;

	for (m=0; m<img->getRow(); m++)
	  for (n=0; n<img->getCol(); n++) {

		double intensity = (double) img->get(m, n);

		if (intensity < 0.1) /* skip insignificant values */
			continue;

		/* increment for all templates */
                for (i=0; i<NumPlane; i++) {
		   int Cx, Cy;
		   SNAXEL *sxptr = Template[i]->getHead();

		   for ( j=0; sxptr; j++, sxptr = sxptr->getNext() ) {

			Cx = n - ROUNDOFF(sxptr->getCol());
			Cy = m - ROUNDOFF(sxptr->getRow());

			if (Qx > 1)
			   Cx = (int) ((double)Cx/(double)Qx) ;
			if (Qy > 1)
			   Cy = (int) ((double)Cy/(double)Qy);

			/* calculate projection angle and inc.  planes */
			if (Cx == CAP(sN, Cx, eN) && Cy == CAP(sM, Cy, eM)) 
			   Planes[i]->put(Cy,Cx,Planes[i]->get(Cy,Cx) + 
			   		intensity);
	
		   }
		}
	} 
}

/*
 * find peak
 */
CONTOUR *GHOUGH::getMax( void )
{
	register short i, j, row, col;
	double maxval = 0;
	double sign;
	CONTOUR *contour;
	

	if (plane_max >= 0 && row_max >= 0 && col_max >= 0) {
		for(i=-1; i<=1; i++)
		  for(j=-1; j<=1; j++)
			Planes[plane_max]->put(
			     CAP(0, i+row_max, Planes[plane_max]->getRow()-1),
			     CAP(0, j+col_max, Planes[plane_max]->getCol()-1), 
			     0);
	 }

	/* search for maximun */
	for (i=0; i<NumPlane; i++) {
	  for (row=0; row<Planes[i]->getRow(); row++) {
	    for (col=0; col<Planes[i]->getCol(); col++) {
		double data = Planes[i]->get(row,col);
		if (fabs(data) > maxval) {
			plane_max = i;
			row_max = row;
			col_max = col;
			maxval = fabs(data);
			sign = SIGN(data);
		}
	    }
	  }
	}

	contour = Template[plane_max]->duplicate();
	contour->putCgRow( (double)(row_max * Qy) + Qy/2);
	contour->putCgCol( (double)(col_max * Qx) + Qx/2);
	contour->putDirection( (short)sign );
	contour->imageCentered();
	return contour;
}