histogram.c

学习跟踪的好程序

  /*****************************************************************************
 * File:	Histogram.c
 * Desc:	functions of Histogram Tracking Algorithm
 * Author:	Xuhui Zhou  @ Carnegie Mellon University
 * Date:	12/09/2003
 *****************************************************************************/

#include "Histogram.h"

//histogram bins
int  RBINS = 10;
int  GBINS = 10;
int  BBINS = 10;

//global histogram
CvHistogram  *gHistInit; //init click hist
CvHistogram  *gHistModel; //model hist without outer ring
CvHistogram	 *gHistRatio;	//ratio of foreground/(foreground+background)
CvHistogram	 *gHistDivide;	//ratio of foreground/(foreground+background)
CvHistogram	 *gHistDiff;	//Difference of foreground-background

//projection image
IplImage	*gImgBackProj;
IplImage	*gImgFgMask;
IplImage	*gImgObjMask;

//global path
char		gExePath[_MAX_PATH];

POINT		gGravityCenter;
int			gBorder = 20; //border width for outer ring box
int			gTrackMethod = 1; //SHIFT, 2:CONV

int			gMaskPixCount;
float		gMaskFgRatio;

///////////////////////////////////////////////////////////////////////////////
void his_PrintHistogram(CvHistogram *hist, char* filename)
{
	FILE	*file;
	float	data;
	int 	h, s, v;

	file = fopen(filename, "w");
	for( h = 0; h < RBINS; h++ ){			
		for( s = 0; s < GBINS; s++ ){			
			fprintf(file, "%02d_%02d_(0~%02d): ", h, s, BBINS-1);
			for(v = 0; v < BBINS; v++ ){						
				data = cvQueryHistValue_3D(hist,h,s,v);
				fprintf(file, "%3d ", (long)data );						
			}			
			fprintf(file, "\n");
		}		
	}
	fclose(file);	
}

/******************************************************************************/
CvHistogram* his_CalcHistogram(IplImage* inImg){
/******************************************************************************/
	CvHistogram *outHist;
	CvSize imgSize = cvGetSize(inImg);	
	IplImage* h_plane = cvCreateImage(imgSize, 8, 1 );
    IplImage* s_plane = cvCreateImage(imgSize, 8, 1 );
    IplImage* v_plane = cvCreateImage(imgSize, 8, 1 );        
    float h_ranges[] = { 0, 255}; /* hue varies from 0 (~0皉ed) to 180 (~360皉ed again) */
    float s_ranges[] = { 0, 255}; /* saturation varies from 0 (black-gray-white) to 255 (pure spectrum color) */
	float v_ranges[] = { 0, 255}; 
    int		hist_size[] = {RBINS, GBINS, BBINS};
	float*	ranges[] = { h_ranges, s_ranges, v_ranges};
	IplImage* planes[] = { h_plane, s_plane, v_plane};
	float max_value = 0;
    
    cvCvtPixToPlane(inImg, h_plane, s_plane, v_plane, 0 );
    outHist = cvCreateHist(3, hist_size, CV_HIST_ARRAY, ranges, 1 );
    cvCalcHist( planes, outHist, 0, 0 );

	//nornalize histogram
	cvNormalizeHist(outHist, 255);
	cvGetMinMaxHistValue( outHist, 0, &max_value, 0, 0 );

	//cvThreshHist(outHist, max_value*0.5);
	//cvGetMinMaxHistValue( outHist, 0, &max_value, 0, 0 );	

	cvReleaseImage(&h_plane);
	cvReleaseImage(&s_plane);	
	cvReleaseImage(&v_plane);	
	return outHist;
}


/******************************************************************************/
CvHistogram* his_HistogramFromArray(float* modelhist){
/******************************************************************************/		
	CvHistogram *outHist;
	int  d[1] = {256};
	int  bin;
	float *p;		
	float range[] = {0,255};
	float *ranges[] = {range};	
	float maxVal,factor; //minVal
	//int		minIdx, maxIdx;

	outHist = cvCreateHist(1,d,CV_HIST_ARRAY,ranges,1);
	
	maxVal = 0;
	for (bin=0; bin<256; bin++){		
		p = cvGetHistValue_1D(outHist, bin);
		//*p = modelhist[bin];
		*p = bin+1.0f;
		if (modelhist[bin]>maxVal){
			maxVal = modelhist[bin];
		}
	}
		
	
	if (maxVal!=0){
		factor = 255;
		cvNormalizeHist(outHist, factor);
 	}	

	return outHist;
}

/******************************************************************************/
CvHistogram* his_HistogramRatio(CvHistogram* hist1_model, CvHistogram *hist2_denom)
//ratio by probDensity
{
	CvHistogram *outHist=NULL;
	
	cvCopyHist(hist1_model, &outHist);
	cvClearHist(outHist);
	cvCalcProbDensity(hist2_denom,hist1_model,outHist,255);
		
	return outHist;
}

/******************************************************************************/
CvHistogram* his_HistogramDivide(CvHistogram* hist1, CvHistogram *hist2)
//divide and normalize
{
	CvHistogram *histRst=NULL;	
	int b1, b2, b3;
	float   *p1, *p2, *p3;
	float	rst;

	cvCopyHist(hist1, &histRst);
	cvClearHist(histRst);

	for( b1 = 0; b1 < RBINS; b1++ ){			
		for( b2 = 0; b2 < GBINS; b2++ ){
			for( b3 = 0; b3 < BBINS; b3++ ){			
				p1 = cvGetHistValue_3D(hist1,b1,b2,b3);
				p2 = cvGetHistValue_3D(hist2,b1,b2,b3);
				p3 = cvGetHistValue_3D(histRst,b1,b2,b3);
				
				if ((*p2)==0) rst = 0;
				else
					rst = (*p1)/(*p2);
				(*p3) = max(rst, 0);				
			}
		}
	}					

	his_HistogramNormalize(histRst);
		
	return histRst;
}

/******************************************************************************/
void his_HistogramDiff(CvHistogram* hist1, CvHistogram *hist2, CvHistogram *histRst)
//inner minus outer; (foreground) - (foreground+background)
{	
	int b1, b2, b3;
	float   *p1, *p2, *p3;
	float	rst;

	for( b1 = 0; b1 < RBINS; b1++ ){			
		for( b2 = 0; b2 < GBINS; b2++ ){
			for( b3 = 0; b3 < BBINS; b3++ ){			
				p1 = cvGetHistValue_3D(hist1,b1,b2,b3);
				p2 = cvGetHistValue_3D(hist2,b1,b2,b3);
				p3 = cvGetHistValue_3D(histRst,b1,b2,b3);
				
				rst = (*p1)*2 - (*p2);
				(*p3) = max(rst, 0);				
			}
		}
	}					
}

/******************************************************************************/
void his_HistogramNormalize(CvHistogram* inHist)
//make the max as 255, multiply by 255/max
{
	int b1, b2, b3;
	float	maxVal, minVal;
	//float	val; //histogram count value
	float   *pix;

	//get max value
	cvGetMinMaxHistValue( inHist, &minVal, &maxVal, 0, 0 );			

	if (maxVal==0) return;

	//get min non-zero value
	minVal = 255;
	for( b1 = 0; b1 < RBINS; b1++ ){			
		for( b2 = 0; b2 < GBINS; b2++ ){
			for( b3 = 0; b3 < BBINS; b3++ ){			
				pix = cvGetHistValue_3D(inHist,b1,b2,b3);
				if ((*pix)>0 && (*pix)<minVal){
					minVal = (*pix);
				}
			}
		}
	}					


	//normalize
	for( b1 = 0; b1 < RBINS; b1++ ){			
		for( b2 = 0; b2 < GBINS; b2++ ){
			for( b3 = 0; b3 < BBINS; b3++ ){			
				pix = cvGetHistValue_3D(inHist,b1,b2,b3);
				if (*pix>0){
					(*pix) = ((*pix)-minVal)*255/(maxVal-minVal);
				}
			}
		}
	}					
}

/******************************************************************************/
void his_HistogramThresholdUp(CvHistogram* inHist, float thresholdVal)
//make pix 255 if above threshold
{
	int b1, b2, b3;
	float	max_value = 0;
	//float	val; //histogram count value
	float   *pix;

	//get max value
	cvGetMinMaxHistValue( gHistModel, 0, &max_value, 0, 0 );			

	//normalize
	for( b1 = 0; b1 < RBINS; b1++ ){			
		for( b2 = 0; b2 < GBINS; b2++ ){
			for( b3 = 0; b3 < BBINS; b3++ ){			
				pix = cvGetHistValue_3D(inHist,b1,b2,b3);
				if (*pix>thresholdVal){
					(*pix) = 255;			
				}
			}
		}
	}					
}

/******************************************************************************/
CvHistogram* his_HistogramAdd(CvHistogram* hist1, CvHistogram *hist2, float weight){
/******************************************************************************/
	CvHistogram *outHist=NULL;	
	int	h,s,v;
	float h1, h2,*p;

	cvCopyHist(hist1, &outHist);

	for( h = 0; h < RBINS; h++ ){			
		for( s = 0; s < GBINS; s++ ){
			for( v = 0; v < BBINS; v++ ){			
				h1 = cvQueryHistValue_3D(hist1,h,s,v);
				h2 = cvQueryHistValue_3D(hist2,h,s,v);
				p = cvGetHistValue_3D(outHist,h,s,v);
				(*p) = h1*(1-weight) + h2*weight;
			}
		}
	}			
		
	return outHist;
}
/******************************************************************************/
void his_BackProject(BYTE *imgin, BYTE *imgout, CvHistogram *hist){
/******************************************************************************/
	int		imgWidth, imgHeight;	
	IplImage* imgSrc[CV_MAX_DIM];
	IplImage* imgRst = cvCreateImage( cvSize(320,243), IPL_DEPTH_8U, 1 );
	IplImage* tmp	=	cvCreateImage( cvSize(320,243), IPL_DEPTH_8U, 1 );

	imgSrc[0] = cvCreateImage( cvSize(320,243), IPL_DEPTH_8U, 1 );
	//copy input image
	imgWidth = 320;
	imgHeight = 243;	
	memcpy(imgSrc[0]->imageData, imgin, imgWidth*imgHeight);	

	//calculate backproject
	cvCalcBackProject(imgSrc, imgRst, hist);

//	//debug display		
//	cvSaveImage("result/imgsrc.bmp", imgSrc[0]);
//	cvSaveImage("result/imgRst.bmp", imgRst);
//	img_Output_C1_float(imgRst->imageData, imgWidth, imgWidth, imgHeight, "result/imgRst.txt");	

	//copy to output
	memcpy(imgout,imgRst->imageData,imgWidth*imgHeight);
	
	//clear memory
	cvReleaseImage( &imgSrc[0]);
	cvReleaseImage( &imgRst);
}

/******************************************************************************/
IplImage* his_CalcBackProject(IplImage *inImgHSV, CvHistogram *inHist){
/******************************************************************************/		
	//calculate back project		
	CvSize imgSize = cvGetSize(inImgHSV);	
	IplImage*	h_plane = cvCreateImage(imgSize, 8, 1 );
	IplImage*	s_plane = cvCreateImage(imgSize, 8, 1 );
	IplImage*	v_plane = cvCreateImage(imgSize, 8, 1 );
	IplImage*	planes[] = { h_plane, s_plane, v_plane };    	
	IplImage*	bpRst = cvCreateImage(imgSize, 8, 1 );	
		
	//seperate HSV planes
	cvCvtPixToPlane(inImgHSV, h_plane, s_plane, v_plane, 0 );
    
	//calculate backproject
	cvCalcBackProject(planes, bpRst, inHist);

	//cvThreshold(outImgBP, outImgBP, 50, 255, CV_THRESH_TOZERO );
	
	//debug display	
//	{
//		IplImage* tmp = cvCreateImage(imgSize, 8, 1 );
//		cvFlip(inImgHSV,NULL,0);
//		cvSaveImage("result/inImgHSV.bmp",inImgHSV);
//		cvFlip(outImgBP,tmp,0);
//		cvSaveImage("result/outImgBP.bmp",tmp);		
//	}
	
	//clear memory
	cvReleaseImage( &h_plane);
	cvReleaseImage( &s_plane);
	cvReleaseImage( &v_plane);
	return bpRst;
}

/******************************************************************************/
void his_CalcForeground2(IplImage *inImgHSV, CvHistogram *inHist, CvRect inRect, IplImage *outImgFg){
/******************************************************************************/		
	int			border;
	CvRect		rectSearch;	 //search window in tracking
	IplImage	*imgsrc = cvCloneImage(inImgHSV);		
	IplImage	*imgrst;		
	double		dRstMax;
	int			row, col;
	double		pixval;
	
	//extent neighboring search window
	border = gBorder;
	rectSearch.x = inRect.x - border;
	rectSearch.y = inRect.y - border;
	rectSearch.width = inRect.width+2*border;
	rectSearch.height = inRect.height+2*border;
	his_CheckBoundRect(&rectSearch, inImgHSV->width, inImgHSV->height);		
	
	//get back project		
	imgrst = his_CalcBackProject(imgsrc, inHist);

	//scale to 0~255
	cvMinMaxLoc(imgrst,NULL, &dRstMax, NULL,NULL,NULL);	
	if (dRstMax!=0){		
		cvConvertScale(imgrst, imgrst, 255/dRstMax,0);
	}
	
	//copy imgrst to outImgFg
	cvSetZero(outImgFg);
	for (col=rectSearch.x; col<rectSearch.x+rectSearch.width;col++){
		for (row=rectSearch.y;row<rectSearch.y+rectSearch.height;row++){
			pixval = cvGetReal2D(imgrst, row, col);
			cvSetReal2D(outImgFg, row, col, pixval);
		}
	}
	
	//release memory
	cvReleaseImage( &imgsrc);
	cvReleaseImage( &imgrst);
}

/******************************************************************************/
void his_CheckBoundRect(CvRect *rect, int imgWidth, int imgHeight){
/******************************************************************************/
	if (rect->x <0) rect->x=0;
	if (rect->x >imgWidth) rect->x = imgWidth;
	if (rect->y <0) rect->y=0;
	if (rect->y >imgHeight) rect->y = imgHeight;

	if (rect->width >=imgWidth - rect->x)
		rect->width = imgWidth -1 - rect->x;
	if (rect->height>=imgHeight - rect->y)
		rect->height = imgHeight -1 - rect->y;
}

///////////////////////////////////////////////////////////////////////////////
IplImage* his_GetBackProjImage()
{
	return gImgBackProj;
}

///////////////////////////////////////////////////////////////////////////////
void his_ReadHistBinsFromFile(){
	FILE *filein;
	char *binfile = "histogram.txt";
	char prompt[10];
	//int i;
	//char buffer[_MAX_PATH];
	
	//_getcwd(buffer, _MAX_PATH);

	filein = fopen(binfile, "r");
	if (filein==NULL){
		return;
	}
		
	fscanf(filein, "%s", &prompt);
	fscanf(filein, "%d", &RBINS);
	fscanf(filein, "%s", &prompt);
	fscanf(filein, "%d", &GBINS);
	fscanf(filein, "%s", &prompt);
	fscanf(filein, "%d", &BBINS);
	
	fclose(filein);

	//in case of wrong setting
	if(RBINS<=0) 
		RBINS = 10;