colormodelycbcr.cpp

color based segment in YCbCr space

#ifdef _CH_
#define WIN32
#error "The file needs cvaux, which is not wrapped yet. Sorry"
#endif

#ifndef _EiC
#include "cv.h"
#include "highgui.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <math.h>
#include <float.h>
#include <limits.h>
#include <time.h>
#include <ctype.h>
#include <cvaux.h>

#include <fstream.h>
#include <limits>
#endif

#ifdef _EiC
#define WIN32
#endif

#ifdef _WIN32
/* use clock() function insted of time() */
#define time( arg ) (((double) clock()) / CLOCKS_PER_SEC)
#endif /* _WIN32 */

#ifndef PATH_MAX
#define PATH_MAX 512
#endif /* PATH_MAX */

void main( int argc, char** argv )
{
	int i,j;
//	int s = 0;    //计算总的象素点数
//	int n = 0;
//
//	FILE* info;
//	char* infoname;
//    char fullname[PATH_MAX];
//    //char detfilename[PATH_MAX];
//    char* filename;
//    //char detname[] = "";
//
//	int num;
//
	CvMat *M_Cov = cvCreateMat( 2, 2, CV_64FC1 );
	CvMat *M_U = cvCreateMat( 2, 1, CV_64FC1 );
	CvMat *T_Cov = cvCreateMat( 2, 2, CV_64FC1 );//协方差的逆
//    
//	for( i = 1; i < argc; i++ )
//    {
//        if( !strcmp( argv[i], "-info" ) )
//        {
//            infoname = argv[++i];
//        }
//        else if( !strcmp( argv[i], "-num" ) )
//        {
//            num = (int) atof( argv[++i] );
//        }
//	}
//
//	strcpy( fullname, infoname );
//	filename = strrchr( fullname, '\\' );
//	if( filename == NULL )
//	{
//		filename = strrchr( fullname, '/' );
//	}
//	if( filename == NULL )
//	{
//		filename = fullname;
//	}
//	else
//	{
//		filename++;
//	}
//
//	info = fopen( infoname, "r" );
//	if ( info != NULL)
//	{
//		
//		//unsigned char R8,B8,G8; //0~255
//		//float R,B,G, r,g,b;
//		//int rr,bb,gg;
//        unsigned char Y,Cb,Cr;
//
	    double U[2]={0};//均值
		

		double Cov[2][2]={0};//协方差
//	
//		float Cb_Cr[256][256]={0};
//		float p_CbCr[256][256]={0};
//		float cb[256]={0};
//		float cr[256]={0};
//		float p_cb[256]={0};
//		float p_cr[256]={0};

//		IplImage* image = NULL;
//		IplImage* Y_image = NULL;
//
//		while ( !feof( info ))
//		{
//			 if ( fscanf( info, "%s", filename) != 1 )
//				 break;
//			 image = cvLoadImage( fullname );//从文件中读入一幅图像
//			 if( !image ) 
//				 continue;
//
//			 s = s+image->height*image->width;
//			 n += 1;
//				 
//			 //把图像由RGB空间转换到YCbCr空间中
//			 IplImage* Y_image = cvCreateImage( cvGetSize(image),IPL_DEPTH_8U, 3 );//若为IPL_DEPTH_32F则报错
//			 cvCvtColor(image, Y_image, CV_RGB2YCrCb);
//
//			 //读取CbCr的值
//			 for ( i = 0; i < Y_image->height; i ++ )
//			 {
//				 for ( j = 0; j < Y_image->width; j ++ )
//				 {
//					 Y = (uchar) (*( Y_image->imageData + Y_image->widthStep*i+j*3 ));
//					 Cb = (uchar) (*( Y_image->imageData + Y_image->widthStep*i+j*3+1 ));
//				     Cr = (uchar) (*( Y_image->imageData + Y_image->widthStep*i+j*3+2 ));
//
//					 //二维直方图(b,g)
//                     Cb_Cr[Cb][Cr]+=1;
//					 cb[Cb]+=1;
//					 cr[Cr]+=1;
//				 }
//			 }
//		}        
//		//每一级的概率
//        for( i=0;i<256;i++ )
//		{
//			for(j=0;j<256;j++)
//			{
//				p_CbCr[i][j]=Cb_Cr[i][j]/s;
//			}
//		}
//	    for(i=0;i<256;i++)
//		{
//			p_cb[i]=cb[i]/s;
//		    p_cr[i]=cr[i]/s;
//		}           
//
//		//二维统计值方图的均值
//		for(i=0;i<256;i++)
//		{
//			U[0]=U[0]+i*p_cb[i];
//            U[1]=U[1]+i*p_cr[i];
//
//        //二维统计值方图的协方差
//			Cov[0][0] += i*i*p_cb[i];
//			Cov[1][1] += i*i*p_cr[i];
//		}
//
//        Cov[0][0]=Cov[0][0]-U[0]*U[0];
//		Cov[1][1]=Cov[1][1]-U[1]*U[1];
//		for(i=0;i<256;i++)
//		{
//			for(j=0;j<256;j++)
//			{
//				Cov[0][1] += i*j*p_CbCr[i][j];
//			}
//		}
//		Cov[0][1]=Cov[0][1]-U[0]*U[1];
//		Cov[1][0]=Cov[0][1];
//		
		
        Cov[0][0]=55.3961;
		Cov[0][1]=-36.4732;
		Cov[1][0]=-36.4732;
		Cov[1][1]=60.4410;

		U[0]=112.0816;
		U[1]=148.8160;
        cvInitMatHeader( M_Cov, 2, 2, CV_64FC1, Cov);//数组变为矩阵
		cvInitMatHeader( M_U, 2, 1, CV_64FC1, U);//数组变为矩阵
	    cvInvert( M_Cov, T_Cov, CV_LU);	//用高斯消去法求逆矩阵,并存在T_Cov中
	
	
    IplImage* Test_image = NULL;
	IplImage* Test_image1 = NULL;
   
    double X[2];
	CvMat *M_X = cvCreateMat( 2, 1, CV_64FC1 );//Creates new matrix,the pointer to the matrix is M_X
	CvMat *Sub_XU = cvCreateMat( 2, 1, CV_64FC1 );
	CvMat *T_Sub_XU = cvCreateMat( 1, 2, CV_64FC1 );
	CvMat *temp = cvCreateMat( 1, 2, CV_64FC1 );
	CvMat *D = cvCreateMat( 1, 1, CV_64FC1 );
    //CString img=2.jpg;
	Test_image = Test_image1 = cvLoadImage("39.jpg" );//从文件中读入一幅图像
   

	IplImage* Test_Y_image = cvCreateImage( cvGetSize(Test_image),IPL_DEPTH_8U, 3 );//Creates header and allocates data
	IplImage* B_image = cvCreateImage( cvGetSize(Test_image),IPL_DEPTH_8U, 1 );

	cvCvtColor(Test_image, Test_Y_image, CV_RGB2YCrCb);//Converts image from RGB space to YCbCr


    
	for ( i = 0; i < Test_Y_image->height; i ++ )
	{
		for ( j = 0; j < Test_Y_image->width; j ++ )
		{
			X[0] = (uchar) (*( Test_Y_image->imageData + Test_Y_image->widthStep*i+j*3+1 ));//Cb  
			//char *imageData;  /* pointer to aligned image data */
			//int  widthStep;   /* size of aligned image row in bytes */

			X[1] = (uchar) (*( Test_Y_image->imageData + Test_Y_image->widthStep*i+j*3+2 ));//Cr
			
            cvInitMatHeader( M_X, 2, 1, CV_64FC1, X);//数组变为矩阵

            //求马氏距离
			cvSub(M_X,M_U,Sub_XU);//Computes per-element difference between two arrays
			cvTranspose( Sub_XU, T_Sub_XU );//转置
			cvMatMul( T_Sub_XU, T_Cov, temp);
            cvMatMul( temp, Sub_XU, D);

		    //double d = CV_MAT_ELEM(*D,double,0,0);
            if (CV_MAT_ELEM(*D,double,0,0)<=8.0)
				*( B_image->imageData + B_image->widthStep*i+j)=255;
			else
				*( B_image->imageData + B_image->widthStep*i+j)=0;
		}
	}
	cvNamedWindow( "binaryimage", CV_WINDOW_AUTOSIZE );
	cvShowImage( "binaryimage", B_image );
	cvWaitKey(10);//显示窗口并停留一段时间
	//cvvSaveImage( "binaryimage.jpg", B_image );

	//形态学滤波
	int size = 3;
	IplConvKernel* Element = cvCreateStructuringElementEx( size, size, 0, 0,CV_SHAPE_ELLIPSE);
	//cvErode( B_image, B_image, Element, 1 );
    //cvDilate( B_image, B_image, Element, 1 );
	cvMorphologyEx( B_image, B_image, NULL , Element, CV_MOP_OPEN, 2 );

	cvNamedWindow( "Morphology", CV_WINDOW_AUTOSIZE );
	cvShowImage( "Morphology", B_image );
	cvWaitKey(1000);//显示窗口并停留一段时间

	// Find ellipes in input_img by using its edge information
    // and draw the found ellipes on the images
    //void FindCircles(IplImage* input_img, IplImage* input_edge, CvBox2D* box2DArray, int* pEllipseNum)
    //{
	// DO SOME DECLARARTIONS
	// FOR ALL TYPES "Cv.." LOOK AT OPENCV-DOC !!
    //椭圆拟合
    CvMemStorage *storage;//Growing memory storage
    int header_size, count;
    CvSeq *contour;//Growable sequence of elements
	CvSeq MyContour;
    CvBox2D * myBox;
	CvBox2D* box2DArray;
    CvPoint *PointArray;
    CvPoint2D32f *PointArray32f;
	double MyContourArea;
	int pEllipseNum = 0;

	// ALLOCATE SOME MEMORY        
	myBox= (CvBox2D *) malloc(sizeof(CvBox2D));
	box2DArray= (CvBox2D *) malloc(sizeof(CvBox2D));
    header_size = sizeof(CvContour);
    storage = cvCreateMemStorage (0); // For FindContours.

    //THIS FUNCTION SEPARATES ALL CONTOURS
    // LOOK AT THE DOC !!!! AGAIN
    cvFindContours (B_image , storage, &contour, header_size, CV_RETR_LIST,  CV_CHAIN_APPROX_SIMPLE);
	
	// GO THROUGH ALL CONTOURS
    while(contour!=NULL)
    {
		
		if(CV_IS_SEQ_CURVE(contour))
		{
            // HOW MANY POITS HAS THE CONTOUR ??
            count = contour ->total;
           
            //ALLOC MEM
            PointArray = (CvPoint *)malloc(count * sizeof(CvPoint));
           
            //COPY THE POINTS TO A ARRAY
            cvCvtSeqToArray(contour, PointArray, CV_WHOLE_SEQ);

            // ALLOC MEM
            PointArray32f = (CvPoint2D32f *) malloc((count + 1) * sizeof(CvPoint2D32f));

            //CONVERT THE ARRAY TO A 2ND (32FLOAT)
            // THIS WE NEED FOR cvFitELLIPSE

            for (i=0; i<count-1; i++)
            {                              
                    PointArray32f[i].x = (float)(PointArray[i].x);
                    PointArray32f[i].y = (float)(PointArray[i].y);
            }
            PointArray32f[i].x = (float)(PointArray[0].x);
            PointArray32f[i].y = (float)(PointArray[0].y);        

            //cvFITELLIPSE NEEDS MIN. 6 POINTS
            if (count >= 7)
            {
				
				memcpy(&MyContour,contour,sizeof(CvSeq));
				MyContour.h_next=NULL;
				MyContour.h_prev=NULL;
				MyContourArea = fabs(cvContourArea(&MyContour,CV_WHOLE_SEQ));

				//FOR DEBUGGING YOU CAN DRAW THE CONTOUR
				//cvDrawContours(Test_image, contour, CV_RGB(255,0,0), CV_RGB(0,255,0),0);

				//cvNamedWindow( "DrawContours", CV_WINDOW_AUTOSIZE );
                //cvShowImage( "DrawContours", Test_image );
                //cvWaitKey(100);

				//FIND THE BEST FITTING ELLIPSE INTHIS CONTOUR
                cvFitEllipse(PointArray32f, count, myBox);

				//double sizeContour = fabs(cvContourArea(contour,CV_WHOLE_SEQ));
				double ellipseArea   = (double)(CV_PI*myBox->size.height*myBox->size.width/4);
                float lwb=0;//The ratio of 
                // IS THE ELLIPSE INSIDE THE IMAGE ??
				// Is it large enough?
                if ((myBox->center.x > 0) && (myBox->center.y >0) && (ellipseArea>100))// && (ellipseArea<20000))
                {  
					box2DArray[pEllipseNum].angle			= myBox->angle;
					box2DArray[pEllipseNum].center.x		= myBox->center.x;
					box2DArray[pEllipseNum].center.y		= myBox->center.y;
					box2DArray[pEllipseNum].size.height	= myBox->size.height;
					box2DArray[pEllipseNum].size.width		= myBox->size.width;
					box2DArray[pEllipseNum].angle			= myBox->angle;

					pEllipseNum++;		// one ellipse was found
                }
			}
			free(PointArray32f );
            free(PointArray );
		}
        // GOT TO THE NEXT CONTOUR (IF ANY)
        contour = contour->h_next;
    }      
	
	
    //free (myBox);            
    free (contour);
    cvReleaseMemStorage(&storage);

    for ( i=0; i<pEllipseNum; i++)
	{
		CvPoint center;
		center.x = (int) box2DArray[i].center.x;
		center.y = (int) box2DArray[i].center.y;
		CvSize size;
		size.width = (int)box2DArray[i].size.width/2;
		size.height = (int)box2DArray[i].size.height/2;
		//DRAW THE ELLIPSE
		cvEllipse(Test_image1, center, size, box2DArray[i].angle, 0, 360, CV_RGB(0,0,255),1);
	}
	cvNamedWindow( "Detect Result", CV_WINDOW_AUTOSIZE );
	cvShowImage("Detect Result", Test_image);
	cvvSaveImage( "39.jpg", Test_image1 );
	cvWaitKey(0);

}