contour.c

snake算法

} 


/*
 * convert a snake from V (image-centered coordinates) to U (object centered)
 */
void CONTOUR::contourCentered()
{
        SNAXEL *sxptr;
	short positive = 1;


	/* check that the co-ordinates are not already in U form */
	/* if in contour centered form, some co-ordinates will be negative */
	sxptr = head;
	
	while (sxptr && positive)
	{
		if (sxptr->getRow()<0) positive =0;
		sxptr=sxptr->getNext();
	};
	
	/* Co-ordinates are in V form */
	if (positive)
	{
        	for(sxptr=head; sxptr; sxptr=sxptr->getNext()) {
                	sxptr->putRow( sxptr->getRow() - cg_row );
                	sxptr->putCol( sxptr->getCol() - cg_col );
        	};
     	};
}

/*
 * convert a snake from U to V coordinate
 */
void CONTOUR::imageCentered()
{
        SNAXEL *sxptr;
	short negative=0;

	/* check that the co-ordinates are not already in V form */
	sxptr = head;
	
	while ( (sxptr) && (!negative) )
	{
		if (sxptr->getRow()<0) negative =1;
		sxptr=sxptr->getNext();
	};

   	if (negative)
   	{

        	for(sxptr=head; sxptr; sxptr=sxptr->getNext()) {
                	sxptr->putRow( cg_row + sxptr->getRow() );
                	sxptr->putCol( cg_col + sxptr->getCol() );
        	};
	};
}

/*
 * compute average power (length)  of snaxels
 */
double CONTOUR::computeAvgLength()
{
        double sumlen = 0;              /* total power */
        short numpoints;                /* number of snaxel points */
        SNAXEL *sxptr, *next;
        register short i;

        /* compute center of gravity */
        computeCG();

        numpoints = (mode == _OPENED) ? numsnaxel - 1 :
                                        numsnaxel;

        for ( i=0, sxptr=getHead(); i<numpoints;
                sxptr=sxptr->getNext(), i++ ) {

                next = sxptr->getNext( mode, head );

                sumlen += hypot( (double) (sxptr->getRow() - next->getRow()),
                                 (double) (sxptr->getCol() - next->getCol()) );
        }

        avglen = sumlen / (double) numpoints;

        return ( avglen );
}

/*
 * compute center of gravity
 */
void CONTOUR::computeCG()
{
        SNAXEL *sxptr;

        /* reset center of gravity */
        cg_row = 0;
        cg_col = 0;

        for(sxptr=getHead(); sxptr; sxptr=sxptr->getNext()) {
                cg_row += sxptr->getRow();
                cg_col += sxptr->getCol();
        }

        cg_row = cg_row / (double) (numsnaxel);
        cg_col = cg_col / (double) (numsnaxel);
}
/*
 * use Monte Carlo estimation to compute Z : the normalizing constant
 */
void CONTOUR::computeZ()
{
	SNAXEL *sptr, *tSptr ;
	double Esum  ;
	int numSample ;
	float M ;
	double oldZ ;
	
	
	CONTOUR *temp = CONTOUR::duplicate() ;
	
	/* get average lambda */
	double avgLam = 0 ;
	for(sptr=getHead(); sptr; sptr=sptr->getNext())
		avgLam += sptr->getLambda() ;
	avgLam /= getNumSnaxel() ;	
	
	/* compute region for importance sampling and total size */
	M = 3*sqrt(sig_nu_sqr*(1.0/avgLam - 1.0))*computeAvgLength() ;
	
	Esum = 0 ;
	Z = 0 ;
	numSample = 1 ;
	do {
	     oldZ = Z ;
	
	     /* generate random assignment of snake */
	     for(sptr=getHead(), tSptr=temp->getHead() ; sptr && tSptr; 
	     		sptr=sptr->getNext(), tSptr=tSptr->getNext() ) {

	 	tSptr->putRow( sptr->getRow() + rand_uniform(-M, M) ) ;
	     	tSptr->putCol( sptr->getCol() + rand_uniform(-M, M) ) ;
	     }
	     
	     /* calculate internal energy */
	     double _Einternal = 0 ;
	     temp->computeAvgLength() ;
	     for(tSptr=temp->getHead(); tSptr; tSptr=tSptr->getNext())
	     	_Einternal += temp->EInternal(tSptr) / 
	     			(sig_nu_sqr*(1.0/tSptr->getLambda() - 1.0)) ;
	     
	     Esum += exp(-_Einternal) ;
	     
	     
	     Z =  Esum / numSample ;
   	     numSample++; 

   	} while( fabs(Z - oldZ)/Z > 0.001 || numSample < 100 ) ;
   	
    	Z =  ( pow((2*M+1), (double)(getNumSnaxel()) ) * Z) + 1.0 ;
      	
        delete temp ;
}



/*
 * expand snake
 * if spacing >  0 : desired spacing of expanded snake
 * if spacing <= 0 : desired number of snaxels inserted
 */
int CONTOUR::expand( double ratio,         /* ratio of expanding */
             	     int spacing)          /* spacing */
{

        SNAXEL *now, *next;
        double *ydata, *xdata;  /* snaxel coordinate */
        double dist;            /* dist. betw. snaxels */
        register short i, j;
        int numpts;             /* number of snaxels */
        int instNum;            /* # of pts inserted betw. snaxels */
        int counter;
        double dy, dx;

	if ( !spacing ) { /* if no new snaxels is inserted */

	    SNAXEL *sxptr;

	    for ( sxptr = getHead(); sxptr; sxptr = sxptr->getNext() ) {
                sxptr->putRow( sxptr->getRow() * ratio );
                sxptr->putCol( sxptr->getCol() * ratio );
	    }

	    return NOERROR;
        }

        /* allocate more memory to hold expanded locations */
        
	numpts = ROUNDOFF(ratio) * numsnaxel * ROUNDOFF( computeAvgLength() );
	
	
        ydata = (double *)_iMemAlloc( numpts * sizeof(double) );
        xdata = (double *)_iMemAlloc( numpts * sizeof(double) );

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

        /* expand source snake and put the results into temporary buffer.
           The intersnaxel distance in the expanded snake is aprrox
           equal to "spacing" */

        now = getHead();
        ydata[0] = (ratio) * now->getRow();       /* starting point */
        xdata[0] = (ratio) * now->getCol();
        numpts = 1;
        counter = (mode==_CLOSED) ? numsnaxel + 1 : numsnaxel;

        for( j = 1; j < counter && now; j++, now = now->getNext() ) {

                /* find out how many points to be inserted between two
                   successive snaxels */

                next = now->getNext(mode, head);
                dy = (ratio) * ( next->getRow() - now->getRow() );
                dx = (ratio) * ( next->getCol() - now->getCol() );

                dist = MAX( fabs(dx), fabs(dy) );

                instNum = (spacing > 0) ?
                          ROUNDOFF( dist / (double) spacing ) - 1 :
                          -spacing;

                instNum = MIN( instNum, (int) (dist - 1) );

                for( i = 0 ; i < instNum  ; i++ ) {  /* insert the points */
                        ydata[numpts] = (ratio) * now->getRow() +
                                ROUNDOFF((i+1) * dy / (double) (instNum + 1));
                        xdata[numpts] = (ratio) * now->getCol() +
                                ROUNDOFF((i+1) * dx / (double) (instNum + 1));
                        numpts++ ;
                }

                /* put int the next point from scrc snake */
                if( j < numsnaxel ) {
                        ydata[numpts] = (ratio) * next->getRow() ;
                        xdata[numpts] = (ratio) * next->getCol() ;
                        numpts++;
                }
        }

        /* copy temporary buffer to expanded snake */
        if ( init(ydata, xdata, numpts) == MEMORYERROR )
                return MEMORYERROR;

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

        return NOERROR;
}


/*
 * copy snaxel value of snake
 */
void CONTOUR::copySnaxel(CONTOUR *target)
{
        SNAXEL *src;    /* source snaxel */
        SNAXEL *dest;   /* destination snaxel */

        src = getHead();                /* get head of src snake */
        dest = target->getHead();     /* get head of dest snake */

        /* start copying */
        for ( ; src; src = src->getNext(), dest = dest->getNext() ) {
                dest->putAlpha(	src->getAlpha() );
                dest->putBeta(	src->getBeta() );
		dest->putEint(	src->getEint() );
		dest->putEext(	src->getEext() );
                dest->putRow(	src->getRow() );
                dest->putCol(	src->getCol() );
                dest->putLambda(src->getLambda() );
                dest->putEsnaxel(src->getEsnaxel() );
        }
}


/*
 * Allows several functions of duplication
 	1. target NULL -> create new CONTOUR. All will be copied.
 	2. flag 1 -> copy snaxels only to target contour. For target NULL,
 	             all will be copied. This flag is ignored. 	
 	3. flag 0 -> copy all contour data. Will reset target contour and
 	             copy existing contour data
  */

CONTOUR *CONTOUR::duplicate(CONTOUR *target, short flag)
{
        CONTOUR *dupContour;

	/* If no existing contour create new one */
	if (!target)
	{

	   if ( !( dupContour = new CONTOUR ) ) {
               fprintf(stderr," <CONTOUR::duplicate> : memory allocation 	
               	               error\n");
               return NULL;
            }
	 }
	 else dupContour=target;   /* else assign pointer to target */


        /* If flag is 0, copy all data.
           Either creating a new contour or copying everything to target */ 

	if ( (!flag) || (!target) ) {
        	
        	dupContour->reset();
        	if ( dupContour->initSnaxel( numsnaxel ) == MEMORYERROR ) {
           	    fprintf(stderr,"\n<CONTOUR::duplicate> : memory allocation 
           	    `		  error\n");
     	 	    return NULL;
        	}
        	
        	/* copy parameter fields */
        	dupContour->putMode(mode);
        	dupContour->putNumSnaxel(numsnaxel);
        	dupContour->putCgRow(cg_row);
        	dupContour->putCgCol(cg_col);
        	dupContour->putDirection(direction); 
        	dupContour->putSigNuSqr(sig_nu_sqr);
        	dupContour->putZ(Z);            	       		
	}
	
              
        /* Check if target snake is the same as source snake. Otherwise
           unable to copy snaxels over. However, dupContour must still be
           returned, otherwise existing target snake will be lost */
           
        if ( dupContour->getNumSnaxel() == numsnaxel ) 
		copySnaxel( dupContour); 
	else {
	
	     printf("\nTarget snake and source snake have different number 
	                of snaxels. Did not copy.");
	}

   	dupContour->computeAvgLength();
 
 	return dupContour;
}



/*
 * adjust head and tail of gsnake
 */
void CONTOUR::adjust()
{
        SNAXEL *sxptr;

        /* adjust head */
        sxptr = head;
        sxptr->putBeta( 1.0/sxptr->getNext()->getAlpha() );
        sxptr->putAlpha( - (sxptr->getNext()->getBeta() / 
		sxptr->getNext()->getAlpha()) );

        /* adjust tail */
        sxptr = tail;
        sxptr->putAlpha( 1.0/sxptr->getPrev()->getBeta() );
        sxptr->putBeta(  - (sxptr->getPrev()->getAlpha() /
		sxptr->getPrev()->getBeta()) );
}

/*
 * affine invariant transform
 */
void CONTOUR::affineTransform( double sx,	/* x-scale changes */
				double sy,	/* y-scale changes */
				double rt,	/* rotation changes */
				double tx,	/* x-translation changes */
				double ty,	/* y-translation changes */
				double idx,	/* x-dilation changes */
				double idy )	/* y-dilation changes */
{
	SNAXEL *sxptr;
	double row, col;

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

		/* rotation and dilation */
		col = (cos(rt) + sin(rt)*idy)   * sxptr->getCol() +
                      (cos(rt)*idx + sin(rt))   * sxptr->getRow();
                row = ( -sin(rt)+ cos(rt)*idy ) * sxptr->getCol() +
                      ( -sin(rt)*idx + cos(rt)) * sxptr->getRow();

		/* scaling */
		col *= sx;
		row *= sy;

		/* translation */
		col += tx;
		row += ty;

		sxptr->putCol( col );
		sxptr->putRow( row );
	}
}

void CONTOUR::print()
{
	register short i;
        SNAXEL *sxptr;

	printf("numsnaxel = %d, mode = %s, avglen = %.3f\n", numsnaxel,
		( mode == _OPENED ) ? "opened" : "closed", avglen) ;
	printf("sig_nu_sqr = %.3e, Z = %.3f\n", sig_nu_sqr, Z) ;
	printf("COG : (%d %d)\n", ROUNDOFF(cg_row), ROUNDOFF(cg_col) ) ;
        for(i=0, sxptr=head; sxptr; i++, sxptr = sxptr->getNext())
                printf("(%d) col:%d row:%d\talpha:%f\tbelta:%f\tlambda:%f\n",
                        i, ROUNDOFF(sxptr->getCol()), 
			ROUNDOFF(sxptr->getRow()), 
			sxptr->getAlpha(), sxptr->getBeta(),
                        sxptr->getLambda());
	printf("\n");
}


void CONTOUR::show(IMAGE *backgnd, unsigned char blowup,
		   int pt_Xoffset, int pt_Yoffset, short expand)
{
	SNAXEL *sxptr;

	xwin_setexpandratio( expand, blowup );

        for ( sxptr = head; sxptr; sxptr=sxptr->getNext() ) {
            
                sxptr->show( backgnd, ROUNDOFF( sxptr->getRow() ), 
                	     ROUNDOFF( sxptr->getCol() ), pt_Xoffset, 
                	     pt_Yoffset ); 

        }
}


/* Method to display a contour on background. Contour will be shown
   centralised in the image */

void CONTOUR::display(short mag)
{
	SNAXEL *sptr;
	IMAGE background;
	int tx,ty;
	int minX, minY,maxX,maxY,length,height;
	tx=0;ty=0;
	short Row,Col;

	maxX=0;
	maxY=0;					
					
	minX=ROUNDOFF(head->getCol());
	minY=ROUNDOFF(head->getRow());		
	for (sptr=head; sptr;sptr=sptr->getNext())
	{
		tx=ROUNDOFF(sptr->getCol());
		ty=ROUNDOFF(sptr->getRow());
		if (tx<minX) minX = tx ;
			else if (tx>maxX) maxX = tx;
		if (ty<minY) minY = ty ;
			else if (ty>maxY) maxY = ty;	
	}
	
	length = maxX-minX;
	height = maxY-minY;
				
	/* Assign minimum image size */
	if (length<10) length= 30;
	if (height<10) height =30;	

	/* Setting all points with reference to minimum.  */
	for (sptr=head;sptr;sptr=sptr->getNext())
	{
		sptr->putRow(sptr->getRow()-minY+ROUNDOFF(height/2));
		sptr->putCol(sptr->getCol()-minX+ROUNDOFF(length/2));
	
	}
	
	background.init(height*2,length*2);
	for(ty=0;ty<background.getRow();ty++)
		for(tx=0;tx<background.getCol();tx++)
			background.put(ty,tx,20);
					
	background.show(mag);
	xwin_setexpandratio( 1, mag );

        for ( sptr = head; sptr; sptr=sptr->getNext() ) {

                Col = ROUNDOFF( sptr->getCol() ) ;
                Row = ROUNDOFF( sptr->getRow() ) ;
                
                
                sptr->show( &background, Row, Col );

        }

	for (sptr=head;sptr;sptr=sptr->getNext())
	{
		sptr->putRow(sptr->getRow()+minY-ROUNDOFF(height/2));
		sptr->putCol(sptr->getCol()+minX-ROUNDOFF(length/2));
	}

}