OpenCV2马拉松第11圈——meanshift与直方图反向投影

收入囊中

  • meanshift图像聚类
  • meanshift object detect

葵花宝典

今天有点累,理论就讲少点吧T_T

meanshift中文是均值飘逸,就是给定一个点,然后会移动到概率密度最大的地方。

对于图像,什么是概率密度最大?

我们可以定义很多要素:

距离

RGB

HSV

下面我有个例子,就是用距离(x,y)和HSV(h,s,v)作图像聚类的。

于是我们有5个要素,当前点与其他点的距离,HSV越接近,则概率密度越高。

假定我们有一点(m,n),如何选择下一个点呢,如何在一个矩形中找到概率密度最高的那个点?

top,down,left,right是我们选的矩形,(m,n)是中心

for(s=top;s<=down;s++){
     				for(t=left;t<=right;t++)
     				{
      					ws=(s-m)*(s-m)+(t-n)*(t-n);//spatial information
      					ws/=(hs*hs);
      					ws=exp(-ws);

      					wr=(data[s*step+t*channels]-data[m*step+n*channels])*(data[s*step+t*channels]-data[m*step+n*channels]);
      					wr+=(data[s*step+t*channels+1]-data[m*step+n*channels+1])*(data[s*step+t*channels+1]-data[m*step+n*channels+1]);
      					wr+=(data[s*step+t*channels+2]-data[m*step+n*channels+2])*(data[s*step+t*channels+2]-data[m*step+n*channels+2]);

      					wr/=(hr*hr);

      					if(wr>1)
       						wr=0.;
      					else
       						wr=exp(-wr);

      					sumw+=wr*ws;
      					for(k=0;k<5;k++)//try
       						y[1][k]+=oridata[s*width+t][k]*wr*ws;
     				}
    			}

    			for(k=0;k<5;k++) //try
     				y[1][k]/=sumw;

    			//下一个要到的点
    			m=(int)(y[1][0]+0.5);
    			n=(int)(y[1][1]+0.5);

假设我们有了反向投影图,再用meanshift,不就能找到图像中最接近当前区域的位置了么~~

图像聚类

用meanshift可以实现图像聚类,因为本次重点是mean shift与back project的结合,所以关于聚类我就直接贴上代码了

算法超级easy,就是循环图像,对点作mean shift,每个点最大迭代100次,一直到点不动(收敛)为止,然后用vector记录下移动过程,将收敛点的颜色赋值给这些点,并记录哪些点已经赋值过了,下次就不用循环直接continue.

hs,hr是寻找的矩形大小,比EPSLON小就收敛。

note:是用x,y和HSV

//#include <stdafx.h>
#include <cxcore.h>
#include <cv.h>
#include <highgui.h>
#include <stdio.h>
#include <vector>
#include <utility>

using namespace cv;
using namespace std;

#define hs 25
#define hr 25
#define maxstep 100
#define EPSLON 0.01
#define width 1920
#define height 1080
float oridata[height*width][5];//try
int visited[height][width];

int main(int argc,char **argv)
{
	IplImage *oriImg,*luvImg,*fltImg,*afterImg; 

	char *filename = argv[1];
	oriImg=cvLoadImage(filename,1); 

	if(!oriImg){
		printf("cannot load the file.\n");
  		return -1;
 	}

 	luvImg=cvCreateImage(cvSize(width,height),oriImg->depth,oriImg->nChannels);
 	fltImg=cvCreateImage(cvSize(width,height),oriImg->depth,oriImg->nChannels);
 	afterImg=cvCreateImage(cvSize(width,height),oriImg->depth,oriImg->nChannels); 

 	cvCvtColor(oriImg,luvImg,CV_RGB2Luv); 

 	uchar *data,*newdata;
 	int channels, step,depth;

 	depth=luvImg->depth;
 	step=luvImg->widthStep;
 	channels=luvImg->nChannels;

 	data=(uchar *) luvImg->imageData;
 	newdata=(uchar *) fltImg->imageData;

 	int i,j,s,t,k,m,n,index;
 	int top,down,left,right;
 	float ws,wr;
 	float mhlength;
 	float sumw;
 	float y[2][5];//try, otherwise the second is 5

 	for(i=0;i<height;i++)
 	{
  		for(j=0;j<width;j++)
  		{
   			oridata[i*width+j][0]=i;
   			oridata[i*width+j][1]=j;
   			visited[i][j]=0;

    		oridata[i*width+j][2]=data[i*step+j*channels];
    		oridata[i*width+j][3]=data[i*step+j*channels+1];
    		oridata[i*width+j][4]=data[i*step+j*channels+2];
  		}
 	}

 	for(i=0;i<height;i++)
 	{
  		for(j=0;j<width;j++)
  		{
  			if(visited[i][j])
  				continue;

   			m=i;n=j;//当前的中心点
   			for(k=0;k<5;k++)//try
    			y[0][k]=oridata[i*width+j][k];//向量的初始值

			vector<pair<int,int> >vss;

   			for(index=0;index<maxstep;index++)//对当前的结点而言,最多迭代100次
   			{
   				pair<int,int>newone;
				newone = make_pair(m,n);
				vss.push_back(newone);

    			for(k=0;k<5;k++)//try
     				y[1][k]=0;

    			mhlength=0.;
    			sumw=0.;

    			top=m-hs;
    			down=m+hs;
    			left=n-hs;
    			right=n+hs;

    			if(top<0) top =0;
    			if(down>height-1) down=height-1;
    			if(left<0) left=0;
    			if(right>width-1) right=width-1;

    			for(s=top;s<=down;s++)
    			{
     				for(t=left;t<=right;t++)
     				{
      					ws=(s-m)*(s-m)+(t-n)*(t-n);//spatial information
      					ws/=(hs*hs);
      					ws=exp(-ws);
      					//ws=1-ws+(ws*ws)/2-(ws*ws*ws)/6+(ws*ws*ws*ws)/24-(ws*ws*ws*ws*ws)/120;

      					wr=(data[s*step+t*channels]-data[m*step+n*channels])*(data[s*step+t*channels]-data[m*step+n*channels]);
      					wr+=(data[s*step+t*channels+1]-data[m*step+n*channels+1])*(data[s*step+t*channels+1]-data[m*step+n*channels+1]);
      					wr+=(data[s*step+t*channels+2]-data[m*step+n*channels+2])*(data[s*step+t*channels+2]-data[m*step+n*channels+2]);

      					wr/=(hr*hr);

      					if(wr>1)
       						wr=0.;
      					else
       						wr=exp(-wr);

      					sumw+=wr*ws;
      					for(k=0;k<5;k++)//try
       						y[1][k]+=oridata[s*width+t][k]*wr*ws;
     				}
    			}

    			for(k=0;k<5;k++) //try
     				y[1][k]/=sumw;

    			//下一个要到的点
    			m=(int)(y[1][0]+0.5);
    			n=(int)(y[1][1]+0.5);

    			if(visited[m][n])
     				break;

    			if(m<hs||m>height-hs||n<hs||n>width-hs)
     				break;
    			else{
     				for(k=0;k<5;k++)//try
     				{
      					mhlength+=(y[1][k]-y[0][k])*(y[1][k]-y[0][k]);
      					y[0][k]=y[1][k];
     				}
     				mhlength=sqrt(mhlength);
     				if(mhlength<EPSLON)//找到极值点
      					break;
    			}
   			}//邻域处理结束

			for(int ii = 0;ii < vss.size();ii++)
			{
				int row = vss[ii].first;
				int line = vss[ii].second;
   				newdata[row*step+line*channels]=int(y[1][2]+0.5);//try
   				newdata[row*step+line*channels+1]=int(y[1][3]+0.5);
   				newdata[row*step+line*channels+2]=int(y[1][4]+0.5);
   				visited[row][line]=1;
   			}
  		}
 	}

 	cvCvtColor(fltImg,afterImg,CV_Luv2RGB);

 	cvNamedWindow("ori",1);
 	cvNamedWindow("filtered",1);
 	cvShowImage("ori",oriImg);
 	cvShowImage("filtered",afterImg);

 	cvWaitKey(0);
 	cvReleaseImage(&oriImg);
 	cvDestroyWindow("image");
 	cvReleaseImage(&afterImg);
 	cvDestroyWindow("filtered");

    return 0;
}

初识API

C++: int meanShift(InputArray probImage,
Rect& window, TermCriteria criteria)
 
  • probImage – 反向投影矩阵
  • window – 一开始的窗口,就是自己用矩形圈起来的地方
  • criteria – Stop criteria for the iterative search algorithm.
   

返回值:收敛时的迭代次数

TermCriteria::TermCriteria

The constructors.

C++: TermCriteria::TermCriteria()
C++: TermCriteria::TermCriteria(int type,
int maxCount, double epsilon)
C++: TermCriteria::TermCriteria(const
CvTermCriteria& criteria)
 
  • type – The type of termination criteria: TermCriteria::COUNTTermCriteria::EPS or TermCriteria::COUNT +TermCriteria::EPS.
  • maxCount – The maximum number of iterations or elements to compute.
  • epsilon – The desired accuracy or change in parameters at which the iterative algorithm stops.
  • criteria – Termination criteria in the deprecated CvTermCriteria format.

举个例子:假设我们对(110,260)坐标处宽35长40的方框感兴趣,最大迭代10次,差值(精确度)为0.01

result是我们的反向投影矩阵,那就可以如下调用

Rect rect(110,260,35,40);

TermCriteria criteria(TermCriteria::MAX_ITER,10,0.01);

meanShift(result,rect,criteria);

荷枪实弹

下面我附上完整代码,并有详细注释

#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/video/tracking.hpp"
#include <iostream>
using namespace cv;
using namespace std;

class ColorHistogram {
private:
    int histSize[3];
	float hranges[2];
    const float* ranges[3];
    int channels[3];

public:
	ColorHistogram() {
		// Prepare arguments for a color histogram
		histSize[0]= histSize[1]= histSize[2]= 256;
		hranges[0]= 0.0;    // BRG range
		hranges[1]= 255.0;
		ranges[0]= hranges; // all channels have the same range
		ranges[1]= hranges;
		ranges[2]= hranges;
		channels[0]= 0;		// the three channels
		channels[1]= 1;
		channels[2]= 2;
	}
	//这里获取HSV的第一个通道的直方图
	Mat getHueHistogram(const cv::Mat &image,int minSaturation=0) {
		Mat hist;
		//色彩空间转换
		Mat hsv;
		cvtColor(image, hsv, CV_BGR2HSV);
		// mask,要注意哦
		Mat mask;
		if (minSaturation>0) {
			// 三通道分开
			vector<Mat> v;
			split(hsv,v);
			// 去除低饱和的值
			threshold(v[1],mask,minSaturation,255,THRESH_BINARY);
		}
		hranges[0]= 0.0;    //hue的range是[0,180]
		hranges[1]= 180.0;
		channels[0]= 0;
		calcHist(&hsv,1,channels,mask,hist,1,histSize,ranges);//注意,mask不为0的地方计算在内,也就是高饱和的Hue计算在内
		return hist;
	}

};

class ContentFinder {
private:
	float hranges[2];
	const float* ranges[3];
	int channels[3];
	float threshold;
	Mat histogram;

public:
	ContentFinder() : threshold(-1.0f) {
		hranges[0]= 0.0;	// range [0,255]
		hranges[1]= 255.0;
		channels[0]= 0;		// the three channels
		channels[1]= 1;
		channels[2]= 2;
		ranges[0]= hranges; // all channels have same range
		ranges[1]= hranges;
		ranges[2]= hranges;
	}

	// Sets the reference histogram
	void setHistogram(const Mat& h) {
		histogram= h;
		normalize(histogram,histogram,1.0);
	}

	cv::Mat find(const cv::Mat& image, float minValue, float maxValue, int *channels, int dim) {
		cv::Mat result;

		hranges[0]= minValue;
		hranges[1]= maxValue;
		ranges[0]= hranges; // all channels have same range
		ranges[1]= hranges;
		ranges[2]= hranges;

		calcBackProject(&image,1,channels,histogram,result,ranges,255.0);
		return result;
	}

};

int main( int, char** argv )
{
	Mat image= cv::imread("baboon1.jpg");
	// Baboon‘s face ROI
	Mat imageROI= image(cv::Rect(110,260,35,40));
	// 获得Hue直方图
	int minSat=65;
	ColorHistogram hc;
	Mat colorhist = hc.getHueHistogram(imageROI,minSat);

        ContentFinder finder;
	finder.setHistogram(colorhist);

	image= cv::imread("baboon3.jpg");
	// 色彩空间转换
	Mat hsv;
	cvtColor(image, hsv, CV_BGR2HSV);

	vector<Mat> v;
	split(hsv,v);
	//除去低饱和度
	threshold(v[1],v[1],minSat,255,cv::THRESH_BINARY);

	// 获得反向投影
	int ch[1]={0};
	Mat result= finder.find(hsv,0.0f,180.0f,ch,1);
	// 除去低饱和度的点
	bitwise_and(result,v[1],result);
	Rect rect(110,260,35,40);
	rectangle(image, rect, Scalar(0,0,255));
	TermCriteria criteria(TermCriteria::MAX_ITER,10,0.01);
	meanShift(result,rect,criteria);

	rectangle(image, rect, cv::Scalar(0,255,0));
	namedWindow("result");
	imshow("result",image);

  	waitKey(0);
  	return 0;
}	

看下效果

举一反三

camshift算法是对meanshift算法的改进,它先用mean shift找到物体,再调整窗的大小,然后可以旋转窗找到最合适的旋转角度

计算机视觉讨论群:162501053

转载请注明:http://blog.csdn.net/abcd1992719g

OpenCV2马拉松第11圈——meanshift与直方图反向投影

时间: 2024-08-27 19:01:25

OpenCV2马拉松第11圈——meanshift与直方图反向投影的相关文章

OpenCV2马拉松第10圈——直方图反向投影(back project)

收入囊中 灰度图像的反向投影 彩色图像的反向投影 利用反向投影做object detect 葵花宝典 什么是反向投影?事实上没有那么高大上! 在上一篇博文学到,图像能够获得自己的灰度直方图. 反向投影差点儿相同是逆过程,由直方图得到我们的投影图. 步骤例如以下: 依据模版图像,得到模版图像的灰度直方图. 对灰度直方图对归一化,归一化后是个概率分布,直方图的积分是1 依据概率分布的直方图,求输入图像的投影图,也就是对每个像素点,我们依据灰度值,能够得到其概率 得到的投影图是介于[0,1]之间的,为

OpenCV2马拉松第8圈——绘制直方图

收入囊中 灰度直方图 彩色直方图 葵花宝典 直方图的理论还是非常丰富的,应用也很多,诸如: 直方图均衡化 直方图匹配(meanshift,camshift) 在这里,我先介绍基础,如何绘制图像的直方图. 拿灰度图像来说,直方图就是不同的灰度对应的个数,横轴(x)就是[0,256), 纵轴(y)就是对应的个数 如下图,分别是灰度直方图和彩色直方图 初识API C++: void calcHist(const Mat* images, int nimages, const int* channels

OpenCV2马拉松第12圈——直方图比較

收入囊中 使用4种不同的方法进行直方图比較 葵花宝典 要比較两个直方图, 首先必需要选择一个衡量直方图相似度的对照标准.也就是先说明要在哪个方面做对照. 我们能够想出非常多办法,OpenCV採用了下面4种 公式也都不难,我们自己就能实现. d越小,表示差异越低,两幅图像越接近,越相似 初识API C++: double compareHist(InputArray H1, InputArray H2, int method) C++: double compareHist(const Spars

OpenCV2马拉松第9圈——再谈对比度(对比度拉伸,直方图均衡化)

收入囊中 lookup table 对比度拉伸 直方图均衡化 葵花宝典 lookup table是什么东西呢? 举个例子,假设你想把图像颠倒一下,f[i] = 255-f[i],你会怎么做? for( int i = 0; i < I.rows; ++i) for( int j = 0; j < I.cols; ++j ) I.at<uchar>(i,j) = 255 - I.at<uchar>(i,j); 大部分人应该都会这么做.或者: for( i = 0; i &

OpenCV2马拉松第12圈——直方图比较

收入囊中 使用4种不同的方法进行直方图比较 葵花宝典 要比较两个直方图, 首先必须要选择一个衡量直方图相似度的对比标准.也就是先说明要在哪个方面做对比. 我们可以想出很多办法,OpenCV采用了以下4种 公式也都不难,我们自己就能实现. d越小,表示差异越低,两幅图像越接近,越相似 初识API C++: double compareHist(InputArray H1, InputArray H2, int method) C++: double compareHist(const Sparse

OpenCV2马拉松第13圈——模版匹配

收入囊中 在http://blog.csdn.net/abcd1992719g/article/details/25505315这里,我们已经学习了如何利用反向投影和meanshift算法来在图像中查找给定模版图片的位置.meanshift针对的是单张图像,在连续图像序列的跟踪中,camshift(Continuously Adaptive Mean-SHIFT)是一种著名的算法.但在这里,我们先不讨论camshift,而是先讨论最简单的模版匹配. 模版匹配算法 opencv normalize

OpenCV2马拉松第27圈——SIFT论文,原理及源码解读

计算机视觉讨论群162501053 转载请注明:http://blog.csdn.net/abcd1992719g/article/details/28913101 简介 SIFT特征描述子是David G. Lowe 在2004年的ijcv会议上发表的论文中提出来的,论文名为<<Distinctive Image Featuresfrom Scale-Invariant Keypoints>>.这是一个很强大的算法,主要用于图像配准和物体识别等领域,但是其计算量相比也比较大,性价

OpenCV2马拉松第16圈——边缘检测(形态学梯度)

计算机视觉讨论群162501053 转载请注明:http://blog.csdn.net/abcd1992719g 收入囊中 利用OpenCV函数进行形态学梯度操作 自定义结构矩阵进行形态学梯度操作 葵花宝典 在此之前,如果你还没接触过灰度图像形态学膨胀与腐蚀,希望你能仔细阅读灰度图像形态学膨胀与腐蚀 本质上,灰度与二值并不差异,二值不过是0与255,膨胀与腐蚀的操作都是一样的 形态学梯度的定义如下: 形态梯度 dst=morph_grad(src,element)=dilate(src,ele

OpenCV2马拉松第17圈——边缘检测(Canny边缘检测)

计算机视觉讨论群162501053 转载请注明:http://blog.csdn.net/abcd1992719g 收入囊中 利用OpenCV Canny函数进行边缘检测 掌握Canny算法基本理论 分享Java的实现 葵花宝典 在此之前,我们先阐述一下canny检测的算法.总共分为4部分. (1)处理噪声 一般用高斯滤波.OpenCV使用如下核 (2)计算梯度幅值 先用如下Sobel算子计算出水平和竖直梯度 我在OpenCV2马拉松第14圈--边缘检测(Sobel,prewitt,robert