gsnake.c

snake算法

		 now->getLambda() );
        }

        return Epoint;
}

/*
 * compute current snaxel energy based on 3 points
 */
double GSNAKE::ESnaxel( SNAXEL *now,	 	/* current snaxel */
                        int pt_ind,       	/* point index */
                        int loc_ind,      	/* location index */
                        int spacing,      	/* desired spacing */
                        LOCATIONS *searchLocs,  /* search locations */
                        int verbose )
{
	/* gaussian edge and image */
        register short i;
        double Emodel;          /* internal energy */
        double Eimage;          /* external energy */
        double Epoint;          /* snaxel energy */

        /* check out of range */
        if (now->getRow() != CAP(0, now->getRow(), GaussImg->getRow()-1) ||
            now->getCol() != CAP(0, now->getCol(), GaussImg->getCol()-1) )
            return 1.0;

        Emodel = EInternal( now );

        double MinEImage = 1.0;

        for (i=-spacing/2; i<=spacing/2; i++) {

            now->putRow( searchLocs->getRow( pt_ind, loc_ind ) + 
			 i*searchLocs->getNormalY( pt_ind ) );
            now->putCol( searchLocs->getCol( pt_ind, loc_ind ) + 
			 i*searchLocs->getNormalX( pt_ind ) );

	    Eimage = EExternal(now);
	   
           MinEImage = MIN( Eimage, MinEImage );

        }

        now->putCol( searchLocs->getCol( pt_ind, loc_ind ) );
        now->putRow( searchLocs->getRow( pt_ind, loc_ind ) );
        Eimage = MinEImage;

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

        return Epoint;
}

/*
 * energy : edge force
 */
double GSNAKE::EEdge( SNAXEL *now,         /* current snaxel */
		      EDGE *edgeMap )     /* pyramid edge map */
{
        double mag;             /* edge magnitude */
        double Amodel;          /* gsnake object angle */
        double Adata;           /* edge map angle */
        double len_W;
        int mulp = 1;
        double uv_x, uv_y, uv0_x, uv0_y, uv1_x, uv1_y;

	short row,col;

	

	 /* check out of range */
        if (now->getRow() != CAP(0, now->getRow(), edgeMap->getRow()-1) ||
            now->getCol() != CAP(0, now->getCol(), edgeMap->getCol()-1) )       
            return 1.0;
   
        	
        /* get edge magnitude */
        mag = edgeMap->getMag( ROUNDOFF( now->getRow() ),
                               ROUNDOFF( now->getCol() ) );

        if ( mag < VERY_SMALL)
                return 1.0;

        /* return edge magnitude if opened snake and at first or last pt */
        if ( mode == _OPENED &&
             ( now->getPrev() == NULL || now->getNext() == NULL ) )
        	return 1 - mag;

        /* compute snaxel tangent direction */
        now->getDirectionVec( mode, head, tail, &uv0_x, &uv0_y, &uv1_x, &uv1_y);

        /* tangent vector at current snaxel */
        uv_x = uv0_x + uv1_x;
        uv_y = uv0_y + uv1_y;

        len_W = hypot( uv_x, uv_y );
        if ( len_W < VERY_SMALL )
                return 1.0;

        /* compute image energy */
        Amodel = ATAN2( uv_x, -1*uv_y );
        Adata = edgeMap->getAng( ROUNDOFF( now->getRow() ),
                                 ROUNDOFF( now->getCol() ) );

        if ( direction )
                mulp = EDGEDIR(-1*uv_y, uv_x, cos(Adata), sin(Adata)) *
                        direction;

        return (        /* compute projection angle */
                1 - mulp*mag*fabs(cos(Amodel-Adata)) );
}

double GSNAKE::EMag( SNAXEL *now,         /* current snaxel */
                     EDGE *edgeMap )     /* pyramid edge map */
{
	

	 /* check out of range */
        if (now->getRow() != CAP(0, now->getRow(), edgeMap->getRow()-1) ||
            now->getCol() != CAP(0, now->getCol(), edgeMap->getCol()-1) )
            return 1.0;

	short row = ROUNDOFF( now->getRow() );
	short col = ROUNDOFF( now->getCol() ); 

	    
	return ( MAXMAG - edgeMap->getMag(row, col) );
}

double GSNAKE::EInten( SNAXEL *now,		/* current snaxel */
		       IMAGE *gaussImg )	/* gaussian image */
{
	 /* check out of range */
        if (now->getRow() != CAP(0, now->getRow(), gaussImg->getRow()-1) ||
            now->getCol() != CAP(0, now->getCol(), gaussImg->getCol()-1) )
            return 1.0;

	short row = ROUNDOFF( now->getRow() );
	short col = ROUNDOFF( now->getCol() ); 


	return ( MAXMAG - gaussImg->get(row, col) );
}

double GSNAKE::EInternal(SNAXEL *now)
{
	if (now)
	{
		Eint = CONTOUR::EInternal(now);
		return Eint;
	
	} else return PARAMERROR;
}


/* Top level external energy method. Will return PARAMERROR if null
   oassed in. */
double GSNAKE::EExternal(SNAXEL *now)
{  
	double Eimage;

     	if (now) {	
     		
		switch( EextType ) {

	   		case _EDGE :
		   	 	Eimage = EEdge(now, EdgeMap);
             	   		break;

	   		case _EDGEMAG :
				Eimage = EMag( now, EdgeMap);			
                		break;

	   		case _INTENSITY :
				Eimage = EInten(now, GaussImg);
				break;
	   	}
	   	
	   	return Eimage;

	} else 
		return PARAMERROR;	
}


/*
 * regularizetion : minmax principle or global/local lambda
 */
void GSNAKE::regularize( SNAXEL *now,	/* current snaxel */
			     double *Emodel,    /* model energy */
                             double *Eimage,    /* image energy */
                             double *Epoint )   /* snaxel energy */
{
        if ( global_lambda == _LOCAL_MINMAX )   /* local minimax */
                *Epoint = MAX(*Emodel, *Eimage) + 0.0001*MIN(*Emodel, *Eimage);

        else if ( global_lambda == _LOCAL_LAMBDA ) {
                *Emodel *= now->getLambda() / ( 1.0 - now->getLambda() );
                *Epoint = *Emodel + *Eimage;
           }

        /* combine Emodel and Eedge by global model parameter */
	else if (( global_lambda < 1.0 ) && (global_lambda > 0.0) )
                *Epoint = global_lambda / ( 1.0 - global_lambda ) * 
			*Emodel + *Eimage;
	else
		*Epoint = 1;
}


/* Performs duplication of snake. If target exists, will copy contour and
   parameters. If target does not exist, will create target, copy pyramid 
   and contours and parameters. 
     */
   
GSNAKE *GSNAKE::duplicate(GSNAKE* target)
{
	GSNAKE *dupSnake;
	
	/* If no existing snake create new one */
	if (!target){
	   if ( !( dupSnake = new GSNAKE ) ) {
               fprintf(stderr," <GSNAKE::duplicate> : memory allocation 	
               	               error\n");
               return NULL;
            }
	    PYRAMID::duplicate(dupSnake);
	 }
	 else dupSnake = target;   /* else assign pointer to target */

	CONTOUR::duplicate(dupSnake);
		
	/* Copy GSNAKE parameters */
	dupSnake->global_lambda = global_lambda;
	dupSnake->EextType      = EextType;
 	dupSnake->Eext          = Eext;
 	dupSnake->Esnake        = Esnake; 

	return dupSnake;
}


/*
 * deform contour by mouse clicking
 */
void GSNAKE::deform(unsigned char blowup, short expand)
{
	register short i;
	SNAXEL *sxptr;
		   
	double xdata[numsnaxel];
	double ydata[numsnaxel];
	
			    
	for ( i=0, sxptr=getHead(); i<numsnaxel && sxptr;
		i++, sxptr=sxptr->getNext() ) {

		xdata[i] = sxptr->getCol()*blowup;
		ydata[i] = sxptr->getRow()*blowup;
	}

	show(blowup,expand,1);								       
	xwin_ShapeContour( rawImg->ximg, ydata, xdata, numsnaxel );
									
	CONTOUR::init( ydata, xdata, numsnaxel, blowup); 
}

/*
 * draw snake with given offset
 */
void GSNAKE::showLine( unsigned char blowup,	/* blow up ratio */
		   int Xoffset,			/* x offset */
		   int Yoffset )		/* y offset */
{
        SNAXEL *sxptr;
	short nowrow, nowcol;
	short prevrow, prevcol;
	short register i;

	prevcol = ROUNDOFF( getHead()->getCol() ) + Xoffset;
	prevrow = ROUNDOFF( getHead()->getRow() ) + Yoffset;

        for ( i = ( mode == _CLOSED ? 0 : 1 ),
		sxptr = getHead()->getNext(); ( i < numsnaxel ) && sxptr; 
		i++, sxptr=sxptr->getNext(mode, head) ) {

                nowcol = ROUNDOFF( sxptr->getCol() ) + Xoffset;
                nowrow = ROUNDOFF( sxptr->getRow() ) + Yoffset;

		xwin_DrawLine( prevcol, prevrow, nowcol, nowrow, blowup );

		prevcol = nowcol;
		prevrow = nowrow;

        }
}

/*
 * draw snake point
 */
void GSNAKE::show(  unsigned char blowup, 	/* blow up ratio */
		   short expand,                /* expand ratio */
               	   int showImg,			/* show image ? */
		   int img_Xoffset, 		/* image offset */
		   int img_Yoffset,
		   int pt_Xoffset,
		   int pt_Yoffset )
{
        SNAXEL *sxptr;
        
        if(!rawImg) {
          	fprintf(stderr, "GSNAKE does not contain a raw image.") ;
          	fprintf(stderr, " Using CONTOUR::display()\n") ;
          	CONTOUR::display(blowup) ;
        }
        
        else {

	 if ( showImg )
		rawImg->show( blowup, 1, img_Xoffset, img_Yoffset );

	 else 	
          for ( sxptr = getHead(); sxptr; sxptr=sxptr->getNext() )
		xwin_clrpoint( rawImg->ximg, 
			     (ROUNDOFF( sxptr->getCol() )+ img_Xoffset)*blowup, 
			     (ROUNDOFF( sxptr->getRow() )+ img_Yoffset)*blowup);


	 CONTOUR::show(rawImg,blowup, pt_Xoffset,pt_Yoffset, expand);
	}

}

/*
 * display new location of snaxel
 */
void GSNAKE::move(  GSNAKE *old_gsnake,	/* old gsnake */
		   int pt_Xoffset,	/* point offset */
		   int pt_Yoffset )
{
	SNAXEL *old_sxp;
	SNAXEL *new_sxp;

	for( new_sxp = getHead(), old_sxp = old_gsnake->getHead(); 
	     new_sxp && old_sxp; new_sxp = new_sxp->getNext(),
			 	 old_sxp = old_sxp->getNext() ) {

		short Col = ROUNDOFF( old_sxp->getCol() + pt_Xoffset) ;
                short Row = ROUNDOFF( old_sxp->getRow() + pt_Yoffset) ;

		new_sxp->show( rawImg, Row, Col, pt_Xoffset,pt_Yoffset );
	}
}