《Optimized contrast enhancement for real-time image and video dehazing》大气光估计策略

最近研究去雾算法中接触到了Jin-HwanKim 2013年的paper，如标题，文中提出一种新的大气光估计策略，现对其进行粗略总结。
Optimized contrast enhancement for real-time image and video dehazing原文
文章主页

大气光A估计
传统大气光A通常使用图像中最亮的像素点，因为雾霾使得图像看起来灰蒙蒙的，导致图像灰度值较高。然而，场景中不可避免地可能出现部分本身灰度值就很高的像素点，如果继续使用传统方法会出现估计。
作者提出了一个更加可靠的估计方法，在雾霾区域像素方差值一般较小。即雾霾区域内图像像素间的差异较小，这样可以排除那些存在高亮像素点但细节较多的区域。因此，作者提出了一种称为四叉树子空间划分的层次搜索方法。
具体步骤是：

1. 输入图像：Input，若其尺寸小于固定阈值（如200），进入步骤2，否则进入步骤5。
2. 将图像等分为四个子块，分别计算每个子块R,G,B三分量的均值与方差。
3. 对于每个子块，将对应分量的均值减去方差得到差值，再将三个分量的差值相加，记为每个子块的得分score。
4. 找到得分最大的子块，对此子块进入步骤1进行递归操作。
5. 选择当前块计算大气光，即找到当前块中与（255,255,255）最接近的像素点。
   $$||(Ir(p),Ig(p),Ib(p))-(255,255,255)||$$

matlab实现

close all;
I = imread('test.png'): %load image
img_rgb = double(I);
[hei, wid, ~] = size(img_rgb);

figure, imshow(I); hold on; %hierarchical search display
line([1 wid], [hei/2 hei/2], 'Color', 'g');
line([wid/2 wid/2], [1 hei], 'Color', 'g');
yoffset = 1;
xoffset = 1;

while(hei * wid > 200)
   % divide the previous image into four equal subblocks
    img_ul = img_rgb(1:hei/2, 1:wid/2, :); 
    img_ur = img_rgb(1:hei/2, wid/2+1:wid, :);
    img_ll = img_rgb(hei/2+1:hei, 1:wid/2, :);
    img_lr = img_rgb(hei/2+1:hei, wid/2+1:wid, :);

    % calculate block mean - std
    dpscore(1) = mean2(img_ul(:,:,1)) - std2(img_ul(:,:,1));
    dpscore(2) = mean2(img_ul(:,:,2)) - std2(img_ul(:,:,2));
    dpscore(3) = mean2(img_ul(:,:,3)) - std2(img_ul(:,:,3));
    afscore(1) = sum(dpscore(:));
    
    dpscore(1) = mean2(img_ur(:,:,1)) - std2(img_ur(:,:,1));
    dpscore(2) = mean2(img_ur(:,:,2)) - std2(img_ur(:,:,2));
    dpscore(3) = mean2(img_ur(:,:,3)) - std2(img_ur(:,:,3));
    afscore(2) = sum(dpscore(:));

    dpscore(1) = mean2(img_ll(:,:,1)) - std2(img_ll(:,:,1));
    dpscore(2) = mean2(img_ll(:,:,2)) - std2(img_ll(:,:,2));
    dpscore(3) = mean2(img_ll(:,:,3)) - std2(img_ll(:,:,3));
    afscore(3) = sum(dpscore(:));
    
    dpscore(1) = mean2(img_lr(:,:,1)) - std2(img_lr(:,:,1));
    dpscore(2) = mean2(img_lr(:,:,2)) - std2(img_lr(:,:,2));
    dpscore(3) = mean2(img_lr(:,:,3)) - std2(img_lr(:,:,3));
    afscore(4) = sum(dpscore(:));
    
    % select the region with the highest score and 
    % divide it further into four smaller regions
    [maxscore, maxind] = max(afscore);
    clear img_rgb;
    if(maxind == 1)
      img_rgb = img_ul;
      [hei, wid, ~] = size(img_ul);
    elseif(maxind == 2)
      img_rgb = img_ur;
      [hei, wid, ~] = size(img_ur);
      xoffset = xoffset + wid;
    elseif(maxind == 3)
      img_rgb = img_ll;
      [hei, wid, ~] = size(img_ll);
      yoffset = yoffset + hei;
    elseif(maxind == 4)
      img_rgb = img_lr;
      [hei, wid, ~] = size(img_lr);
      yoffset = yoffset + hei;
      xoffset = xoffset + wid;
    end
    line([xoffset, xoffset+wid], [yoffset+hei/2, yoffset+hei/2], 'Color', 'g');
    line([xoffset+wid/2, xoffset+wid/2], [yoffset, yoffset+hei], 'Color', 'g');
end

if(maxind == 1)
  rectangle('Position', [xoffset, yoffset, wid, hei], 'FaceColor', 'r', 'EdgeColor', 'g');
elseif(maxind == 2)
  rectangle('Position', [xoffset+wid/2, yoffset, wid, hei], 'FaceColor', 'r', 'EdgeColor', 'g');
elseif(maxind == 3)
   rectangle('Position', [xoffset, yoffset+hei/2, wid, hei], 'FaceColor', 'r', 'EdgeColor', 'g');
elseif(maxind == 4)
   rectangle('Position', [xoffset+wid/2, yoffset+hei/2, wid, hei], 'FaceColor', 'r', 'EdgeColor', 'g');
end

% choose the color vector which minimizes the distance as the atmospheric light
nDiff = abs(img_rgb - 255*ones(hei, wid, 3));
nDistance = nDiff(:,:,1) + nDiff(:,:,2) + nDiff(:,:,3);  %here abs-sum can be replaced by square-sum
[rowind, colind] = find(nDistance == max(nDistance(:)));

Airlight(1) = img_rgb(rowind(1));
Airlight(2) = img_rgb(rowind(2));
Airlight(3) = img_rgb(rowind(3));
'