gsnake.c

snake算法

/****************************************************************************
 File Name  : gsnake.c
 Purpose    : provides g-snake manipulation routines, includes energy
	      minimization, minimax regularization and initialization
	      by 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"
#include "xwindow.h"

GSNAKE::GSNAKE(EEXTTYPE etype, double glambda)
{
	EextType = etype;
	global_lambda = glambda;
	EdgeMap = NULL;
	GaussImg = NULL;
}

GSNAKE::~GSNAKE(void)
{
	reset();
}

void GSNAKE::reset(void)
{
	 EdgeMap = NULL;
	 GaussImg = NULL;
}

GSNAKE::genPyramid(short level, short verbose)
{
	if(EextType == _INTENSITY) 
	  return PYRAMID::generate(level, verbose, 1, 1, 1, 1, 1, _INTENSITY);
	else
	  return PYRAMID::generate(level, verbose) ;
}

/* localize gsnake at desired place. 
   19 May - For original version, localise will crash if called when pyramid 
            is not built.
            Implemented to check for this, will use root image if no pyramid.
            If no root image, will exit 
   
 */
int GSNAKE::localize(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)	
{
	CONTOUR *contour;
	GHOUGH GHT_OBJ(_Qx, _Qy, _NR, _QR, _Nrxy, _Qrxy,
               	_Nt, _Qt, _Ndx, _Qdx, _Ndy, _Qdy, _R, _THETA);

        short level = PYRAMID::getLevel() - 1;
       
        if (level<0) {
        	printf("\nAuto generating pyramid to level 1 to localize.\n");
        	PYRAMID::generate(1,1);
        	level = 0;
        }
        
  	expand( LEVEL(-level) );
  	
  	GaussImg = getGauss(level) ;
	EdgeMap = getEdge(level) ;

 
  	/* generalize hough transform */
	switch(EextType) {
	    case _INTENSITY :
  		contour = GHT_OBJ.localize(this, GaussImg) ;
		break ;
	    case _EDGE :
		contour = GHT_OBJ.localize(this, EdgeMap) ;
		break ;
	    case _EDGEMAG :
		contour = GHT_OBJ.localize(this, EdgeMap->getMagImg());
		break ;
	}
        
	contour->duplicate( (CONTOUR *) this); 

	if ( _Qx > 1 || _Qy > 1 ) 
		fineLocalize( _Qx, _Qy); 

	GaussImg = getGauss(0) ;
        EdgeMap = getEdge(0) ;


	/* Expand contour to original size */
	CONTOUR::expand(LEVEL(level));
	
  
	delete contour ;

	return NOERROR;
}

/*
 *  perform fine correlation around center of gravity
    
    Edit : 22 June tmp now duplicates a contour instead of the whole gsnake
 
 */
 
void GSNAKE::fineLocalize(  int _qx,	/* translation resolution */
			    int _qy,
			    int _cg_col,/* centre of gravity */
			    int _cg_row )
{
	register short i, j ;
	SNAXEL *sxptr;
	int max_i, max_j ;
	double max_match, match ;

	/* at least correlate in 3x3 window */
	_qx = MAX(_qx, 3) ;
	_qy = MAX(_qy, 3) ;

	 /* get gaussian pyramid */
        GaussImg = PYRAMID::getGauss( 0 );
        EdgeMap  = PYRAMID::getEdge( 0 );

	/* move centre of gravity if neccessary */
	if ( _cg_col || _cg_row ) {
		translate( ROUNDOFF( _cg_col - cg_col ) , 
			   ROUNDOFF( _cg_row - cg_row ) );
	}
						    
	CONTOUR *tmp = CONTOUR::duplicate();
	
	max_i = max_j = 0 ;
	max_match = -getNumSnaxel() ;
	for (i = - (_qy/2) ; i <= _qy/2 ; i++)
	   for(j = - (_qx/2) ; j <= _qx/2 ; j++) {

		tmp->duplicate(this, 1);
		translate( j, i );
		match = 0;

		/* compute external energy */
		for (sxptr=getHead(); sxptr; sxptr=sxptr->getNext())
			match += EExternal( sxptr ) ;

		match = -match;	

		/* select maximun energy */
		if ( match > max_match ) {
			max_match = match;
			max_i = i;
			max_j = j;
		} 
	}

	tmp->duplicate(this, 1);
	translate( max_j, max_i );
	
	delete tmp;  
}

/*
 * Use a gaussian pyramid to perform coarse to fine minimization
 */
int  GSNAKE::minimize( 
		int segmentSpacing,   /* coarse search spacing */
                int numSearchSegment, /* region to move around */
                int snaxelSpacing,   /* snaxel spacing */
                unsigned char blowup, /* blow up ratio */
                int verbose,
                double theLambda,     /* lambda */
                int showImg,		/* show image ? */
		int img_Xoffset,	/* offset from image */
		int img_Yoffset)
{
        short curLevel;                 /* current gaussian pyramid level */
        int not_minimize = 1;		/* Initialise to not minimized */ 
	double old_lambda;
	
        /* get gaussian pyramid level */
        short level = PYRAMID::getLevel() - 1;

        if (level<0) {
        	printf("\nAuto generating pyramid to level 1 to minimize.\n");
        	PYRAMID::generate(1,1);
        	level = 0;
        }
        
        old_lambda = global_lambda;
        if( theLambda != _DEFINED_LAMBDA) 
             	global_lambda = theLambda;
        
	
	/* Reduce contour back to highest level */
	expand(LEVEL(-level)); 
       
        /* beginning from top-most level, do coarse to fine minimization */
        for (curLevel = level; curLevel >= 0; --curLevel) {

                not_minimize = mn_dymprog(curLevel,
                        segmentSpacing, numSearchSegment, blowup, verbose,
			showImg, img_Xoffset, img_Yoffset);

                 /* go down the levels if neccessary */
                if ( curLevel ) {

		        /* expand snake */
                        expand( LEVEL(1), snaxelSpacing );
                     
                        computeShape() ;

                        /* all verification step must use local minimax */
                        if(snaxelSpacing) global_lambda =  _LOCAL_MINMAX;
                }

        }
	global_lambda = old_lambda;

        return not_minimize;
}

/* 
 * marginalize gsnake to get probability 
 */
double GSNAKE::marginalize(int nhoodTangent, int nhoodNormal)
{
	LOCATIONS searchLocs ;
	nhoodTangent = 2*((int) nhoodTangent/2) + 1;
	nhoodNormal  = 2*((int) nhoodNormal/2)  + 1 ;
	
	/* must use lowest level */
        EdgeMap = getEdge(0) ;
        GaussImg = getGauss(0) ;
        
        if( searchLocs.init(numsnaxel, nhoodTangent*nhoodNormal) == MEMORYERROR)
        {
            fprintf(stderr, 
            	"ERROR : <GSNAKE::dymprog> memory allocation error\n");
	    return MEMORYERROR;
	}

	/* generate square window for marginalization */
	searchLocs.generateLoc(mode, head, tail, 
		nhoodNormal/2, nhoodTangent/2, 1, 0, 1) ;
	
	return dymprog(1, &searchLocs, 1) * 
			 exp(getNumSnaxel()/sig_nu_sqr) / getZ() ;		
}


/*
 * compute total energy of a gsnake object
 */
double GSNAKE::ESnake( short level, int verbose )    /* pyramid edge map */
{
        SNAXEL *pnow;           /* current snaxel */
        
        switch(EextType) {

            case _INTENSITY :
            	if( !(GaussImg = getGauss(level)))
			return PARAMERROR ;
            		break ;

            case _EDGE :
            case _EDGEMAG :
        	if ( !(EdgeMap = getEdge(level))  || 
		     !(GaussImg = getGauss(level)) )
			return PARAMERROR ;
			break;
       }
	
        /* compute average distant of snaxel coordinate */
        computeAvgLength() ;

        /* get initial value */
        Esnake = 0;     /* snake energy */
        Eext = 0;       /* average gradient power on snake */
    
        /* compute snake energy */

        for ( pnow = getHead(); pnow; pnow = pnow->getNext() ) {

                Esnake += ESnaxel( pnow,  _TRUE, verbose );
                Eext += pnow->getEext();                /* gradient power */
        }

        /* compute data strength information */
        Eext    /= (double) numsnaxel;  /* ave. gradient power */
  
        /* compute ave. snake energy */
        Esnake /= (double) numsnaxel;

	
        return Esnake;
}

/*
 * compute  total energy at current snaxel
 */
double GSNAKE::ESnaxel( SNAXEL *now,	/* current snaxel */
                        BOOLEAN store,
                        int verbose )
{
        double Eimage;  /* edge energy */
        double Emodel;  /* model energy */
        double Epoint;  /* total energy at current snaxel */

	double row = now->getRow();
	double col = now->getCol();

        /* check out of range */
        if (row != CAP(0, row, GaussImg->getRow()-1) ||
            col != CAP(0, col, GaussImg->getCol()-1) ) {

                if ( store == _TRUE ) {
                        now->putEext( 1 );
                        now->putEsnaxel( 1 );
                }
                return 1.0;
        }
        
        Eimage = EExternal(now);

        /* compute internal (model) energy */
        Emodel = EInternal( now );

        /* regularize */
        regularize( now, &Emodel, &Eimage, &Epoint);

        /* store energy term if neccessary */
        if (store == _TRUE) {
                now->putEint( Emodel );
                now->putEext( Eimage );
                now->putEsnaxel( Epoint );
        }

        if (verbose){
        	
           fprintf(stdout, "[%d, %d]\t%1.3f\t%1.3f\t%1.3f\t%1.3f\n",
                 ROUNDOFF(col), ROUNDOFF(row), Emodel, Eimage, Epoint ,