C++(OpenCV) 计算图像局部方差

图像局部方差计算公式：

其中，I(x+i, y+j)表示以(x, y)为中心的坐标点的像素；Ws表示窗口内像素的个数；M(x)表示局部窗口内的均值

具体实现代码：

 float* ImageProcess::getLocalVarV1(Mat& img, int r, float* localVarImg)
{
    // 对图像进行边界填充
	localVarImg = new float[img.cols * img.rows];
	Mat borderImg(Size(img.cols + 2*r, img.rows + 2*r), CV_8UC1, -1);
	cv::copyMakeBorder(img, borderImg, r, r, r, r,BORDER_DEFAULT);

	for(int i = r; i < img.rows + r; i++)
	{
		for(int j = r; j < img.cols + r; j++)
		{
			if(borderImg.data[i*img.cols + j] == 0)
			{
			     localVarImg[(i-r)*img.cols + (j-r)] = 0;
			}
			else
			{
			     double localSum = 0;
			     float  localAve = 0;
			     float  localVar = 0;

			     //局部窗口内像素求和localSum
		         for( int m = i-r; m <= i+r; m++)
			     {
			    	 for(int n = j-r; n <= j+r; n++)
				     {
					      localSum = localSum + borderImg.at<short>(m, n);
	
	     			 }
		    	 }

			     // 局部窗口内像素的均值计算localAve;
			     localAve = localSum / POW(2*r+1);
			     localSum = 0;

			     // 局部窗口内每个像素与均值差求和
			     for(int m = i-r; m <= i+r; m++)
			     {
				     for(int n = j-r; n <= j+r; n++)
				     {
					     localSum = localSum + POW(borderImg.at<short>(m, n) - localAve);
				     }
			     }

			      // 局部方差
			     localVar = localSum*1.0 / (POW(2*r+1));
			     localVarImg[(i-r)*img.cols + (j-r)] = localVar;
			 }
			 
		}
	}
	return localVarImg;
}

函数调用实例：

int main(int argc, char* argv)
{

	// 图像索引编号
	int imgIndex = 9;
	string imageIndex;
	
	//int类型 转 string类型
	stringstream ss;
	ss<<imgIndex;
	ss>>imageIndex;

	string path   = "..//data//";
	string IRFile = "IR_";
	string imgFormat = ".bmp";
	string IRImgFile = path + IRFile + imageIndex + imgFormat;

	Mat IRImg = imread(IRImgFile, 0);
	if(IRImg.empty())
	{

		 CV_Error(CV_StsBadArg, "read image wrong, pease check ");
	}
	if(IRImg.channels() == 3)
	{
		cv::cvtColor(IRImg, IRImg, CV_RGB2GRAY);
	}
	imshow("IRImg", IRImg);

	float* IRWeight  = 0;
	int r = 3;
	IRWeight  = new float[IRImg.cols * IRImg.rows];
	getLocalVar(IRImg, r, IRWeight);
 
    return 0；

}

 

 

 

修改版：

以窗口的方式去遍历每个像素的局部均值，再计算局部方差的方法耗时较长，无法满足实时性的需求。

目前常用的积分图的方式进行加速。积分图是一种与窗口尺寸大小无关的计算方式，在处理局部像素上有明显的优势。

 

Mat_<short> myBoxFilter(Mat_<short> img, int r, Mat_<short>& result)
{

	int iHeight = 0;
	int iWidth  = 0;

    int num  = 0;
	int num2 = 0;
	int num3 = 0;
	int num4 = 0;

	iHeight = img.rows;
	iWidth  = img.cols;

	Mat_<short> rowSum = Mat::zeros(iHeight, iWidth, CV_8UC1);
	Mat_<short> colSum = Mat::zeros(iHeight, iWidth, CV_16SC1);		// 防止数据溢出
	Mat_<short> dstImg = Mat::zeros(iHeight, iWidth, CV_16SC1);

	// 行累加操作 
	// imCum = cumsum(img, 1);
	for(int i = 0; i < iHeight; i++)
	{
		if(i == 0)
		{
			img.row(i).copyTo(rowSum.row(i));
		}
		else
		{
			rowSum.rowRange(i, i+1) = img.rowRange(i, i+1) + rowSum.rowRange(i-1, i);
			cout<<rowSum.rowRange(i, i+1)<<endl;
		}
	}
	//cout<<"rowSum: "<<endl<<rowSum<<endl;

	for(int i = 0; i < r+1; i++)
	{
		rowSum.row(i+r).copyTo(dstImg.row(i));
	}

    // imDst(r + 2 : hei - r , : ) = imCum(2 * r + 2 : hei, : ) - imCum(1 : hei - 2 * r - 1, :);
	Mat temp = rowSum.rowRange(2*r+1, iHeight) - rowSum.rowRange(0, iHeight - 2*r - 1);	

	for(int i = r+1; i < iHeight-r; i++)
	{
		temp.row(num).copyTo(dstImg.row(i));
		num++;
	}

	// 对应matlab Code
	// imDst(hei - r + 1:hei, :) = repmat(imCum(hei, :), [r, 1]) - imCum(hei - 2 * r: hei - r - 1, :);
	Mat temp2;
	temp2 = repeat(rowSum.row(iHeight-1), r, 1) - rowSum.rowRange(iHeight - 2*r -1, iHeight-r-1);
	for(int i = iHeight - r; i < iHeight; i++)
	{
		temp2.row(num2).copyTo(dstImg.row(i));
		num2++;
	}
	
	// imCum = cumsum(imDst, 2);
	// 再次基础上进行列累加操作
	for(int j = 0; j < iWidth; j++)
	{
		if(j == 0)
		{
			dstImg.col(0).copyTo(colSum.col(0));
		}
		else
		{	
			colSum.colRange(j, j+1) = dstImg.colRange(j, j+1) + colSum.colRange(j-1, j);
		}
	}

	Mat_<short> imCum; 
	colSum.copyTo(imCum);

    for(int j = 0; j < r + 1; j++)
	{
		imCum.col(j+r).copyTo(dstImg.col(j));
	}

	Mat_<short> temp3;
	temp3 = imCum.colRange(2*r + 1, iWidth) - imCum.colRange(0, iWidth - 2*r - 1);
	for(int j = r + 1; j < iWidth -r; j++)
	{
		temp3.col(num3).copyTo(dstImg.col(j));
		num3++;
	}

	Mat_<short> temp4;
	temp4 = repeat(imCum.col(iWidth-1), 1, r) - imCum.colRange(iWidth - 2*r -1, iWidth -r-1);
	for(int j = iWidth -r; j < iWidth; j++)
	{
		temp4.col(num4).copyTo(dstImg.col(j));
		num4++;
	}
    dstImg.copyTo(result);

	return result;
}

测试调用实例：

int main(int argc, char* argv[])
{
	Mat img = (Mat_<uchar>(7, 5) << 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 12, 13, 14, 15, 16 );


	int r = 2;
	Mat_<short> imgimg;      // 表示img*img
	Mat_<short> imgimgMean ; // 表示img*img的局部均值图像 
	Mat_<short> imgMean;     // 表示img的局部均值图像
	Mat_<short> rImgMean;    // 表示以半径r为像素值的矩阵的局部均值

	Mat_<short> imgMean_imgMean;
	Mat_<short> localVarImg;

	imgimg = Mat_<short>(img).mul(Mat_<short>(img));
	cout<<"imgimg: "<<imgimg<<endl;
	myBoxFilter(imgimg, r, imgimgMean);	// 输入图像为imgimg， 窗口半径为r, imgimgMean为图像imgimg局部均值
	cout<<"imgimgMean: "<<endl<<imgimgMean<<endl;
	myBoxFilter(img, r, imgMean);
	cout<<"imgMean: "<<imgMean<<endl;
	Mat oneImg = Mat::ones(Size(img.cols, img.rows), CV_8UC1);
	myBoxFilter(oneImg, r, rImgMean);
	cout<<"rImgMean: "<<rImgMean<<endl;

	cv::divide(imgimgMean, rImgMean, imgimgMean);		 
	cv::divide(imgMean, rImgMean, imgMean);
    cout<<"imgMean: "<<endl<<imgMean<<endl;
	cout<<"imgimgMean: "<<endl<<imgimgMean<<endl;
	imgMean_imgMean = imgMean.mul(imgMean);
	cout<<"imgMean_imgMean: "<<imgMean_imgMean<<endl;
	localVarImg = imgimgMean - imgMean_imgMean;

	cout<<"localVarImg: "<<endl<<localVarImg<<endl;

	return 0;
}

OpenCV中Mat的默认类型为uchar类型，其最大值为255。由于存在矩阵的乘法，为了防止在计算中出现溢出的情况，这里使用short类型。该结果都是取整，精度上会缺失。