【图像算法】彩色图像分割专题八：基于MeanShift的彩色分割

-------------------------------------------------------------------------------------------------------------------------------

SkySeraph July 14th 2011 HQU

Email：zgzhaobo@gmail.com QQ：452728574

Latest Modified Date：July 14th 2011 HQU

-------------------------------------------------------------------------------------------------------------------------------

》原理

不啰嗦了，看这：http://people.csail.mit.edu/sparis/#cvpr07

-------------------------------------------------------------------------------------------------------------------------------

》源码

核心代码(参考网络)

//============================Meanshift==============================//
void MyClustering::MeanShiftImg(IplImage * src , IplImage * dst , float r , int Nmin ,int Ncon )
{
	int i , j , p ,k=0,run_meanshift_slec_number=0;
	int pNmin;								//mean shift产生的特征的搜索框内的特征数
	IplImage * temp , * gray;						//转换到Luv空间的图像
	CvMat * distance , * result , *mask;				//
	CvMat * temp_mat ,*temp_mat_sub ,*temp_mat_sub2 ,* final_class_mat;			//Luv空间的图像到矩阵，图像矩阵与随机选择点之差，
	CvMat * cn ,* cn1 , * cn2 , * cn3;
	double /*covar_img[3] ,*/ avg_img[3];		//图像的协方差主对角线上的元素和，各个通道的均值
	double r1;			//搜索半径
	int	temp_number;
	meanshiftpoint meanpoint[25];		//存储随机产生的25点
	CvScalar	cvscalar1,cvscalar2;
	int	order[25];
	Feature	feature[100];			//特征
	double	shiftor;
	CvMemStorage * storage=NULL;
	CvSeq * seq=0 , * temp_seq=0 , *prev_seq;
//---------------------------------------------RGB to Luv空间,初始化----------------------------------------------
	temp			=	cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U, src->nChannels);
	gray			=	cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U, 1);
	temp_mat		=	cvCreateMat(src->height,src->width,CV_8UC3);
	final_class_mat =   cvCreateMat(src->height,src->width,CV_8UC3);
	mask			=	cvCloneMat(temp_mat);
	temp_mat_sub	=	cvCreateMat(src->height,src->width,CV_32FC3);
	temp_mat_sub2	=	cvCreateMat(src->height,src->width,CV_32FC3);
	cvZero(temp);
	cvCvtColor(src,temp,CV_RGB2Luv);					//RGB to Luv空间
	distance		=	cvCreateMat(src->height,src->width,CV_32FC1);
	result			=	cvCreateMat(src->height,src->width,CV_8UC1);
	cvConvert(temp,temp_mat);							//IplImage to Mat
	cn	=	cvCreateMat(src->height,src->width,CV_32FC1);
	cn1	=	cvCloneMat(cn);
	cn2	=	cvCloneMat(cn);
	cn3	=	cvCloneMat(cn);
	storage = cvCreateMemStorage(0);
//-------------------------------------------计算搜索窗口半径 r --------------------------------------------
	if(r!=NULL)
		r1=r;
	else
	{
		cvscalar1	=	cvSum(temp_mat);
		avg_img[0]	=	cvscalar1.val[0]/(src->width * src->height);
		avg_img[1]	=	cvscalar1.val[1]/(src->width * src->height);
		avg_img[2]	=	cvscalar1.val[2]/(src->width * src->height);
		cvscalar1	=	cvScalar(avg_img[0],avg_img[1],avg_img[2],NULL);
		cvScale(temp_mat,temp_mat_sub,1.0,0.0);
		cvSubS(temp_mat_sub , cvscalar1 , temp_mat_sub ,NULL);
		cvMul(temp_mat_sub , temp_mat_sub , temp_mat_sub2);
		cvscalar1	=	cvSum(temp_mat_sub2);
		r1			=	0.4*cvSqrt( (cvscalar1.val[0] + cvscalar1.val[1] + cvscalar1.val[2])/(src->width * src->height));;
	}
	//初始化随机数生成种子
	srand((unsigned)time(NULL));
	
//--------------------循环，使用meanshift进行特征空间分析，终止条件是Nmin--------------------------------------
	do
	{
//--------------------------------------------初始化搜索窗口位置-------------------------------------------
		run_meanshift_slec_number++;
		cvSet(distance,cvScalar(r1*r1,NULL,NULL,NULL),NULL);
		for( i = 0 ; i < 25 ; i++)
		{
			meanpoint[i].pt.x = rand()%src->width;
			meanpoint[i].pt.y = rand()%src->height;
		}
		cvScale(temp_mat,temp_mat_sub,1.0,0.0);
		for( i = 0 ; i < 25 ; i++)
		{
			/*cvSubS(temp_mat_sub ,cvScalar(cvGetReal3D(temp_mat,meanpoint[i].pt.x,meanpoint[i].pt.y,0),
				cvGetReal3D(temp_mat,meanpoint[i].pt.x,meanpoint[i].pt.y,1),
				cvGetReal3D(temp_mat,meanpoint[i].pt.x,meanpoint[i].pt.y,2),
				NULL),temp_mat_sub,NULL);*/
			cvSplit(temp_mat_sub,cn,cn1,cn2,NULL);
			cvSubS(temp_mat_sub,cvScalar(cvmGet(cn,meanpoint[i].pt.y,meanpoint[i].pt.x),
				cvmGet(cn1,meanpoint[i].pt.y,meanpoint[i].pt.x),
				cvmGet(cn2,meanpoint[i].pt.y,meanpoint[i].pt.x),NULL),temp_mat_sub,NULL);
			cvMul(temp_mat_sub,temp_mat_sub,temp_mat_sub2,1);
			cvSplit(temp_mat_sub2,cn,cn1,cn2,NULL);
			cvAdd(cn,cn1,cn3,NULL);
			cvAdd(cn2,cn3,cn3,NULL);			//cn3中存放着，当前随机点与空间中其它点距离的平方。
			cvCmp(cn3,distance,result,CV_CMP_LE);		//距离小于搜索半径则result相应位为1
			cvAndS(result,cvScalar(1,NULL,NULL,NULL),result,NULL);
			cvscalar1	=	cvSum(result);
			meanpoint[i].con_f_number = (int)cvscalar1.val[0];
		}
		for(i = 0 ; i < 25 ; i++)
		{
			order[i]=i;
		}
		for(i = 0 ; i < 25 ; i++)
			for(j = 0 ; j < 25-i-1; j++)
			{
				if(meanpoint[order[j]].con_f_number < meanpoint[order[j+1]].con_f_number)
				{
					temp_number=order[j];
					order[j]=order[j+1];
					order[j+1]=temp_number;
				}
			}
//--------------------------------------------meanshift算法------------------------------------------------	
		double	temp_mean[3];

		for( i = 0 ; i < 25 ; i++)
		{
			cvScale(temp_mat,temp_mat_sub,1.0,0.0);
			cvSplit(temp_mat_sub,cn,cn1,cn2,NULL);
			temp_mean[0]	=	cvmGet(cn  , meanpoint[order[i]].pt.y , meanpoint[order[i]].pt.x);
			temp_mean[1]	=	cvmGet(cn1 , meanpoint[order[j]].pt.y , meanpoint[order[i]].pt.x);
			temp_mean[2]	=	cvmGet(cn2 , meanpoint[order[j]].pt.y , meanpoint[order[i]].pt.x);

			//meanshift过程
			do
			{
				//计算出在搜索窗口内的特征点，并且生成对应的模板，即对应的点置一的矩阵表示对应的点在搜索框内
				cvScale(temp_mat,temp_mat_sub,1.0,0.0);
				cvSubS(temp_mat_sub,cvScalar(temp_mean[0],temp_mean[1],temp_mean[2],NULL),temp_mat_sub,NULL);
				cvMul(temp_mat_sub,temp_mat_sub,temp_mat_sub2,1);
				cvSplit(temp_mat_sub2 , cn , cn1 , cn2 , NULL );
				cvAdd(cn,cn1,cn3,NULL);
				cvAdd(cn2,cn3,cn3,NULL);			//cn3中存放着，当前随机点与空间中其它点距离的平方。
				cvCmp(cn3,distance,result,CV_CMP_LE);		//距离小于搜索半径则result相应位为0XFF
				
				
				//计算shiftor
				cvCopy(temp_mat , final_class_mat ,NULL);				//
				cvMerge(result , result ,result ,NULL,mask);
				cvAnd(final_class_mat , mask ,final_class_mat ,NULL);	//与mask（3通道,0XFF）做与操作，把搜索半径外的点置零
				cvScale(final_class_mat,temp_mat_sub,1.0,0.0);			//搜索半径内的点从8U转换成32F

				cvAndS(result,cvScalar(1,NULL,NULL,NULL),result,NULL);		//相应位set 1
				cvscalar1	=	cvSum(result);				//reslut 作为 模板 ,返回搜索窗口内的特征数

				cvSubS(temp_mat_sub,cvScalar(temp_mean[0],temp_mean[1],temp_mean[2],NULL),temp_mat_sub,result);
				cvscalar2	=	cvSum(temp_mat_sub);
				cvscalar2.val[0] = cvscalar2.val[0]/cvscalar1.val[0] ;
				cvscalar2.val[1] = cvscalar2.val[1]/cvscalar1.val[0] ;
				cvscalar2.val[2] = cvscalar2.val[2]/cvscalar1.val[0] ;
				shiftor		=	cvSqrt(pow(cvscalar2.val[0], 2) + pow(cvscalar2.val[1], 2) +	pow(cvscalar2.val[2], 2));
				temp_mean[0]=temp_mean[0]+cvscalar2.val[0];
				temp_mean[1]=temp_mean[1]+cvscalar2.val[1];
				temp_mean[2]=temp_mean[2]+cvscalar2.val[2];
				/*cvCopy(temp_mat , final_class_mat ,NULL);	//
				cvMerge(result , result ,result ,NULL,mask);
				cvAnd(final_class_mat , mask ,final_class_mat ,NULL);	//与result做与操作，把搜索半径外的点置零
				cvScale(final_class_mat,temp_mat_sub,1.0,0.0);			//搜索半径内的点从8U转换成32F
				cvSplit(temp_mat_sub,cn,cn1,cn2,NULL);
				cvSubS(cn , cvScalar(temp_mean[0],NULL,NULL,NULL),cn,result);
				cvSubS(cn1, cvScalar(temp_mean[1],NULL,NULL,NULL),cn1,result);
				cvSubS(cn2, cvScalar(temp_mean[2],NULL,NULL,NULL),cn2,result);
				cvMerge(cn,cn1,cn2,NULL,temp_mat_sub);
				cvscalar2	=	cvSum(temp_mat_sub);
				shiftor		=	cvSqrt(pow(cvscalar2.val[0] , 2) + pow(cvscalar2.val[1] , 2) +	pow(cvscalar2.val[2] , 2));
				temp_mean[0]=temp_mean[0]+cvscalar2.val[0];
				temp_mean[1]=temp_mean[1]+cvscalar2.val[1];
				temp_mean[2]=temp_mean[2]+cvscalar2.val[2];*/
			}
			while(shiftor>0.1);	//meanshift算法过程
//--------------------------------------------去除不重要特征-----------------------------------------------
			if(k==0)
			{
				feature[k].pt.x = temp_mean[0];
				feature[k].pt.y = temp_mean[1];
				feature[k].pt.z = temp_mean[2];
				feature[k].number= (int)cvscalar1.val[0];	//因为小于等于的情况成立时，result对应位置是0XFF，不成立时对应位置为0
				pNmin	= (int)cvscalar1.val[0];				//此特征搜索窗口内，特征空间的向量个数
				feature[k].result=cvCreateMat(src->height,src->width,CV_8UC1);
				cvAndS(result,cvScalar(1,NULL,NULL,NULL),result,NULL);
				cvCopy(result,feature[k].result,NULL);
				k++;
			}
			else
			{
				int flag = 0;
				for(j = 0 ; j < k ; j++)
				{
					if(pow(temp_mean[0]-feature[j].pt.x , 2) + pow(temp_mean[1]-feature[j].pt.y ,2) + pow(temp_mean[2]-feature[j].pt.z, 2)
						< r1*r1)
					{
						flag = 1;
						break;
					}
				}
				if(flag==0)
				{
					feature[k].pt.x = temp_mean[0];
					feature[k].pt.y = temp_mean[1];
					feature[k].pt.z = temp_mean[2];
					feature[k].number=(int)cvscalar1.val[0];
					pNmin	= (int)cvscalar1.val[0];				//此特征搜索窗口内，特征空间的向量个数
					feature[k].result=cvCreateMat(src->height,src->width,CV_8UC1);
					cvCopy(result,feature[k].result,NULL);
					k++;
					//if(pNmin < Nmin )
					//	break;
				}
			}//去除不重要特征
			//if(pNmin < Nmin)
			//	break;
		}	//

	}while(pNmin > Nmin || run_meanshift_slec_number>60 );

	//------------------------------------------------后处理---------------------------------------------------------
	cvSetZero(result);
	for( i = 0 ; i < k ; i ++)
	{
		cvOr(result,feature[i].result,result,NULL);
	}

	cvScale(temp_mat,temp_mat_sub,1.0,0.0);
	cvSplit(temp_mat_sub,cn,cn1,cn2,NULL);

	for(i = 0 ; i < src->width ; i++)
		for( j = 0 ; j < src->height ; j++)
		{
			if(cvGetReal2D(result,j,i)==0)		//未分类的像素点，进行分类，为最近的特征中心
			{
				double unclass_dis , min_dis;
				int min_dis_index;
				for( p = 0 ; p < k ; p++ )
				{
					unclass_dis = pow(feature[p].pt.x - cvmGet(cn,j,i),2)	//(temp_mat,i,j,0) ,2)
						+ pow(feature[p].pt.y - cvmGet(cn1,j,i),2) //(temp_mat,i,j,1) ,2)
						+ pow(feature[p].pt.z - cvmGet(cn2,j,i),2);//(temp_mat,i,j,2) ,2);
					if(p==0)
					{
						min_dis = unclass_dis;
						min_dis_index = p;
					}
					else
					{
						if(unclass_dis < min_dis)
						{
							min_dis = unclass_dis;
							min_dis_index = p;
						}
					}
				}// end for 与特征比较
				cvSetReal2D(feature[min_dis_index].result ,j  ,i ,1);
			}
		}//完成未分类的像素点的分类
	cvSetZero(final_class_mat);
	for( i = 0 ; i < k ; i++)
	{
		cvSet(temp_mat, cvScalar(rand()%255,rand()%255,rand()%255,rand()%255), feature[i].result);
		cvCopy(temp_mat,final_class_mat,feature[i].result);
	}
	cvConvert(final_class_mat,dst);
	//删除小于Ncon大小的区域
	for( i = 0 ; i < k ; i++)
	{
		cvClearMemStorage(storage);
		if(seq) cvClearSeq(seq);
		cvConvert( feature[i].result , gray);
		cvFindContours( gray , storage , & seq ,sizeof(CvContour) , CV_RETR_LIST);
		for(temp_seq = seq ; temp_seq ; temp_seq = temp_seq->h_next)
		{
			CvContour * cnt = (CvContour*)seq;
			if(cnt->rect.width * cnt->rect.height < Ncon)
			{
				prev_seq = temp_seq->h_prev;
				if(prev_seq)
				{
					prev_seq->h_next = temp_seq->h_next;
					if(temp_seq->h_next) temp_seq->h_next->h_prev = prev_seq ;
				}
				else
				{
					seq = temp_seq->h_next ;
					if(temp_seq->h_next ) temp_seq->h_next->h_prev = NULL ;
				}
			}
		}//
		cvDrawContours(src, seq , CV_RGB(0,0,255) ,CV_RGB(0,0,255),1);
	}

	//----------------释放空间-------------------------------------------------------	
	cvReleaseImage(& temp);
	cvReleaseImage(& gray);
	cvReleaseMat(&distance);
	cvReleaseMat(&result);
	cvReleaseMat(&temp_mat);
	cvReleaseMat(&temp_mat_sub);
	cvReleaseMat(&temp_mat_sub2);
	cvReleaseMat(&final_class_mat);
	cvReleaseMat(&cn);
	cvReleaseMat(&cn1);
	cvReleaseMat(&cn2);
	cvReleaseMat(&cn3);
}