2016-12-08 19:48:28 u011630458 阅读数 5075
  • Java经典算法讲解

    在面试中,算法题目是必须的,通过算法能够看出一个程序员的编程思维,考察对复杂问题的设计与分析能力,对问题的严谨性都能够体现出来。一个算法的好坏,直接影响一个方法调用的性能,进而影响软件的整体性能。算法是学习所有编程语言的基础,在Java的学习过程中首先也会选择以算法起步,本次课程重点讲解Java开发中常用的基本算法。

    29956 人正在学习 去看看 张中强

简介

  本篇主要记录下一个图像自适应对比度增强算法实现。参考论文:a_fast_and_adaptive_method_for_image_contrast_enhancement

实现流程

  详细算法原理请参考论文资料。
  1、拿到待处理图像,以每个待处理像素为中心,3x3或者其他大小的windows,计算出对应窗口下最大值、最小值、平均值。
       注意:直接对原图计算最大、最小、平均值,在结果图像中容易出现块状显像。可以通过对原图做高斯模糊之后在计算处理。
      或者对结果值做如下之类权重操作:
             
             
     2、根据前面的Imin,Imax,Iave,计算出每个像素的权重强度w
             
     3、根据每个像素的权重强度w,计算出Inew和Avenew,再利用这两个参数值,计算出每个像素点对应的参数a
        
     4、基于前面计算结果,再通过如下公式得到结果图像。
                   
               

效果演示

    左边为原始图像,右边为处理后结果对比。
 注:该结果图像,是在原论文基础上做了些变种优化和参数调整后结果。
         
2016-11-04 12:09:02 baimafujinji 阅读数 16714
  • Java经典算法讲解

    在面试中,算法题目是必须的,通过算法能够看出一个程序员的编程思维,考察对复杂问题的设计与分析能力,对问题的严谨性都能够体现出来。一个算法的好坏,直接影响一个方法调用的性能,进而影响软件的整体性能。算法是学习所有编程语言的基础,在Java的学习过程中首先也会选择以算法起步,本次课程重点讲解Java开发中常用的基本算法。

    29956 人正在学习 去看看 张中强

图像去雾哪家强?之前我们已经讨论过了著名的基于暗通道先验的图像去雾(Kaiming He, 2009)算法,如果你用兴趣可以参考:

此外,网上也有很多同道推荐了一篇由韩国学者所发表的研究论文《Optimized contrast enhancement for real-time image and video dehazing》(你也可以从文末参考文献【1】给出的链接中下载到这篇经典论文),其中原作者就提出了一个效果相当不错的图像去雾算法。最近有朋友在我们的图像处理算法研究学习群中也提到了该算法,恰巧想到主页君也很久未发图像相关之文章了,今天遂心血来潮,向大家科普一下这个新算法。

欢迎关注白马负金羁的博客 http://blog.csdn.net/baimafujinji,为保证公式、图表得以正确显示,强烈建议你从该地址上查看原版博文。本博客主要关注方向包括:数字图像处理、算法设计与分析、数据结构、机器学习、数据挖掘、统计分析方法、自然语言处理。


算法核心1:计算大气光值

通常,图像去雾问题的基本模型可以用下面这个公式来表示(这一点在基于暗通道先验的图像去雾中我们也使用过):

I(p)=t(p)J(p)+(1t(p))AI(p)=t(p)J(p)+(1-t(p))A

其中,J(p)=(Jr(p),Jg(p),Jb(p))TJ(p)=(J_r(p),J_g(p),J_b(p))^T 表示原始图像(也就是没有雾的图像);$ I§= (I_r§,I_g§,I_b§)^T$ 表示我们观察到的图像(也就是有雾的图像)。rrggbb 表示位置 pp 处的像素的三个分量。$ A=(A_r, A_g, A_b)^T $ 是全球大气光,它表示周围环境中的大气光。
此外,t(p)[0,1]t(p)\in[0,1] 是反射光的透射率, 由场景点到照相机镜头之间的距离所决定。因为光传播的距离越远,那么通常光就约分散而且越发被削弱。所以上面这个公式的意思就是,本来没有被雾所笼罩的图像 JJ 与大气光 AA 按一定比例进行混合后就得到我们最终所观察到的有雾图像。

大气光 AA 通常用图像中最明亮的颜色来作为估计。因为大量的灰霾通常会导致一个发亮(发白)的颜色。然而,在这个框架下,那些颜色比大气光更加明亮的物体通常会被选中,因而便会导致一个本来不应该作为大气光参考值的结果被用作大气光的估计。为了更加可靠的对大气光进行估计,算法的作者利用了这样一个事实:通常,那些灰蒙蒙的区域(也就是天空)中像素的方差(或者变动)总体来说就比较小。

基于这个认识,算法的作者提出了一个基于四叉树子空间划分的层次搜索方法。如下图所示,我们首先把输入图像划分成四个矩形区域。然后,为每个子区域进行评分,这个评分的计算方法是“用区域内像素的平均值减去这些像素的标准差”(the average pixel value subtracted by the standard deviation
of the pixel values within the region)。记下来,选择具有最高得分的区域,并将其继续划分为更小的四个子矩形。我们重复这个过程直到被选中的区域小于某个提前指定的阈值。例如下图中的红色部分就是最终被选定的区域。在这被选定的区域里,我们选择使得距离 (Ir(p),Ig(p),Ib(p))(255,255,255)||(I_r(p),I_g(p),I_b(p))-(255,255,255)|| 最小化的颜色(包含 r,g,br,g,b 三个分量)来作为大气光的参考值。注意,这样做的意义在于我们希望选择那个离纯白色最近的颜色(也就是最亮的颜色)来作为大气光的参考值。

我们假设在一个局部的小范围内,场景深度是相同的(也就是场景内的各点到相机镜头的距离相同),所以在一个小块内(例如32×3232\times32)我们就可以使用一个固定的透射率 tt,所以前面给出的有雾图像与原始(没有雾的)图像之间的关系模型就可以改写为
J(p)=1t(I(p)A)+AJ(p) = \frac{1}{t}(I(p)-A)+A

可见,在求得大气光 AA 的估计值之后,我们希望复原得到的原始(没有雾的)图像 J(p)J(p) 将依赖于散射率 tt
总的来说,一个有雾的块内,对比度都是比较低的,而被恢复的块内的对比度则随着 tt 的估计值的变小而增大,我们将设法来估计一个最优的 tt 值,从而使得去雾后的块能够得到最大的对比度。

下面是原作者给出的大气光计算函数代码(C/C++版)

void dehazing::AirlightEstimation(IplImage* imInput)
{
	int nMinDistance = 65536;
	int nDistance;

	int nX, nY;

	int nMaxIndex;
	double dpScore[3];
	double dpMean[3];
	double dpStds[3];

	float afMean[4] = {0};
	float afScore[4] = {0};
	float nMaxScore = 0;

	int nWid = imInput->width;
	int nHei = imInput->height;

	int nStep = imInput->widthStep;

	// 4 sub-block
	IplImage *iplUpperLeft = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 3);
	IplImage *iplUpperRight = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 3);
	IplImage *iplLowerLeft = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 3);
	IplImage *iplLowerRight = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 3);

	IplImage *iplR = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 1);
	IplImage *iplG = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 1);
	IplImage *iplB = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 1);

	// divide 
	cvSetImageROI(imInput, cvRect(0, 0, nWid/2, nHei/2));
	cvCopyImage(imInput, iplUpperLeft);
	cvSetImageROI(imInput, cvRect(nWid/2+nWid%2, 0, nWid, nHei/2));
	cvCopyImage(imInput, iplUpperRight);
	cvSetImageROI(imInput, cvRect(0, nHei/2+nHei%2, nWid/2, nHei));
	cvCopyImage(imInput, iplLowerLeft);
	cvSetImageROI(imInput, cvRect(nWid/2+nWid%2, nHei/2+nHei%2, nWid, nHei));
	cvCopyImage(imInput, iplLowerRight);

	// compare to threshold(200) --> bigger than threshold, divide the block
	if(nHei*nWid > 200)
	{
		// compute the mean and std-dev in the sub-block
		// upper left sub-block
		cvCvtPixToPlane(iplUpperLeft, iplR, iplG, iplB, 0);

		cvMean_StdDev(iplR, dpMean, dpStds);
		cvMean_StdDev(iplG, dpMean+1, dpStds+1);
		cvMean_StdDev(iplB, dpMean+2, dpStds+2);
		// dpScore: mean - std-dev
		dpScore[0] = dpMean[0]-dpStds[0];
		dpScore[1] = dpMean[1]-dpStds[1];
		dpScore[2] = dpMean[2]-dpStds[2];

		afScore[0] = (float)(dpScore[0]+dpScore[1]+dpScore[2]);

		nMaxScore = afScore[0];
		nMaxIndex = 0;

		// upper right sub-block
		cvCvtPixToPlane(iplUpperRight, iplR, iplG, iplB, 0);

		cvMean_StdDev(iplR, dpMean, dpStds);
		cvMean_StdDev(iplG, dpMean+1, dpStds+1);
		cvMean_StdDev(iplB, dpMean+2, dpStds+2);
		
		dpScore[0] = dpMean[0]-dpStds[0];
		dpScore[1] = dpMean[1]-dpStds[1];
		dpScore[2] = dpMean[2]-dpStds[2];

		afScore[1] = (float)(dpScore[0]+dpScore[1]+dpScore[2]);
		
		if(afScore[1] > nMaxScore)
		{
			nMaxScore = afScore[1];
			nMaxIndex = 1;
		}

		// lower left sub-block
		cvCvtPixToPlane(iplLowerLeft, iplR, iplG, iplB, 0);

		cvMean_StdDev(iplR, dpMean, dpStds);
		cvMean_StdDev(iplG, dpMean+1, dpStds+1);
		cvMean_StdDev(iplB, dpMean+2, dpStds+2);
		
		dpScore[0] = dpMean[0]-dpStds[0];
		dpScore[1] = dpMean[1]-dpStds[1];
		dpScore[2] = dpMean[2]-dpStds[2];

		afScore[2] = (float)(dpScore[0]+dpScore[1]+dpScore[2]);
		
		if(afScore[2] > nMaxScore)
		{
			nMaxScore = afScore[2];
			nMaxIndex = 2;
		}

		// lower right sub-block
		cvCvtPixToPlane(iplLowerRight, iplR, iplG, iplB, 0);

		cvMean_StdDev(iplR, dpMean, dpStds);
		cvMean_StdDev(iplG, dpMean+1, dpStds+1);
		cvMean_StdDev(iplB, dpMean+2, dpStds+2);
		
		dpScore[0] = dpMean[0]-dpStds[0];
		dpScore[1] = dpMean[1]-dpStds[1];
		dpScore[2] = dpMean[2]-dpStds[2];

		afScore[3] = (float)(dpScore[0]+dpScore[1]+dpScore[2]);

		if(afScore[3] > nMaxScore)
		{
			nMaxScore = afScore[3];
			nMaxIndex = 3;
		}

		// select the sub-block, which has maximum score
		switch (nMaxIndex)
		{
		case 0:
			AirlightEstimation(iplUpperLeft); break;
		case 1:
			AirlightEstimation(iplUpperRight); break;
		case 2:
			AirlightEstimation(iplLowerLeft); break;
		case 3:
			AirlightEstimation(iplLowerRight); break;
		}
	}
	else
	{
		// select the atmospheric light value in the sub-block
		for(nY=0; nY<nHei; nY++)
		{
			for(nX=0; nX<nWid; nX++)
			{
				// 255-r, 255-g, 255-b
				nDistance = int(sqrt(float(255-(uchar)imInput->imageData[nY*nStep+nX*3])*float(255-(uchar)imInput->imageData[nY*nStep+nX*3])
					+float(255-(uchar)imInput->imageData[nY*nStep+nX*3+1])*float(255-(uchar)imInput->imageData[nY*nStep+nX*3+1])
					+float(255-(uchar)imInput->imageData[nY*nStep+nX*3+2])*float(255-(uchar)imInput->imageData[nY*nStep+nX*3+2])));
				if(nMinDistance > nDistance)
				{
					nMinDistance = nDistance;
					m_anAirlight[0] = (uchar)imInput->imageData[nY*nStep+nX*3];
					m_anAirlight[1] = (uchar)imInput->imageData[nY*nStep+nX*3+1];
					m_anAirlight[2] = (uchar)imInput->imageData[nY*nStep+nX*3+2];
				}
			}
		}
	}
	cvReleaseImage(&iplUpperLeft);
	cvReleaseImage(&iplUpperRight);
	cvReleaseImage(&iplLowerLeft);
	cvReleaseImage(&iplLowerRight);

	cvReleaseImage(&iplR);
	cvReleaseImage(&iplG);
	cvReleaseImage(&iplB);
}


或者你也可以使用下面这个Matlab的实现,显而易见代码更加简单(源码由图像算法研究学习群里“薛定谔的猫”提供,略有修改):


function airlight = est_airlight(img)
%compute atmospheric light A through hierarchical
%searching method based on the quad-tree subdivision
global best;
[w,h,z] = size(img);
img = double(img);

if w*h > 200
    lu = img(1:floor(w/2),1:floor(h/2),:);
    ru = img(1:floor(w/2),floor(h/2):h,:);
    lb = img(floor(w/2):w,1:floor(h/2),:);
    rb = img(floor(w/2):w,floor(h/2):h,:);
    
    lu_m_r = mean(mean(lu(:,:,1)));
    lu_m_g = mean(mean(lu(:,:,2)));
    lu_m_b = mean(mean(lu(:,:,3)));
 
    ru_m_r = mean(mean(ru(:,:,1)));
    ru_m_g = mean(mean(ru(:,:,2)));
    ru_m_b = mean(mean(ru(:,:,3)));
    
    lb_m_r = mean(mean(lb(:,:,1)));
    lb_m_g = mean(mean(lb(:,:,2)));
    lb_m_b = mean(mean(lb(:,:,3)));
    
    rb_m_r = mean(mean(rb(:,:,1)));
    rb_m_g = mean(mean(rb(:,:,2)));
    rb_m_b = mean(mean(rb(:,:,3)));
    
    lu_s_r = std2(lu(:,:,1));
    lu_s_g = std2(lu(:,:,2));
    lu_s_b = std2(lu(:,:,3));
 
    ru_s_r = std2(ru(:,:,1));
    ru_s_g = std2(ru(:,:,2));
    ru_s_b = std2(ru(:,:,3));
    
    lb_s_r = std2(lb(:,:,1));
    lb_s_g = std2(lb(:,:,2));
    lb_s_b = std2(lb(:,:,3));
    
    rb_s_r = std2(rb(:,:,1));
    rb_s_g = std2(rb(:,:,2));
    rb_s_b = std2(rb(:,:,3));   
    
    score0 = lu_m_r + lu_m_g + lu_m_b - lu_s_r - lu_s_g - lu_s_b;
    score1 = ru_m_r + ru_m_g + ru_m_b - ru_s_r - ru_s_g - ru_s_b;    
    score2 = lb_m_r + lb_m_g + lb_m_b - lb_s_r - lb_s_g - lb_s_b;    
    score3 = rb_m_r + rb_m_g + rb_m_b - rb_s_r - rb_s_g - rb_s_b;    
    x = [score0,score1,score2,score3];
    if max(x) == score0
        est_airlight(lu);
    elseif max(x) == score1
        est_airlight(ru);
    elseif max(x) == score2
        est_airlight(lb);
    elseif max(x) == score3
        est_airlight(rb);   
    end
else
    for i = 1:w
        for j = 1:h
            nMinDistance = 65536;
            distance = sqrt((255 - img(i,j,1)).^2 +  (255 - img(i,j,2)).^2 + (255 - img(i,j,3)).^2);
            if nMinDistance > distance
					nMinDistance = distance;
					best = img(i,j,:);
            end 
        end
    end        
end
 airlight =best;
end


算法核心2:透射率的计算

我们首先给出图像对比度度量的方法(论文中,原作者给出了三个对比度定义式,我们只讨论其中第一个):
CMSE=p=1N=(Jc(p)Jˉc)2NC_{MSE}=\sum^N_{p=1}=\frac{(J_c(p)-\bar{J}_c)^2}{N}

其中 c{r,g,b}c\in\{r,g,b\} 是颜色通道的索引标签,Jˉc\bar{J}_cJc(p)J_c(p)的平均值,并且p=1,,Np=1,\cdots, NNN 是块中像素的数量。

根据之前给出的有雾图像与原始(没有雾的)图像之间的关系模型
J(p)=1t(I(p)A)+AJ(p) = \frac{1}{t}(I(p)-A)+A
我们可以把上述对比度定义式重新为
CMSE=p=1N=(Ic(p)Iˉc)2t2NC_{MSE}=\sum^N_{p=1}=\frac{(I_c(p)-\bar{I}_c)^2}{t^2N}
其中 Iˉc\bar{I}_cIc(p)I_c(p)的平均值,而且你会发现上述式子也告诉我们对比度是关于 tt 的递减函数。

既然我们希望通过增强对比度的方法来去雾,那么不妨将一个区块B内三个颜色通道上的MSE对比度加总,然后再取负,如下

Econtrast=c{r,g,b}pB(Jc(p)Jˉc)2NB=c{r,g,b}pB(Ic(p)Iˉc)2t2NBE_{contrast}=-\sum_{c\in\{r,g,b\}}\sum_{p\in B}\frac{(J_c(p)-\bar{J}_c)^2}{N_B}=-\sum_{c\in\{r,g,b\}}\sum_{p\in B}\frac{(I_c(p)-\bar{I}_c)^2}{t^2N_B}

由于加了负号,所以取对比度最大就等同于取上式最小。

另外一方面,因为对比度得到增强,可能会导致部分像素的调整值超出了0和255的范围,这样就会造成信息的损失以及视觉上的瑕疵。所以算法作者又提出了一个信息量损失的计算公式:

下面是基于OpenCV实现的示例代码:

float dehazing::NFTrsEstimationColor(int *pnImageR, int *pnImageG, int *pnImageB, int nStartX, int nStartY, int nWid, int nHei)
{
	int nCounter;	
	int nX, nY;		
	int nEndX;
	int nEndY;

	int nOutR, nOutG, nOutB;		
	int nSquaredOut;				
	int nSumofOuts;					
	int nSumofSquaredOuts;			
	float fTrans, fOptTrs;			
	int nTrans;						
	int nSumofSLoss;				
	float fCost, fMinCost, fMean;	
	int nNumberofPixels, nLossCount;

	nEndX = __min(nStartX+m_nTBlockSize, nWid); 
	nEndY = __min(nStartY+m_nTBlockSize, nHei); 

	nNumberofPixels = (nEndY-nStartY)*(nEndX-nStartX) * 3;	

	fTrans = 0.3f;	
	nTrans = 427;

	for(nCounter=0; nCounter<7; nCounter++)
	{
		nSumofSLoss = 0;
		nLossCount = 0;
		nSumofSquaredOuts = 0;
		nSumofOuts = 0;
	
		for(nY=nStartY; nY<nEndY; nY++)
		{
			for(nX=nStartX; nX<nEndX; nX++)
			{
				
				nOutB = ((pnImageB[nY*nWid+nX] - m_anAirlight[0])*nTrans + 128*m_anAirlight[0])>>7; // (I-A)/t + A --> ((I-A)*k*128 + A*128)/128
				nOutG = ((pnImageG[nY*nWid+nX] - m_anAirlight[1])*nTrans + 128*m_anAirlight[1])>>7;
				nOutR = ((pnImageR[nY*nWid+nX] - m_anAirlight[2])*nTrans + 128*m_anAirlight[2])>>7;		

				if(nOutR>255)
				{
					nSumofSLoss += (nOutR - 255)*(nOutR - 255);
					nLossCount++;
				}
				else if(nOutR < 0)
				{
					nSumofSLoss += nOutR * nOutR;
					nLossCount++;
				}
				if(nOutG>255)
				{
					nSumofSLoss += (nOutG - 255)*(nOutG - 255);
					nLossCount++;
				}
				else if(nOutG < 0)
				{
					nSumofSLoss += nOutG * nOutG;
					nLossCount++;
				}
				if(nOutB>255)
				{
					nSumofSLoss += (nOutB - 255)*(nOutB - 255);
					nLossCount++;
				}
				else if(nOutB < 0)
				{
					nSumofSLoss += nOutB * nOutB;
					nLossCount++;
				}
				nSumofSquaredOuts += nOutB * nOutB + nOutR * nOutR + nOutG * nOutG;;
				nSumofOuts += nOutR + nOutG + nOutB;
			}
		}
		fMean = (float)(nSumofOuts)/(float)(nNumberofPixels);  
		fCost = m_fLambda1 * (float)nSumofSLoss/(float)(nNumberofPixels) 
			- ((float)nSumofSquaredOuts/(float)nNumberofPixels - fMean*fMean); 

		if(nCounter==0 || fMinCost > fCost)
		{
			fMinCost = fCost;
			fOptTrs = fTrans;
		}

		fTrans += 0.1f;
		nTrans = (int)(1.0f/fTrans*128.0f);
	}
	return fOptTrs; 
}

在原文中,作者还提供了并行实现的算法(如果你参考作者的源码,你会发现他使用了OpenMP)用于实时对视频进行去雾,这已经超出本文的研究范围,我们不做深究。


实验效果

作者的项目主页(文献[2])中提供有一个基于OpenCV的实现,不过由于OpenCV代码配置起来太麻烦,所以我并未采用。非常感谢图像算法研究群里的同学分享给我MATLAB版的代码(我略有修改)。其实这个算法用MATLAB来写真的非常方便,大约不超过150行即可搞定,特别是MATLAB2014之后集成了导向滤波算法,所以一个函数调用直接搞定,省去很多麻烦。下面是我基于这个MATLAB版程序给出一些测试结果,可供参考(由于我和原作者所使用的参数可能存在差异,所以效果并不保证完全相同)。

主要结论

  • 算法对于天空部分的处理相当到位,更优于Kaiming He的暗通道先验算法;
  • 对比度过大时,图像很容易发暗,可以后期将图片稍微调亮一些。
  • 算法本身是从对比增强的角度来进行去雾操作的,所以你可以看出结果自动带对比度增强加成,所以这个算法所取得的结果通常更加鲜亮。

===========
无冥冥之志者,无昭昭之明,无惛惛之事者,无赫赫之功。


参考文献与推荐阅读材料

[1] Kim, J.H., Jang, W.D., Sim, J.Y. and Kim, C.S., 2013. Optimized contrast enhancement for real-time image and video dehazing. Journal of Visual Communication and Image Representation, 24(3), pp.410-425. (百度云下载链接
[2] 原作者的项目主页——可以下载到基于OpenCV实现的程序源代码(*貌似该网站已经失效)
[3] laviewpbt更早之前发布的一篇研究该算法的文章——该博客上同时包含有很多其他图像去雾方面的文章,推荐有兴趣的读者参考

2016-12-08 19:50:39 u011630458 阅读数 694
  • Java经典算法讲解

    在面试中,算法题目是必须的,通过算法能够看出一个程序员的编程思维,考察对复杂问题的设计与分析能力,对问题的严谨性都能够体现出来。一个算法的好坏,直接影响一个方法调用的性能,进而影响软件的整体性能。算法是学习所有编程语言的基础,在Java的学习过程中首先也会选择以算法起步,本次课程重点讲解Java开发中常用的基本算法。

    29956 人正在学习 去看看 张中强

简介

  本篇主要是记录:基于局域对比度增强的单帧图像插值算法。
 参考论文:1、a_fast_and_adaptive_method_for_image_contrast_enhancement
           2、IMAGE INTERPOLATION USING CONSTRAINED ADAPTIVE CONTRAST ENHANCEMENT TECHNIQUES


具体实现

  1、通过双线性插值,将原图像插值到需要的图像尺寸。
    2、遍历大图,并缩放映射每个像素坐标,到原图像对应像素,以原图像对应像素为中心,3x3之类窗口遍历整个图像,找到每个窗口对应最大值Imax,最小值Imin以及局域平均强度Iavg。
    3、求出参数a=(Imax-Imin)/255。
    4、如果大图当前像素值小于Iavg,则:
       如果大图当前像素值大于Iavg,则: 
    5、计算结果In替换掉大图当前像素旧值。
    6、对新得到的结果大图,进行滤波处理,避免边缘锯齿之类情况。
            
       权重W可以为高斯滤波核之类。

结果显示

  从左往右,依次为原始图像,opencv双线性插值结果,本文结果。


2018-11-15 20:35:45 just_sort 阅读数 563
  • Java经典算法讲解

    在面试中,算法题目是必须的,通过算法能够看出一个程序员的编程思维,考察对复杂问题的设计与分析能力,对问题的严谨性都能够体现出来。一个算法的好坏,直接影响一个方法调用的性能,进而影响软件的整体性能。算法是学习所有编程语言的基础,在Java的学习过程中首先也会选择以算法起步,本次课程重点讲解Java开发中常用的基本算法。

    29956 人正在学习 去看看 张中强

前言

这是OpenCV图像处理专栏的第六篇文章,我们一起来看看何凯明博士这篇获得CVPR 2009的最佳论文。这篇论文的灵感来自于作者两个个观察,第一个是在3D游戏中的雾使得作者坚信人眼有特殊的东西去感知雾,而不仅仅是靠对比度。第二个是作者阅读了之前的一篇去雾方面的论文《Single Image Dehazing》,发现这篇论文中的Dark Object Subtraction可以处理均匀的雾,但是非均匀的就处理不好,所以作者尝试在局部使用了Dark Object Subtraction,然后得到了惊人的效果。

原理

  • 暗通道先验:首先说在绝大多数非天空的局部区域里,某一些像素总会有至少一个颜色通道具有很低的值,也就是说该区域光强是一个很小的值。所以给暗通道下了个数学定义,对于任何输入的图像J,其暗通道可以用下面的公式来表示:在这里插入图片描述
    其中JCJ^C表示彩色图像每个通道,Ω(x)Ω(x)表示以像素X为中心的一个窗口。要求暗通道的图像是比较容易的,先求出每个像素在3个通道的最小值,存到一个二维Mat中(灰度图),然后做一个最小值滤波,滤波的半径由窗口大小决定,这里窗口大小为WindowSizeWindowSize,公式表示为WindowsSize=2Radius+1WindowsSize=2*Radius+1,其中RadiusRadius表示滤波半径。
  • 在这里插入图片描述
    暗通道先验理论得出的结论,这个我不知道如何证明,不过论文给出了几个原因:
  • a)汽车、建筑物和城市中玻璃窗户的阴影,或者是树叶、树与岩石等自然景观的投影;
  • b)色彩鲜艳的物体或表面,在RGB的三个通道中有些通道的值很低(比如绿色的草地/树/植物,红色或黄色的花朵/叶子,或者蓝色的水面);
  • c)颜色较暗的物体或者表面,例如灰暗色的树干和石头。

总之,自然景物中到处都是阴影或者彩色,这些景物的图像的暗原色总是很灰暗的。作者在论文中,统计了5000多副图像的特征,也都基本符合这个先验。因此,我们可以认为它是一条定理。

  • 基于这个先验,就是该论文中最核心的部分了。首先,在计算机视觉和图像处理中,下面这个雾生成模型被广泛的应用:I(x)=J(x)t(x)+A(1t(x))I(x)=J(x)t(x)+A(1-t(x)),其中I(x)I(x)是我们待处理的图像,J(x)J(x)是我们要恢复的没有雾的图像,AA是全球大气光成分,t(x)t(x)为透射率。现在已知了I(X)I(X),我们需要求取J(X)J(X),显然这个不定方程有无数解,所以还需要定义一些先验。
  • 将上式处理变形得到:Ic(x)Ac=t(x)Jc(x)Ac+1t(x)\frac{I^c(x)}{A^c}=t(x)\frac{J^c(x)}{A^c}+1-t(x),其中上标cc代表RGBR、G、B三个通道。然后假设在每一个窗口中透射率t(x)t(x)是一个常数,定义为t(x)^\hat{t(x)}并且AA值已经给定,然后对这个式子左右两边同时取2次最小值,得到下面的式子:
  • 在这里插入图片描述
    其中t(x)^\hat{t(x)}就是公式(8)中那个t(x)t(x)部分,因为我不知道怎么用markdown语法写这个符号。
    上式中,JJ是待求的无雾的图像,根据前述的暗原色先验理论有:在这里插入图片描述
    因此,可推导出:
    在这里插入图片描述
    把式(10)带入式(8)中,得到:
    在这里插入图片描述
    这就是透射率t(x)^\hat{t(x)}的预估值。
    在现实生活中,即使是晴天白云,空气中也存在着一些颗粒,因此,看远处的物体还是能感觉到雾的影响,另外,雾的存在让人类感到景深的存在,因此,有必要在去雾的时候保留一定程度的雾,这可以通过在式(11)中引入一个在[0,1] 之间的系数,则式(11)被修正为:
    在这里插入图片描述
    本推文中所有的测试结果依赖于: ω=0.95ω=0.95
  • 上述的推导是基于A已知的情况下,然而事实是A还不知道呢?A怎么计算呢?在实际中,我们可以借助于暗通道图来从有雾图像中获取该值。具体步骤如下:(1)从暗通道图中按照亮度的大小取前0.1%的像素。(2)在这些位置中,在原始有雾图像II中寻找对应的具有最高亮度的点的值,作为A值。 到这一步,我们就可以进行无雾图像的恢复了。由I(x)=J(x)t(x)+A(1t(x))I(x)=J(x)t(x)+A(1-t(x)),推出 J(x)=(I(x)A)/t(x)+AJ(x)=(I(x)-A)/t(x)+A,现在IAtI、A、t都已经求得了,因此,完全进行出J,也就是去雾后的图像了。当投射图tt的值很小时,会导致JJ的值偏大,从而使得图像整体向白场过度,因此一般可设置一阈值t0t_0,当tt值小于t0t_0时,令t=t0t=t_0,本推文中所有效果图均以t0=0.1t_0=0.1为标准计算得来。
  • 最终的结果计算表示为:
    在这里插入图片描述
    按照上面的公式复现了论文,给几张图片测试结果,都是原图和算法处理后的图这样的顺序:
    原图1去雾后额图片在这里插入图片描述

在这里插入图片描述原图3在这里插入图片描述

代码实现

#include <opencv2/opencv.hpp>
#include <iostream>
#include <algorithm>
#include <vector>
using namespace cv;
using namespace std;

int rows, cols;
//获取最小值矩阵
int **getMinChannel(cv::Mat img){
    rows = img.rows;
    cols = img.cols;
    if(img.channels() != 3){
        fprintf(stderr, "Input Error!");
        exit(-1);
    }
    int **imgGray;
    imgGray = new int *[rows];
    for(int i = 0; i < rows; i++){
        imgGray[i] = new int [cols];
    }
    for(int i = 0; i < rows; i++){
        for(int j = 0; j < cols; j++){
            int loacalMin = 255;
            for(int k = 0; k < 3; k++){
                if(img.at<Vec3b>(i, j)[k] < loacalMin){
                    loacalMin = img.at<Vec3b>(i, j)[k];
                }
            }
            imgGray[i][j] = loacalMin;
        }
    }
    return imgGray;
}

//求暗通道
int **getDarkChannel(int **img, int blockSize = 3){
    if(blockSize%2 == 0 || blockSize < 3){
        fprintf(stderr, "blockSize is not odd or too small!");
        exit(-1);
    }
    //计算pool Size
    int poolSize = (blockSize - 1) / 2;
    int newHeight = rows + blockSize - 1;
    int newWidth = cols + blockSize - 1;
    int **imgMiddle;
    imgMiddle = new int *[newHeight];
    for(int i = 0; i < newHeight; i++){
        imgMiddle[i] = new int [newWidth];
    }
    for(int i = 0; i < newHeight; i++){
        for(int j = 0; j < newWidth; j++){
            if(i < rows && j < cols){
                imgMiddle[i][j] = img[i][j];
            }else{
                imgMiddle[i][j] = 255;
            }
        }
    }
    int **imgDark;
    imgDark = new int *[rows];
    for(int i = 0; i < rows; i++){
        imgDark[i] = new int [cols];
    }
    int localMin = 255;
    for(int i = poolSize; i < newHeight - poolSize; i++){
        for(int j = poolSize; j < newWidth - poolSize; j++){
            localMin = 255;
            for(int k = i-poolSize; k < i+poolSize+1; k++){
                for(int l = j-poolSize; l < j+poolSize+1; l++){
                    if(imgMiddle[k][l] < localMin){
                        localMin = imgMiddle[k][l];
                    }
                }
            }
            imgDark[i-poolSize][j-poolSize] = localMin;
        }
    }
    return imgDark;
}

struct node{
    int x, y, val;
    node(){}
    node(int _x, int _y, int _val):x(_x),y(_y),val(_val){}
    bool operator<(const node &rhs){
        return val > rhs.val;
    }
};

//估算全局大气光值
int getGlobalAtmosphericLightValue(int **darkChannel, cv::Mat img, bool meanMode = false, float percent = 0.001){
    int size = rows * cols;
    std::vector <node> nodes;
    for(int i = 0; i < rows; i++){
        for(int j = 0; j < cols; j++){
            node tmp;
            tmp.x = i, tmp.y = j, tmp.val = darkChannel[i][j];
            nodes.push_back(tmp);
        }
    }
    sort(nodes.begin(), nodes.end());
    int atmosphericLight = 0;
    if(int(percent*size) == 0){
        for(int i = 0; i < 3; i++){
            if(img.at<Vec3b>(nodes[0].x, nodes[0].y)[i] > atmosphericLight){
                atmosphericLight = img.at<Vec3b>(nodes[0].x, nodes[0].y)[i];
            }
        }
    }
    //开启均值模式
    if(meanMode == true){
        int sum = 0;
        for(int i = 0; i < int(percent*size); i++){
            for(int j = 0; j < 3; j++){
                sum = sum + img.at<Vec3b>(nodes[i].x, nodes[i].y)[j];
            }
        }
    }
    //获取暗通道在前0.1%的位置的像素点在原图像中的最高亮度值
    for(int i = 0; i < int(percent*size); i++){
        for(int j = 0; j < 3; j++){
            if(img.at<Vec3b>(nodes[i].x, nodes[i].y)[j] > atmosphericLight){
                atmosphericLight = img.at<Vec3b>(nodes[i].x, nodes[i].y)[j];
            }
        }
    }
    return atmosphericLight;
}

//恢复原图像
// Omega 去雾比例 参数
//t0 最小透射率值
cv::Mat getRecoverScene(cv::Mat img, float omega=0.95, float t0=0.1, int blockSize=15, bool meanModel=false, float percent=0.001){
    int** imgGray = getMinChannel(img);
    int **imgDark = getDarkChannel(imgGray, blockSize=blockSize);
    int atmosphericLight = getGlobalAtmosphericLightValue(imgDark, img, meanModel=meanModel, percent=percent);
    float **imgDark2, **transmission;
    imgDark2 = new float *[rows];
    for(int i = 0; i < rows; i++){
        imgDark2[i] = new float [cols];
    }
    transmission = new float *[rows];
    for(int i = 0; i < rows; i++){
        transmission[i] = new float [cols];
    }
    for(int i = 0; i < rows; i++){
        for(int j = 0; j < cols; j++){
            imgDark2[i][j] = float(imgDark[i][j]);
            transmission[i][j] = 1 - omega * imgDark[i][j] / atmosphericLight;
            if(transmission[i][j] < 0.1){
                transmission[i][j] = 0.1;
            }
        }
    }
    cv::Mat dst(img.rows, img.cols, CV_8UC3);
    for(int channel = 0; channel < 3; channel++){
        for(int i = 0; i < rows; i++){
            for(int j = 0; j < cols; j++){
                int temp = (img.at<Vec3b>(i, j)[channel] - atmosphericLight) / transmission[i][j] + atmosphericLight;
                if(temp > 255){
                    temp = 255;
                }
                if(temp < 0){
                    temp = 0;
                }
                dst.at<Vec3b>(i, j)[channel] = temp;
            }
        }
    }
    return dst;
}

int main(){
    cv::Mat src = cv::imread("/home/zxy/CLionProjects/Acmtest/4.jpg");
    rows = src.rows;
    cols = src.cols;
    cv::Mat dst = getRecoverScene(src);
    cv::imshow("origin", src);
    cv::imshow("result", dst);
    cv::imwrite("../zxy.jpg", dst);
    waitKey(0);
}

参考文章

论文原文:https://ieeexplore.ieee.org/document/5567108
参考博客:https://www.cnblogs.com/Imageshop/p/3281703.html
我的github链接:https://github.com/BBuf/Image-processing-algorithm

后记

关于何凯明博士的暗通道去雾算法就介绍到这里了,希望对你有帮助。代码也可以在我的github获取哦。


欢迎关注我的微信公众号GiantPandaCV,期待和你一起交流机器学习,深度学习,图像算法,优化技术,比赛及日常生活等。
图片.png

2019-06-26 14:01:26 jgj123321 阅读数 535
  • Java经典算法讲解

    在面试中,算法题目是必须的,通过算法能够看出一个程序员的编程思维,考察对复杂问题的设计与分析能力,对问题的严谨性都能够体现出来。一个算法的好坏,直接影响一个方法调用的性能,进而影响软件的整体性能。算法是学习所有编程语言的基础,在Java的学习过程中首先也会选择以算法起步,本次课程重点讲解Java开发中常用的基本算法。

    29956 人正在学习 去看看 张中强
  • 背景介绍

opencv中有逼近多边形曲线的函数approxPolyDP,但是不能一直做一个调包侠,因此我用opencv实现了一篇论文中介绍的两种多边形拟合算法。论文名称:《图像处理中多边形拟合的快速算法》,作者:张帆等,点此下载。在写代码之前,首先要了解两个简单的数学知识。

1、已知两个点的坐标(x_{1},y_{1})(x_{2},y_{2}),如何确定直线方程ax+by+c=0 中a、b、c的数值?

             \tiny y=kx+m\Rightarrow \left\{\begin{matrix} y_{1}=kx_{1}+m\\ y_{2}=kx_{2}+m \end{matrix}\right.\Rightarrow \left\{\begin{matrix} k=\frac{y_{2}-y_{1}}{x_{2}-x_{1}}\\ m=y_{1}-\frac{y_{2}-y_{1}}{x_{2}-x_{1}} x_{1}\end{matrix}\right.\Rightarrow y=\frac{y_{2}-y_{1}}{x_{2}-x_{1}}x+(y_{1}-\frac{y_{2}-y_{1}}{x_{2}-x_{1}} x_{1})

整理上式得:

                                         (y_{1}-y_{2})x+(x_{2}-x_{1})y+(y_{2}x_{1}-x_{2}y_{1})=0

因此:

                                         a=y_{1}-y_{2}b=x_{2}-x_{1}c=y_{2}x_{1}-x_{2}y_{1}

2、点到直线的距离公式:

                                                           d=\left | \frac{Ax_{0}+By_{0}+C}{\sqrt{A^{2}+B^{2}}}\right |

  • Opencv实例

该论文中作者提出的多边形拟合算法,我用opencv实现以后发现找到的边缘并不是目标真实的边缘(本人能力水平有限,可能并没有理解作者的真正思想)。因此我对其算法步骤进行了修改,将第六步中的阈值T变为一个比T小的数值,得到下图所示结果。根据结果图可以看出,还是opencv官方的算法效果比较好,不仅速度快而且精度高。其它两种算法速度较慢,且得到的轮廓有些许变形。

效果图:

                                                                                                   原图

             

                     opencv自带方法                                          迭代端点拟合法                                            多边形拟合法

代码:

/////////////////////////
//opencv4.1.0
/////////////////////////

#include <opencv2/opencv.hpp>
#include <iostream>
#include <math.h>

using namespace cv;
using namespace std;

Mat src, binary, src_gray;
TickMeter tm;


/***** Opencv自带方法*****/
void OpencvFitting(vector<vector<Point>>& contours);

/***** 迭代端点拟合法*****/
float T = 0.9;    //距离阈值
int step = 20;    //步长
vector<Point>ResultContours;
void IterativeEndpointFitting(vector<vector<Point>>& contours);
void fun(int top, int sum, int j, vector<vector<Point>>& contours);
float getDist_P2L(Point pointP, Point pointA, Point pointB);

/***** 多边形拟合法*****/
void PolygonFitting(vector<vector<Point>>& contours);


int main(int argc,char** argv )
{
	//读取图片
	src = imread("3.jpg");
	if (src.empty()){
		printf("读取图像失败\n");
		return -1;
	}
	imshow("输入图片",src);

	//二值化
	cvtColor(src,src_gray,COLOR_BGR2GRAY);
	threshold(src_gray,binary,0,255,THRESH_BINARY|THRESH_OTSU);

	//寻找轮廓,并获取坐标点集
	vector<Vec4i>hireachy;
	vector<vector<Point>> contours;
	findContours(binary,contours,hireachy,RETR_TREE,CHAIN_APPROX_NONE,Point());

	//Opencv自带方法
	tm.start();
	OpencvFitting(contours);
	tm.stop();
	cout << "Opencv自带方法运行时间:" << tm.getTimeMilli() << endl;

	//迭代端点拟合法
	tm.start();
	IterativeEndpointFitting(contours);
	tm.stop();
	cout << "迭代端点拟合法运行时间:" << tm.getTimeMilli() << endl;

	//多边形拟合法
	tm.start();
	PolygonFitting(contours);
	tm.stop();
	cout << "多边形拟合法运行时间:" << tm.getTimeMilli() << endl;

	waitKey(0);
	return 0;
}


void OpencvFitting(vector<vector<Point>>& contours) {

	//逼近多边形曲线
	vector<vector<Point>> poly(contours.size());
	for (size_t t = 0; t<contours.size();t++){
		approxPolyDP(contours[t], poly[t], 5, true);
	}

	//画出轮廓
	Mat OpencvResult = Mat::zeros(src_gray.size(), CV_8UC3);
	for (size_t t = 0; t<contours.size();t++){
		drawContours(OpencvResult, poly, t, Scalar(0, 255, 255), 1, 8, Mat(), 0, Point());
	}
	imshow(" OpencvResult", OpencvResult);
}

void IterativeEndpointFitting(vector<vector<Point>>& contours) {

	for (int j = 0; j<contours[0].size() - 10; j = j + step)
	{
		fun(step, step, j,contours);
	}
	vector<vector<Point>> contours1;
	contours1.push_back(ResultContours);//一维向量点集转换为二维向量点集:vector<Point>  ——  vector<vector<Point>>,否则drawContours函数崩溃

	//画出轮廓
	Mat MyResult = Mat::zeros(src_gray.size(), CV_8UC3);
	for (size_t t = 0; t<contours1.size();t++)
	{
		drawContours(MyResult, contours1, t, Scalar(0, 255, 255), 1, 8, Mat(), 0, Point());
	}
	imshow(" 迭代端点拟合法", MyResult);
}

//top:在步长范围内,距离最大点的序号
//sum:参与计算的点的个数
//j:  步数*步长
void fun(int top, int sum, int j, vector<vector<Point>>& contours)
{
	Point LeftE, RightE;                       //左右两侧点距离直线距离最大点的坐标
	float LeftH = 0, RightH = 0;                   //左右两侧点距离直线距离最大值
	int MaxValueIndexLeft, MaxValueIndexRight; //左右两侧点距离直线距离最大值对应的索引
	for (int i = 0; i <= top;i++)
	{
		//左侧距离计算
		float DistanceLeft = getDist_P2L(contours[0][i + j], contours[0][0 + j], contours[0][top + j]);
		if (DistanceLeft>LeftH)
		{
			LeftH = DistanceLeft;
			LeftE = contours[0][i + j];
			MaxValueIndexLeft = i;
		}

		//右侧距离计算
		if (sum != top && i <= (sum - top))
		{
			float DistanceRight = getDist_P2L(contours[0][i + j + top], contours[0][top + j], contours[0][sum + j]);
			if (DistanceRight>RightH)
			{
				RightH = DistanceRight;
				RightE = contours[0][i + j + top];
				MaxValueIndexRight = i;
			}
		}
	}
	if (LeftH>T)
	{
		ResultContours.push_back(LeftE);
		fun(MaxValueIndexLeft, top, j, contours);
	}
	if (RightH>T && sum != top)
	{
		ResultContours.push_back(RightE);
		fun(MaxValueIndexRight, top, j, contours);
	}
}

/***** 点P到直线AB的距离*****/
float getDist_P2L(Point pointP, Point pointA, Point pointB)
{
	//求直线方程
	int A = 0, B = 0, C = 0;
	A = pointA.y - pointB.y;
	B = pointB.x - pointA.x;
	C = pointA.x*pointB.y - pointA.y*pointB.x;
	//代入点到直线距离公式
	float distance = 0;
	distance = ((float)abs(A*pointP.x + B*pointP.y + C)) / ((float)sqrtf(A*A + B*B));
	return distance;
}

void PolygonFitting(vector<vector<Point>>& contours) {

	int T = 10;
	Point E;
	vector<Point>ResultContours;
	for (int j = 0; j<contours[0].size() / (4 * T); j++) {
		float H = 0;
		for (int i = 2 * T + 4 * T*j; i <= 3 * T + 4 * T*j;i++) {
			float Distance = getDist_P2L(contours[0][i], contours[0][0 + 4 * T*j], contours[0][4 * T + 4 * T*j]);
			if (Distance>H) {
				H = Distance;
				E = contours[0][i];
			}
		}
		if (H>1) {  //按照论文中描述,此处应该是if(H>T),但我做了修改
			ResultContours.push_back(E);
		}
	}

	vector<vector<Point>> contours1;
	contours1.push_back(ResultContours);//一维向量点集转换为二维向量点集:vector<Point>  ——  vector<vector<Point>>,否则drawContours函数崩溃

	//画出轮廓
	Mat MyResult = Mat::zeros(src_gray.size(), CV_8UC3);
	for (size_t t = 0; t<contours1.size();t++)
	{
		drawContours(MyResult, contours1, t, Scalar(0, 255, 255), 1, 8, Mat(), 0, Point());
	}
	imshow(" 多边形拟合法", MyResult);
}

 

 

 

 

 

 

 

 

没有更多推荐了,返回首页