activecontour.cpp

刚上传内容的相关CODEC不能单独上传。于是

// ActiveContour.cpp: implementation of the CActiveContour class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "magicscissors.h"
#include "ActiveContour.h"

#include "EdgeList.h"
#include "EdgePoint.h"
#include "Image/Image.h"
//#include "Image/EdgeDetector.h"

#include <math.h>

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CActiveContour::CActiveContour()
{
	m_pImage = NULL;
	m_pGradientImage = NULL;
	m_pPointList = NULL;
	m_pPolarList = NULL;

	m_uNumofPoint = 0;
	m_dExternalScaler = 0;
	m_dInternalScaler = 0;
	m_dCentripedralScaler = 0;
}

CActiveContour::~CActiveContour()
{
	if( m_pImage )
		delete m_pImage;
	if( m_pPointList )
		delete m_pPointList;
	if( m_pPolarList )
		delete m_pPolarList;
}

//////////////////////////////////////////////////////////////////////
// Public Member Functions Implementation
//////////////////////////////////////////////////////////////////////

BOOL CActiveContour::Create( CEdgeList* pEdgeList )
{
	m_uNumofPoint = (UINT)(pEdgeList->GetLength()) + 2;

	m_pPointList = new CARTCoord[m_uNumofPoint];
	m_pPolarList = new POLARCoord[m_uNumofPoint];

	if( !m_pPointList )	return FALSE;
	if( !m_pPolarList )	return FALSE;

	CEdgePoint* point;
	CEdgePoint* temp;

	point = pEdgeList->GetPoint(m_uNumofPoint-3);	// Get the last point
	m_pPointList[0].x = point->x;
	m_pPointList[0].y = point->y;

	for( UINT i = 1 ; i < m_uNumofPoint-1 ; i++ )
	{
		temp = pEdgeList->GetPoint(i-1);
		m_pPointList[i].x = temp->x;
		m_pPointList[i].y = temp->y;
	}

	point = pEdgeList->GetPoint(0);	// Get the first point
	m_pPointList[m_uNumofPoint-1].x = point->x;
	m_pPointList[m_uNumofPoint-1].y = point->y;

	// Polar transform
	if( !CartToPolar( m_pPointList, m_pPolarList, &m_Centroid, m_uNumofPoint ) )
		return FALSE;
 
	return TRUE;
}

BOOL CActiveContour::FeedEdgeImage( CImage* pImage )
{
	if( !pImage )
		return FALSE;

	if( m_pImage )
	{
		delete m_pImage;
		m_pImage = NULL;
	}
	m_pImage = new CImage( pImage );

	if( !m_pImage )
		return FALSE;

	// Lowpass filtering to the input edge image
	Smoothing( 1.0 );

	// Calculate the garadient of smoothed edge image
	if( !Gradient() )
		return FALSE;

	return TRUE;
}

BOOL CActiveContour::Evolve( double dErrorCapability )
{
	if( !m_pImage )
		return FALSE;

	if( !m_pGradientImage )
		return FALSE;

	if( !m_pPolarList )
		return FALSE;

	CARTCoord temp;
	int		nX, nY;
	double	dTempAngle, dTempRadius;
	double	angle, radius;
	double	dX, dY;
	double	dError = 0;

	do
	{
		dError = 0;
		for( UINT j = 1 ; j < m_uNumofPoint-1 ; j++ )
		{
			// Elastic force
			if( fabs( m_pPolarList[j-1].angle - m_pPolarList[j+1].angle ) > 180 )
			{
				angle = 0.5*m_pPolarList[j-1].angle + 0.5*m_pPolarList[j+1].angle - 180;
				//angle = m_pPolarList[j].angle + m_dInternalScaler*angle;

				if( angle > 360 )
					angle = angle - 360;
				else if( angle < 0 )
					angle = 360 + angle;

				dTempAngle = angle;
			}
			else
			{
				angle = 0.5*m_pPolarList[j-1].angle + 0.5*m_pPolarList[j+1].angle;
				dTempAngle = m_pPolarList[j].angle + m_dInternalScaler*(angle - m_pPolarList[j].angle);
			}

			radius = 0.5*m_pPolarList[j-1].radius + 0.5*m_pPolarList[j+1].radius;
			dTempRadius = m_pPolarList[j].radius + m_dInternalScaler*(radius - m_pPolarList[j].radius);

			// Centripedral force
			dTempRadius -= m_dCentripedralScaler;

			// External force
			PolarToCart( &m_pPolarList[j], &temp, m_Centroid, 1 );
			Round( temp, &nX, &nY );
			dX = temp.x + m_dExternalScaler*m_pGradientImage[nY][nX].x;
			dY = temp.y + m_dExternalScaler*m_pGradientImage[nY][nX].y;

			radius = sqrt( pow(dX - m_Centroid.x,2) + pow(dY - m_Centroid.y,2) );
			angle = atan2( dY - m_Centroid.y , dX - m_Centroid.x );
			angle = angle*180/3.141592;
			if( angle < 0 )
				angle = 360 + angle;

			dTempAngle = dTempAngle + angle - m_pPolarList[j].angle;
			dTempRadius = dTempRadius + radius - m_pPolarList[j].radius;

			if( dTempAngle < 0 )
				dTempAngle = 360 + dTempAngle;
			else if( dTempAngle > 360 )
				dTempAngle = dTempAngle - 360;

			dError += sqrt( pow(m_pPolarList[j].angle - dTempAngle, 2)
				+ pow(m_pPolarList[j].radius - dTempRadius, 2) );
			m_pPolarList[j].angle = dTempAngle;
			m_pPolarList[j].radius = dTempRadius;

			dTempAngle = 0;
			dTempRadius = 0;
		}
		m_pPolarList[0].angle = m_pPolarList[m_uNumofPoint-2].angle;
		m_pPolarList[m_uNumofPoint-1].angle = m_pPolarList[1].angle;
		m_pPolarList[0].radius = m_pPolarList[m_uNumofPoint-2].radius;
		m_pPolarList[m_uNumofPoint-1].radius = m_pPolarList[1].radius;

	}while( dError > dErrorCapability );

	return TRUE;
}

BOOL CActiveContour::SetParameter( double dInternalScaler, double dExternalScaler, double dCentripedralScaler )
{
	m_dInternalScaler = dInternalScaler;
	m_dExternalScaler = dExternalScaler;
	m_dCentripedralScaler = dCentripedralScaler;

	return TRUE;
}

BOOL CActiveContour::GetActiveContour( CEdgeList* pEdgeList )
{
	if( !pEdgeList )
		return FALSE;

	if( !m_pPolarList )
		return FALSE;

	if( m_pPointList )
	{
		delete m_pPointList;
		m_pPointList = NULL;
	}

	m_pPointList = new CARTCoord[m_uNumofPoint];
	if( !m_pPointList )
		return FALSE;

	int temp = pEdgeList->GetLength();

	if( !PolarToCart( m_pPolarList, m_pPointList, m_Centroid, m_uNumofPoint ) )
		return FALSE;

	CEdgePoint*	point;
	for( UINT i = 1 ; i < m_uNumofPoint - 1 ; i++ )
	{
		point = pEdgeList->GetPoint(i-1);
		point->x = (int)m_pPointList[i].x;
		point->y = (int)m_pPointList[i].y;		
	}

	return TRUE;
}

//////////////////////////////////////////////////////////////////////
// Protected Member Functions Implementation
//////////////////////////////////////////////////////////////////////

BOOL CActiveContour::Smoothing( double param )
{
	if( !m_pImage )	return FALSE;
	int w = (int)m_pImage->GetXSize();
	int h = (int)m_pImage->GetYSize();

	int i,j;
	
	double** y1;
	double** y2;
	double** rx;
	double** ry;

	y1 = new double*[h];
	y2 = new double*[h];
	rx = new double*[h];
	ry = new double*[h];

	for( i = 0 ; i < h ; i++ )
	{
		y1[i] = new double[w];		
		y2[i] = new double[w];
		rx[i] = new double[w];
		ry[i] = new double[w];
	}

	BYTE** img = (BYTE**)m_pImage->GetPtr();
	
	for( i = 0 ; i < w ; i++ )
	{

		y1[0][i] = 0.f;		y1[1][i] = 0.f;
		y1[h-1][i] = 0.f;	y1[h-2][i] = 0.f;
		y2[0][i] = 0.f;		y2[1][i] = 0.f;
		y2[h-1][i] = 0.f;	y2[h-2][i] = 0.f;		
	}	
	for( i = 0 ; i < h ; i++ )
	{
		y1[i][0] = 0.f;		y1[i][1] = 0.f;
		y1[i][w-1] = 0.f;	y1[i][w-2] = 0.f;
		y2[i][0] = 0.f;		y2[i][1] = 0.f;
		y2[i][w-1] = 0.f;	y2[i][w-2] = 0.f;		
	}

	double k;
	k = (double)((1-exp(-param))*(1-exp(-param))/(1+2*param*exp(-param)-exp(-2*param)));

	double a1, a2, a3, a4, a5, a6, a7, a8, c1, c2, b1, b2;

	a1  = a5 = k;
	a2 = a6 = k*exp(-param)*(param-1);
	a3 = a7 = k*exp(-param)*(param+1);
	a4 = a8 = -k*exp(-2*param);
	c1 = c2 = 1;

	b1 = (double)(2*exp(-param));
	b2 = (double)(-exp(-2*param));

	int m,n;
	for( j = 2 ; j < h-2 ; j++ )
	{
		for( i = 2 ; i < w-2 ; i++ )
		{	
			y1[j][i] = a1*img[j][i]+a2*img[j][i-1]+b1*y1[j][i-1]+b2*y1[j][i-2];
			n = w-1-i;
			y2[j][n] = a3*img[j][n+1]+a4*img[j][n+2]+b1*y2[j][n+1]+b2*y2[j][n+2];
		}
	}

	for( j = 0 ; j < h ; j++ )
	{
		for( i = 0 ; i < w ; i++ )
			rx[j][i] = y1[j][i]+y2[j][i];
	}

	for( i = 2 ; i < w-2 ; i++ )
	{
		for( j = 2 ; j < h-2 ; j++ )
		{
			y1[j][i] = a5*rx[j][i]+a6*rx[j-1][i]+b1*y1[j-1][i]+b2*y1[j-2][i];
			m = h-1-j;
			y2[m][i] = a7*rx[m+1][i]+a8*rx[m+2][i]+b1*y2[m+1][i]+b2*y2[m+2][i];
		}
	}

	for( j = 0 ; j < h ; j++ )
	{
		for( i = 0 ; i < w ; i++ )
			rx[j][i] = y1[j][i]+y2[j][i];
	}

	for( j = 0 ; j < h ; j++ )
	{
		for( i = 0 ; i < w ; i++ )
		{
			if( rx[j][i] < 0.f )
				img[j][i] = 0;
			else if( rx[j][i] > 255 )
				img[j][i] = 255;
			else
				img[j][i] = (BYTE)rx[j][i];
		}
	}

	for( i = 0 ; i < h ; i++ )
	{
		delete [] y1[i];
		delete [] y2[i];
		delete [] rx[i];
		delete [] ry[i];
	}
	delete [] y1;
	delete [] y2;
	delete [] rx;
	delete [] ry;

	return TRUE;
}

BOOL CActiveContour::Gradient()
{
	if( !m_pImage )	return FALSE;

	int i, j;
	int	nX = m_pImage->GetXSize();
	int	nY = m_pImage->GetYSize();

	BYTE** ptr = m_pImage->GetPtr();

	// 涝仿登绰 捞固瘤绰 亲惑 鞍篮 农扁扼绊 啊沥茄促.
	if( !m_pGradientImage )
	{
		m_pGradientImage = new CARTCoord*[nY];
		if( !m_pGradientImage )	return FALSE;

		for( i = 0 ; i < nY ; i++ )
		{
			m_pGradientImage[i] = new CARTCoord[nX];
			ASSERT( m_pGradientImage[i] != NULL );
		}

		// Initialize the gradient map
		for( i = 0 ; i < nY ; i++ )
		{
			for( j = 0 ; j < nX ; j++ )
			{
				if( i == 0 || i == nY-1 || j == 0 || j == nX-1 )
				{
					m_pGradientImage[i][j].x = 0;
					m_pGradientImage[i][j].y = 0;
				}
			}
		}
	}

	// Horizontal gradient filter : [ -0.5 0 0.5 ]
	// Vertical gradient filter : transpose of [ -0.5 0 0.5 ]
	for( i = 1 ; i < nY-1 ; i++ )
	{
		for( j = 1 ; j < nX-1 ; j++ )
		{
			m_pGradientImage[i][j].x = -0.5*ptr[i][j-1] + 0.5*ptr[i][j+1];
			m_pGradientImage[i][j].y = -0.5*ptr[i-1][j] + 0.5*ptr[i+1][j];			
		}
	}

	return TRUE;
}

BOOL CActiveContour::CartToPolar( CARTCoord* pPointList, POLARCoord* pPolarList, CARTCoord* pCentroid, UINT num )
{
	if( !pPointList )	// Access denied
		return FALSE;

	if( !pPolarList )	// Access denied
		return FALSE;

	// Calculate the centroid of the Active contour
	double	tempX, tempY;
	double	radius;
	double	angle;
	UINT	i;

	tempX = tempY = 0;
	for( i = 0 ; i < num ; i++ )
	{
		tempX += pPointList[i].x;
		tempY += pPointList[i].y;
	}

	tempX /= num;
	tempY /= num;

	pCentroid->x = tempX;
	pCentroid->y = tempY;

	for( i = 0 ; i < num ; i++ )
	{
		radius = sqrt(pow(pPointList[i].x - pCentroid->x,2) + pow(pPointList[i].y - pCentroid->y,2));
		angle = atan2( pPointList[i].y - pCentroid->y, pPointList[i].x - pCentroid->x );
		// radian to degree
		angle = angle*180/3.141592;
		if( angle < 0 )
			angle = 360 + angle;

		pPolarList[i].radius = radius;
		pPolarList[i].angle = angle;
	}

	return TRUE;
}

BOOL CActiveContour::PolarToCart( POLARCoord* pPolarList, CARTCoord* pPointList, CARTCoord Centroid, UINT num )
{
	if( !pPointList )	// Access denied
		return FALSE;

	if( !pPolarList )	// Access denied
		return FALSE;

	UINT	i;
	double	radian, angle;

	for( i = 0 ; i < num ; i++ )
	{
		angle = pPolarList[i].angle;
		if( angle > 180 )
			angle = angle - 360;

		radian = angle*3.141592/180;
		pPointList[i].x = pPolarList[i].radius*cos(radian) + Centroid.x;
		pPointList[i].y = pPolarList[i].radius*sin(radian) + Centroid.y;
	}

	return TRUE;
}

BOOL CActiveContour::Round( CARTCoord coord, int* X, int* Y )
{
	double x = coord.x;
	double y = coord.y;

	if( x >= 0 )
	{
		if( x - (int)x >= 0.5 )
			*X = (int)x + 1;
		else
			*X = (int)x;
	}
	else
	{
		if( x - (int)x <= 0.5 )
			*X = (int)x -1;
		else
			*X = (int)x;
	}

	if( y >= 0 )
	{
		if( y - (int)y >= 0.5 )
			*Y = (int)y + 1;
		else
			*Y = (int)y;
	}
	else
	{
		if( y - (int)y <= 0.5 )
			*Y = (int)y -1;
		else
			*Y = (int)y;
	}

	if( *X < 0 )
		*X = 0;
	if( *Y < 0 )
		*Y = 0;

	return TRUE;
}