snaxel.c

snake算法

/****************************************************************************
 File Name  : snaxel.c
 Purpose    : provide basic routines for manipulation of snaxel angle
              and energy
 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"
#include "random.h"

SNAXEL::SNAXEL()
{
	next = NULL;
	prev = NULL;
	row = col = alpha = beta = lambda = Eint = Eext = Esnaxel = 0 ;
}

/*
 * get the mean position (Eint = 0) for a point on a snake
 */
void SNAXEL::meanPosition( SNAKEMODE mode,      /* open / closed snake */
                           double cg_row,       /* center of gravity */
                           double cg_col,
                           SNAXEL *head,        /* head of snake */
                           SNAXEL *tail,        /* tail of snake */
                           double *meanrow,     /* mean position */
                           double *meancol)
{
        SNAXEL *Prev = getPrev(mode, tail);
        SNAXEL *Next = getNext(mode, head);

        /* compute mean position */
        *meanrow = getAlpha() * ( Prev->getRow() - cg_row  ) +
               getBeta()  * ( Next->getRow() - cg_row  ) +  cg_row;
        *meancol = getAlpha() * ( Prev->getCol() - cg_col  ) +
               getBeta()  * ( Next->getCol() - cg_col  ) + cg_col;
}

/*
 * compute angle at current snaxel
 */
double SNAXEL::getNormalAng( SNAKEMODE mode,        /* opened / closed gsnake */
                         SNAXEL *head,          /* head of snake */
                         SNAXEL *tail )         /* tail of snake */
{
        double uv0_x, uv0_y;    /* unit vector */
        double uv1_x, uv1_y;    /* unit vector */
        double angle;

        if ( mode == _OPENED && prev == NULL )
                angle = ATAN2( (getCol() - next->getCol()),
                                  (getRow() - next->getRow()) );


        else if ( mode == _OPENED && next == NULL )
                angle = ATAN2( (getCol() - prev->getCol()),
                                  (getRow() - prev->getRow()) );

        else {
                getDirectionVec(           /* compute tangent vector */
                         mode, head, tail,
                         &uv0_x, &uv0_y, &uv1_x, &uv1_y);

                angle = ATAN2(     /* compute gradient */
                        (uv0_x + uv1_x),
                       -1*(uv0_y + uv1_y) );
        }

        return angle;
}



/*
 * compute tangent vector of snaxel
 */
void SNAXEL::getDirectionVec(
                SNAKEMODE mode,                 /* opened / closed gsnake */
                SNAXEL *head,                   /* head of snake */
                SNAXEL *tail,                   /* tail of snake */
                double *Uv0_x,double *Uv0_y,    /* unit vectors */
                double *Uv1_x,double *Uv1_y )
{
        double v0_x, v0_y, v1_x, v1_y;          /* directional vectors */
        double len_v0, len_v1;                  /* length of v0 and v1 */
        double m0, m1;                          /* scalar factor */

        SNAXEL *Next = getNext( mode, head);
        SNAXEL *Prev = getPrev( mode, tail);

        /* compute directional vector */
        v0_x = getCol() - Prev->getCol();
        v0_y = getRow() - Prev->getRow();
        v1_x = Next->getCol() - getCol();
        v1_y = Next->getRow() - getRow();

        /* compute vector length */
        len_v0 = hypot( v0_x, v0_y );
        len_v1 = hypot( v1_x, v1_y );

        /* compute unit vector */
        m0 = ((len_v0) ? 1.0/len_v0 : 0.0);
        m1 = ((len_v1) ? 1.0/len_v1 : 0.0);
        *Uv0_x = m0 * v0_x;
        *Uv0_y = m0 * v0_y;
        *Uv1_x = m1 * v1_x;
        *Uv1_y = m1 * v1_y;
}

/*
 * compute normal vector of current snaxel
 */
void SNAXEL::getNormalVec( SNAKEMODE mode,         /* opened / closed gsnake */
                        SNAXEL *head,           /* head of snake */
                        SNAXEL *tail,           /* tail of snake */
                        double *nv_x,           /* normal vectors */
                        double *nv_y )
{
        double uv0_x, uv0_y;    /* unit vector */
        double uv1_x, uv1_y;    /* unit vector */
        double dx, dy;
        double len ;

 
        getDirectionVec(           /* compute tangent vector */
                         mode, head, tail,
                         &uv0_x, &uv0_y, &uv1_x, &uv1_y);
        dx = uv0_x + uv1_x;
        dy = uv0_y + uv1_y;
        
        /* use infinity norm */
        len = MAX(fabs(dx), fabs(dy));

        if ( len < 0.001 ) {
                double angle = rand_uniform( 0.0, 2*M_PI );
                dx = cos( angle );
                dy = sin( angle );
                len = MAX( fabs(dx), fabs(dy) );
        }

 
        /* rotate tangent vectors by 90 degrees */
        *nv_x = -dy/len;
        *nv_y = dx/len;
}

/*
 * get next snaxel
 */
SNAXEL *SNAXEL::getNext( SNAKEMODE mode,        /* opened / closed gsnake */
                         SNAXEL *head )         /* head of snake */
{
        SNAXEL *Next;

        if ( mode == _OPENED ) { /* if opened snake */
                Next = ( next == NULL ) ?
                getPrev()->getPrev() :
                getNext();
        }
        else { /* if closed snake */
                Next = ( next == NULL) ?
                head : getNext();
        }

        return Next;
}

/*
 * get prev snaxel
 */
SNAXEL *SNAXEL::getPrev( SNAKEMODE mode,        /* opened / closed gsnake */
                         SNAXEL *tail )         /* tail of snake */
{
        SNAXEL *Prev;

        if ( mode == _OPENED ) { /* if opened snake */
                Prev = ( prev == NULL ) ?
                getNext()->getNext() :
                getPrev();
        }
        else { /* if closed snake */
                Prev = ( prev == NULL ) ?
                tail : getPrev();
        }

        return Prev;
}


/* Comments : Point offsets are necessary due to different expansion levels.
              Image offset is fixed but the point offset depends on level 
              This offset will be added within the X windows move.*/

void SNAXEL::show( IMAGE *img,                     /* ref. image */
                   short Old_Row, short Old_Col, /* Destination location */
                   int pt_Xoffset, int pt_Yoffset) /* Point offsets  */
{	

	xwin_MovePoint(img->ximg, Old_Col, Old_Row, 
			ROUNDOFF(col), ROUNDOFF(row),
			pt_Xoffset, pt_Yoffset);
	
}