contour.c

snake算法

/****************************************************************************
 File Name  : contour.c
 Purpose    : provides basic contour manipulation routines and serves as 
	      internal energy input for g-snake
 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 "random.h"
#include "xwindow.h"

CONTOUR::CONTOUR( void ) 
{
	head = NULL;
	tail = NULL;
	numsnaxel = 0;
	sig_nu_sqr = 1;
	Z = 1 ;
	direction = 0 ;
}

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

/*
 * initialize snaxels ( link list )
 */
int CONTOUR::initSnaxel( short numpts )
{
        register short i;
	SNAXEL *sxptr;

	/* delete all snaxel if exists */
	for ( sxptr = head; sxptr; sxptr = sxptr->getNext() )
		delete sxptr;

        /* number of snaxel */
        numsnaxel = numpts;

        if ( !(sxptr =  new SNAXEL) )
                return MEMORYERROR;

        sxptr->putPrev( NULL );

        /* set the forward and backward pointers */
        for(i=1, head=sxptr; i<numpts; i++) {

		SNAXEL *next;

		if( !(next = new SNAXEL) )
			return MEMORYERROR;

                sxptr->putNext( next );
                sxptr->getNext()->putPrev( sxptr );
                sxptr = sxptr->getNext();
                sxptr->putNext( NULL );
        }

        tail = sxptr; /* tail of snake */

        return NOERROR;
}

/*
 * _iMemFree all snaxel 
 */
void CONTOUR::reset( void )
{
	SNAXEL *sxptr, *next;

	for ( sxptr = head; sxptr; sxptr = next) {
		next = sxptr->getNext() ;
		delete sxptr ;
	}
		
	head = tail = NULL;
}

/*
 * init a round (closed) snake
 */
int CONTOUR::init( short _cg_row, short _cg_col, double radius, short numpts )
{
        double step, angle;
        SNAXEL *sxptr;
        int Imgrow = 2*_cg_row;
        int Imgcol = 2*_cg_col;
	int i=0;

        /* set gsnake parameters */
        mode = _CLOSED;
        numsnaxel = numpts;
	cg_row = _cg_row;
	cg_col = _cg_col;

        /* generate snaxels */
        if ( initSnaxel(numpts) == MEMORYERROR ) {
           fprintf(stderr," <CONTOUR::init> : memory allocation error\n");
           return MEMORYERROR;
        }

        step = 2*M_PI / (double) numpts;

        /* defaulte values for shape matrix */
        for (angle=0,sxptr=head; sxptr; angle+=step,sxptr=sxptr->getNext()) {

                sxptr->putCol( CAP(0,_cg_col+radius*cos(angle),Imgcol) );
                sxptr->putRow( CAP(0,_cg_row+radius*sin(angle),Imgrow) );
                sxptr->putAlpha( 0.5 / cos(2*M_PI/(double) numpts) );
                sxptr->putBeta( 0.5 / cos(2*M_PI/(double) numpts) );
                sxptr->putLambda( 0.99 );
        }

        return NOERROR;
}

/*
 * init a line (opened) snake
 */
int CONTOUR::init( short sx, short sy,     /* starting points */
                   short ex, short ey,     /* end points */
                   short numpts )           /* # of snaxels on this snake */
{
        double len_x, len_y;
        double dx, dy;
        register short i;
        SNAXEL *sxptr;

        /* calculate number of point needed  */
        len_x = (double) (ex - sx) ;
        len_y = (double) (ey - sy) ;

        /* set gsnake parameters */
        mode = _OPENED;
        numsnaxel = numpts;

        /* generate snaxels */
        if ( initSnaxel(numpts) == MEMORYERROR ) {
           fprintf(stderr," <CONTOUR::init> : memory allocation error\n");
           return MEMORYERROR;
        }

        dx = len_x/(double) (numsnaxel -1);
        dy = len_y/(double) (numsnaxel -1);

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

                sxptr->putCol( dx * i + sx );
                sxptr->putRow( dy * i + sy );
                sxptr->putAlpha( 0.5 );
		sxptr->putBeta( 0.5 );
                sxptr->putLambda( 0.99 );
        }

        /* adjust head and tail of snake */
        adjust();

        return NOERROR;
}

/*
 * init snake by drag / click mouse
 */
int CONTOUR::init( IMAGE *ximg,        /* reference image */
                unsigned char blowup,   /* blowup ratio */
                INITMODE  theINIT,      /* click/drag mode */
                SNAKEMODE _mode,
                int spacing )       	/* desired spacing of snake */
{
        short numpts;
        short *xdata, *ydata;

        /* allocate space for coordinates */
        ydata = (short *)_iMemAlloc( MAX_SNAKEPTS * sizeof(short) );
        xdata = (short *)_iMemAlloc( MAX_SNAKEPTS * sizeof(short) );

        if ( !ydata || !xdata ) {
           fprintf(stderr, "<CONTOUR::init> : memory allocation error\n");
                return MEMORYERROR;
        }

        /* show reference image */
        ximg->show(blowup);
        ximg->show(blowup);

        /* read mouse input */
        ( theINIT == _CLICKMOUSE ) ?
                xwin_InitSnakeClick(ydata, xdata, &numpts) :
                xwin_InitSnakeDrag(ydata, xdata, &numpts, 2*spacing);

        /* set gsnake parameters */
        mode = _mode;

        /* init snaxels */
        if ( init( ydata, xdata, numpts, blowup) == MEMORYERROR ) {
           fprintf(stderr, "<CONTOUR::init> : memory allocation error\n");
           return MEMORYERROR;
        }

        /* compute average power of snaxels */
        computeAvgLength();

        /* learn shape matrix */
        computeShape();

        /* _iMemFree memory */
        _iMemFree ( ydata );
        _iMemFree ( xdata );

        return NOERROR;
}

/*
 * init gsnake by fill in snaxel value
 */
int CONTOUR::init( short *ydata, 	/* row data */
                   short *xdata,   	/* col data */
                   short numpts,   	/* number of snaxel */
                   unsigned char ratio)    /* blow up ratio */
{
        short register i;
        SNAXEL *sxptr;

	if ( numpts != numsnaxel )
		reset();

        /* generate snaxels */
	if ( !getHead() )
          if ( initSnaxel(numpts) == MEMORYERROR ) {
             fprintf(stderr," <CONTOUR::init> : memory allocation error\n");
             return MEMORYERROR;
          }

        for ( i=0, sxptr = getHead(); sxptr; i++, sxptr = sxptr->getNext() ) {
                sxptr->putRow( *(ydata + i) / ratio );
                sxptr->putCol( *(xdata + i) / ratio );
                /* sxptr->putLambda( 0.0 ); */
        }

        return NOERROR;
}

/*
 * init gsnake by fill in snaxel value
 */
int CONTOUR::init( double *ydata,		/* row data */
                   double *xdata,  		/* col data */
                   short numpts,		/* number of snaxel */
                   unsigned char ratio )	/* blow up ratio */
{
        short register i;
        SNAXEL *sxptr;

	if ( numpts != numsnaxel )
		reset();

        /* generate snaxels */
	if ( !getHead() )
          if ( initSnaxel(numpts) == MEMORYERROR ) {
             fprintf(stderr," <CONTOUR::init> : memory allocation error\n");
             return MEMORYERROR;
          }

        for ( i=0, sxptr = getHead(); sxptr; i++, sxptr = sxptr->getNext() ) {
                sxptr->putRow( *(ydata + i) / ratio );
                sxptr->putCol( *(xdata + i) / ratio );
                /* sxptr->putLambda( 0.0 ); */
        }

        return NOERROR;
}

/*
 * read from file
 */
int CONTOUR::read(char *filename)
{
        FILE *fid;      /* file identifier */
        SNAXEL *sxptr;  /* snaxel */
        int magic;      /* magic number */
        int _numsnaxel;
	float _sig_nu_sqr;
        SNAKEMODE _mode;
        float _col, _row, _alpha, _beta, _lambda, _z;
        short i;

        /* open file */
        if ( !( fid = fopen( filename, "r" ) ) ) {
           fprintf(stderr, "<CONTOUR :: read> : %s cannot be read\n", filename);
           return FILEIOERROR;
        }

        /* read magic number */
        fscanf(fid, "%x\n", &magic);

        if ( magic != CONTOURMAGICNUMBER ) {
            fprintf(stderr, "<CONTOUR :: read> : %s wrong format \n", filename);
            return FILEIOERROR;
        }

        /* read gsnake mode */
        fscanf(fid,"%d\n", &_mode);

        /* read number of snaxels */
        fscanf(fid,"%d\n", &_numsnaxel);

	/* read white noise var */
        fscanf(fid,"%f\n", &_sig_nu_sqr);

        /* init snaxel */
        mode =  _mode ;
        numsnaxel =  _numsnaxel ;
	sig_nu_sqr = (double) _sig_nu_sqr;
        if ( initSnaxel(_numsnaxel) == MEMORYERROR ) {
           fprintf(stderr," <CONTOUR::read> : memory allocation error\n");
           return MEMORYERROR;
        }

        /* read snaxel parameters */
        for ( sxptr = getHead(), i = 0; sxptr && i < _numsnaxel;
                sxptr = sxptr->getNext(), i++ ) {

                fscanf(fid,"%f ", &_col);
                fscanf(fid,"%f ", &_row);
                fscanf(fid,"%f ", &_alpha);
                fscanf(fid,"%f ", &_beta);
                fscanf(fid,"%f \n", &_lambda);

                sxptr->putCol( (double)_col );
                sxptr->putRow( (double)_row );
                sxptr->putAlpha( (double)_alpha );
                sxptr->putBeta( (double)_beta );
                sxptr->putLambda( (double)_lambda );
        }
        if ( fscanf(fid,"%f\n", &_z) > 0 )
		Z = _z ;
	else	Z = 1.0 ;

        /* close file */
        fclose(fid);
        
        computeAvgLength() ;
        
        return NOERROR;
}

/*
 * write to file
 */
int CONTOUR::write(char *filename)
{
        FILE *fid;      /* file identifier */
        SNAXEL *sxptr;  /* snaxel */

        /* open file */
        if ( !( fid = fopen( filename, "w" ) ) ) {
                fprintf(stderr, "<CONTOUR :: write> : %s cannot be write\n",
                        filename);
                return FILEIOERROR;
        }

        /* write magic number */
        fprintf(fid,"%x\n", CONTOURMAGICNUMBER);

        /* write to file */
        fprintf(fid,"%d\n", getMode());
        fprintf(fid,"%d\n", getNumSnaxel());
	fprintf(fid,"%e\n", getSigNuSqr());

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

                fprintf(fid,"%f ",sxptr->getCol());
                fprintf(fid,"%f ",sxptr->getRow());
                fprintf(fid,"%f ",sxptr->getAlpha());
                fprintf(fid,"%f ",sxptr->getBeta());
                fprintf(fid,"%f\n",sxptr->getLambda());
        }

	fprintf(fid,"%e\n", getZ());

        /* close file */
        fclose(fid);
        return NOERROR;
}

/*
 * energy : model constraint
 */
double CONTOUR::EInternal( SNAXEL *now )	/* current snaxel */
{
        double l_inf;           /* infinity norm */
        double meanrow, meancol;

        if ( fabs(now->getAlpha() ) < VERY_SMALL &&
             fabs(now->getBeta() )  < VERY_SMALL )
                return 1.0;

        /* compute shape coeff. (alpha & beta) with infinity norm */
        now->meanPosition(           /* mean position */
                mode, cg_row, cg_col,
                head, tail, &meanrow, &meancol);
        
        l_inf = MAX( fabs( meanrow - now->getRow() ),
                     fabs( meancol - now->getCol() ) );
        return( SQR( l_inf / avglen ) );
}

/*
 * learn A matrix coefficients for a snake
 * input is in U coordinate
 */
void CONTOUR::computeShape()
{
        SNAXEL *now, *prev, *next;
        double Xi, Yi, Xalpha, Yalpha, Xbeta, Ybeta, denom;

  
        computeAvgLength();

        contourCentered();         /* convert to U coordinate */
    
 
        for ( now=getHead(); now; now=now->getNext() ) {

                prev = now->getPrev( mode, tail );
                next = now->getNext( mode, head );

                /* compute shape coefficient (alpha & beta) */
                Xi = now->getCol() / avglen;
                Yi = now->getRow() / avglen;
                Xalpha = prev->getCol() / avglen;
                Yalpha = prev->getRow() / avglen;
                Xbeta  = next->getCol() / avglen;
                Ybeta  = next->getRow() / avglen;

                denom = Xalpha * Ybeta - Xbeta * Yalpha;

                if ( fabs(denom) < VERY_SMALL ) {
                        now->putAlpha( 0.5 ); 
			now->putBeta(  0.5 );
		}
                else {
                        now->putAlpha( (Xi*Ybeta - Yi*Xbeta) / denom );
                        now->putBeta(  (Yi*Xalpha - Xi*Yalpha) /denom );
                }

        }

    
        imageCentered();         /*  convert to V coordinate */

        if ( mode == _OPENED ) adjust();
}

/*
 * using shape matrix and the last 2 contour points, regenerate the entire
 * contour
 */
CONTOUR::regenerate()
{
	MATRIX Ar ;	/* reduced shape matrix */
	register int i, j ;
	int status ;
	
	if(numsnaxel < 3) return PARAMERROR ;	/* cannot regenerate */
	
	if( (status = Ar.init(numsnaxel-2, numsnaxel-2) ) != NOERROR ) 
	    return status ;
	
	SNAXEL *sptr ;
	for(i=0, sptr = getHead()->getNext(); 
			i<numsnaxel-2; i++, sptr=sptr->getNext() ) {
	
		Ar.put(i, i, -sptr->getAlpha()) ;
		if(i < numsnaxel-3) Ar.put(i, i+1, 1.0) ;
		if(i < numsnaxel-4) Ar.put(i, i+2, -sptr->getBeta() ) ;
	}
	
	/* regeneration algorithm */
	MATRIX *ArInv = Ar.inverse() ;
	SNAXEL *Un = getTail() ;
	SNAXEL *Un1 = Un->getPrev();
	SNAXEL Vn, Vn1;
	
	computeCG();
	contourCentered() ;
	
	Vn.putRow(Un1->getRow() - Un1->getBeta()*Un->getRow());
	Vn.putCol(Un1->getCol() - Un1->getBeta()*Un->getCol());
	Vn1.putRow(-Un1->getPrev()->getBeta()*Un1->getRow());
	Vn1.putCol(-Un1->getPrev()->getBeta()*Un1->getCol());

	for(i=0, sptr = getHead(); i<numsnaxel-2; i++, sptr=sptr->getNext() ) {
	
		sptr->putRow(-ArInv->get(i, numsnaxel-4)*Vn1.getRow() -
			     ArInv->get(i, numsnaxel-3)*Vn.getRow() ) ;
		sptr->putCol(-ArInv->get(i, numsnaxel-4)*Vn1.getCol() -
			     ArInv->get(i, numsnaxel-3)*Vn.getCol() ) ;
	}
	
	delete ArInv ;
	
	imageCentered() ;
	return NOERROR ;