mrf.cpp

有指导的马尔可夫随机场（MRF）的图像分割代码

// mrf.cpp: implementation of the Cmrf class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "MRF_imgSeg.h"
#include "mrf.h"

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

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

Cmrf::Cmrf()
{
	CvRect rect;
	rect.x=0;rect.y=0;rect.height=0;rect.width=0;
	m_img=NULL;
	m_classNum=0;
	for(int i=0;i<30;i++)
		m_rect[i]=rect;
	
}

Cmrf::~Cmrf()
{
	if (m_pClass!=NULL) {
		delete []m_pClass;
	}
	if (m_pEnergy!=NULL) {
		delete []m_pEnergy;
	}

}

void Cmrf::CalculMeanAndDeviation()
{
	int k,i,j;

	/* Compute mean and variance for given regions
	*/
		
	for( k=0;k<m_classNum;k++)
	{
		double sum=0;
		int n=0;
		for(i=m_rect[k].x;i<m_rect[k].x+m_rect[k].width;i++)
		{
			for( j=m_rect[k].y;j<m_rect[k].y+m_rect[k].height;j++)
			{
				uchar* ptr1=&CV_IMAGE_ELEM(m_img,uchar,j,i);
				sum+=ptr1[0];
				n++;
			}
		}
		m_mean[k]=sum/n;
		sum=0;
		for(i=m_rect[k].x;i<m_rect[k].x+m_rect[k].width;i++)
		{
			for(j=m_rect[k].y;j<m_rect[k].y+m_rect[k].height;j++)
			{
				uchar* ptr1=&CV_IMAGE_ELEM(m_img,uchar,j,i);
				sum+=(ptr1[0]-m_mean[k])*(ptr1[0]-m_mean[k]);
			}
		}
		m_deviation[k]=sum/n;		
	}

		
}

void Cmrf::SetMember(IplImage *image, CvRect rect[30], int number)
{
	int i=0;
	m_classNum=number;
	for(i=0;i<number;i++)
	{
		m_rect[i]=rect[i];
	}
	if (m_img!=NULL) {
		cvReleaseImage( &m_img ); 
	}
	m_img=cvCreateImage(cvSize(image->width,image->height),8,1);
	m_result=cvCreateImage(cvSize(image->width,image->height),8,1);
	if (image->nChannels==3) {
		cvCvtColor(image,m_img,CV_RGB2GRAY);
	}
	else
	{
		cvCopy(image,m_img,NULL);
	}
//	cvShowImage("Result",m_img);
	m_pClass=new int[image->width*image->height];//(uchar *)malloc(image->width*image->height*8);
	m_pEnergy=new float[image->width*image->height];//float*)malloc(image->width*image->height*8);
	m_beta=0.9;
	CalculMeanAndDeviation();
	// srand((unsigned)time(NULL));
//	for(i=0;i<image->width*image->height;i++)
	//	{
	//		//int classNo=rand()%number;
	//		//		m_pClass[i]=classNo;
	//		
	//		
	//	}
	
}

//DEL float Cmrf::CalculEnergy(CvPoint pix)
//DEL {
//DEL 	float u1,u2;
//DEL 	
//DEL }

double Cmrf::Singleton(int i, int j, int label)
{
	uchar* ptr1=&CV_IMAGE_ELEM(m_img,uchar,j,i);
	return log(sqrt(2.0*3.141592653589793*m_deviation[label])) +
		pow((double)ptr1[j*m_img->width+i]-m_mean[label],2)/(2.0*m_deviation[label]+0.0000001);
	
}

double Cmrf::Doubleton(int i, int j, int label)
{
	double energy = 0.0;
	int h=m_img->height;
	int w=m_img->width;
	
	if (i!=h-1) // south
    {
		if (label == m_pClass[(i+1)+j*w]) energy -= m_beta;
		else energy += m_beta;
    }
	if (j!=w-1) // east
    {
		if (label == m_pClass[i+(j+1)*w]) energy -= m_beta;
		else energy += m_beta;
    }
	if (i!=0) // nord
    {
		if (label == m_pClass[(i-1)+j*w]) energy -= m_beta;
		else energy += m_beta;
    }
	if (j!=0) // west
    {
		if (label == m_pClass[i+(j-1)*w]) energy -= m_beta;
		else energy += m_beta;
    }
	return energy;
	
}

double Cmrf::CalculateEnergy()
 {
 	double singletons = 0.0;
 	double doubletons = 0.0;
	int i, j, k=0;
 	for (i=0; i<m_img->width; ++i)
 		for (j=0; j<m_img->height; ++j)
 		{
 			k = m_pClass[i+j*m_img->width];
			// singleton
 			singletons += Singleton(i,j,k);
 			// doubleton
 			doubletons += Doubleton(i,j,k); // Note: here each doubleton is
			// counted twice ==> divide by
 			// 2 at the end!
 		}
 		return singletons + doubletons/2; 
		
 }

double Cmrf::LocalEnergy(int i, int j, int label)
{
	return Singleton(i,j,label) + Doubleton(i,j,label);
}

void Cmrf::Gibbs()
{
	InitialData();
	int i, j;
	double *Ek;		       // array to store local energies
	int s;
	double summa_deltaE;
	double sumE;
	double z;
	double r;
	double E_old,E;
	int K = 0;
	double T = 4;
	double c = 0.98;
	double t=0.05;
	int w = m_img->width;
	int h = m_img->height;
	//TRandomMersenne rg(time(0)); // make instance of random number generator
	srand((unsigned)time( NULL ));


	Ek = new double[m_classNum];
	
	
	E_old = CalculateEnergy();
	
	do
    {
		summa_deltaE = 0.0;
		for (i=0; i<w; i++)
			for (j=0; j<h; j++)
			{
				sumE = 0.0;
				for (s=0; s<m_classNum; s++)
				{
					Ek[s] = exp(-LocalEnergy(i, j, s)/T);
					sumE += Ek[s];
				}
				r = (float)rand()/RAND_MAX;	// r is a uniform random number
				z = 0.0;
				for (s=0; s<m_classNum; s++)
				{
					z += Ek[s]/sumE; 
					if (z > r) // choose new label with probabilty exp(-U/T).
					{
						m_pClass[j*w+i] = s;
						break;
					}
				}
			}
			E = CalculateEnergy();
			summa_deltaE += fabs(E_old-E);
			E_old = E;
			
			T *= c;         // decrease temperature
			++K;	      // advance iteration counter
			CreateOutput(); // display current labeling
		
    } while (summa_deltaE > t); // stop when energy change is small
	
	delete Ek;
	
}

void Cmrf::InitialData()
{
	/* initialize using Maximum Likelihood (~ max. of singleton energy)
	*/
	int i,j,k;
	double e,e2;
	for(i=0;i<m_img->width;i++)
	{
		for(j=0;j<m_img->height;j++)
		{
			e=Singleton(i,j,0);
			m_pClass[i+j*m_img->width]=0;
			for(k=0;k<m_classNum;k++)
			{
				if ((e2=Singleton(i, j, k)) < e)
				{
					e = e2;
					m_pClass[i+j*m_img->width] = k;
				}
			}
						
		}		
	}
}

void Cmrf::CreateOutput()
{
	int i,j;
	int h=0;
	for(i=0;i<m_result->width;i++)
	{
		for(j=0;j<m_result->height;j++)
		{
			uchar* ptr1=&CV_IMAGE_ELEM(m_result,uchar,j,i);
			ptr1[0]=(uchar)(m_pClass[i+j*m_result->width]*255/m_classNum);
		}
	}
	cvShowImage("Result",m_result);
	cvWaitKey(1);
//	AfxGetMainWnd()->SendMessage(WM_REDRAW);
}

IplImage* Cmrf::GetResultImg()
{
	IplImage * img;
	img=cvCreateImage(cvGetSize(m_result),m_result->depth,m_result->nChannels);
	cvCopy(m_result,img);
	return img;
}