图像加强综述

时间 2019-12-12 标签图像加强综述

做者：方阳, 转载请注明地址。
文件和代码能够在Github下载, markdown推荐用typora打开。
这篇文章是DIP的第二次做业，对图像加强技术进行综述，目录以下：html

Point Operations
- Image Negative
- Contrast Stretching
- Compression of dynamic range
- Grey level slicing
- Image Subtraction
- Image Averaging
- Histogram
  - Histogram Equalization
  - adaptive histogram equalization
  - Contrast Limited Adaptive Hitogram Equalization（CLAHE）
Mask Operations
- Smoothing operations
- Median Filtering
- sharpening operations
- Derivative operations
Transform operations
- Low pass filtering
- High pass filtering
- Band pass filtering
- Homomorphic filtering
Coloring Operations
- False coloring
- Full color processing
Retinex
- SSR
- MSR
- MSRCR
- Experiment
Dark Channel Prior
Reference

1. Point Operations

1.1 Image Negative

图像反转（Image Negative）在许多应用中都颇有用，例如显示医学图像和用单色正片拍摄屏幕，其想法是将产生的负片用做投影片。python

转换方程：\(T: G(x, y) = L - F(x,y)\)，其中\(L\)是最大强度值，灰度图像\(L\)为255。git

实验结果：github

source code:算法

%% 图像反转
L = 255;
img_orig = imread('mammogram.png');
figure();
subplot(1, 2, 1);
imshow(img_orig);
title('Origin');
img_negative = L - img_orig;
subplot(1, 2, 2);
imshow(img_negative);
title('Negative');

1.2 Contrast Stretching

低对比度的图像多是因为光照不足，图像传感器缺少动态范围，甚至在图像采集过程当中透镜孔径设置错误。对比度拉伸（Contrast Stretching）背后的想法是增长图像处理中灰度的动态范围。markdown

转换公式：\(G(x, y) = g_1 + (g_2 - g_1)/ (f_2 - f_1)[F(x, y) - 1]\)，其中这里\([f1, f2]\)为灰度在新范围\([g1, g2]\)上的映射，这里f1为图像的最小强度值，f2为图像的最大强度值。该函数加强了图像的对比度，显示了均匀的强度分布。app

实验结果：jsp

source code:函数

%% 对比度拉伸
imgContrastStretch=imadjust(imgOrig, [0.1 0.4], [0.3 0.5 ])
figure();
subplot(1, 2, 1);
imshow(imgOrig);
title('Origin');
subplot(1, 2, 2);
imshow(imgContrastStretch);
title('ContrastStretch');

1.3 Compression of dynamic range

有时，处理后的图像的动态范围远远超过显示设备的能力，在这种状况下，只有图像最亮的部分在显示屏上可见. 则须要对图像进行动态范围压缩（Compression of dynamic range）。ui

转换公式：\(s = c log(1 + |r|)\)， \(c\)是度量常数，是对数函数执行所需的压缩。\(r\)表示当前像素的灰度，\(s\)为转换后该像素的灰度。例如将如下图像的\([0. 255]\)压缩到\([0, 150]\)，其中\(c = \frac{150}{log(1 + 255)}\).

实验结果：

source code:

%% 动态范围压缩
img_orig1 = rgb2gray(imread('circle.png'));
c = 150 / log(1 + 255);
img_size = size(img_orig1);
for i = 1:img_size(1)
    for j = 1:img_size(2)
        s(i, j) = c * log(double(1 + img_orig1(i, j)));
    end
end
figure();
subplot(1, 2, 1);
imshow(img_orig1);
title('Origin');
subplot(1, 2, 2);
imshow(uint8(s));
title('Compression');

1.4 Grey level slicing

灰度切片(Grey level slicing)至关于带通滤波的空间域。灰度切片函数既能够强调一组灰度值而减小其余全部灰度值，也能够强调一组灰度值而不考虑其余灰度值。例如对图像中灰度值为[100, 180]的区域进行强调，对其余区域进行抑制。

实验结果：

source code:

%% 灰度切片
img_orig2 = rgb2gray(imread('GLS.png'));
figure();
subplot(1, 2, 1);
imshow(img_orig2);
title('Origin');
img_size = size(img_orig2);
for i = 1:img_size(1)
    for j = 1:img_size(2)
        if img_orig2(i, j) > 100 || img_orig2(i, j) > 180
            img_orig2(i, j) = 150;
        else
            img_orig2(i, j) = 25;
        end
    end
end
subplot(1, 2, 2);
imshow(img_orig2);
title('Gray level slicing');

1.5 Image Subtraction

图像相减是从另外一个图像中减去一个像素或整个图像的数字数值的过程。这主要是出于两个缘由——调平图像的不均匀部分和检测两幅图像之间的变化。一个经常使用的方法是从场景中减去背景光照的变化，以便更容易地分析其中的前景对象.例如在捕获过程当中光照不好的文本，以便在整个图像中有很强的光照梯度，若是咱们但愿将前景文本从背景页面中分离出来，也许咱们不能调整光照，可是咱们能够在场景中放入不一样的东西。例如，显微镜成像一般就是这样。因此咱们用一张白纸替换文本，在不改变任何东西的状况下，咱们捕捉到一个新的图像。咱们能够从原始图像中减去光场图像，试图消除背景强度的变化。
实验结果以下：

source code:

%% 图像相减
img_orig3 = rgb2gray(imread('son1.png'));
img_add = uint16(img_orig3) + 100;
whitepaper = rgb2gray(imread('son2.png'));
img_sub = uint8(img_add - uint16(whitepaper));
figure();
subplot(1, 3, 1);
imshow(img_orig3);
title('Origin');
subplot(1, 3, 2);
imshow(whitepaper);
title('whitepaper');
subplot(1, 3, 3);
imshow(img_sub);
title('Subtraction');

1.6 Image Averaging

图像平均是经过找到K个图像的平均值来得到的。应用于图像去噪。
转换公式：\(\overline{g}(x, y)=\frac{1}{K} \sum_{i=1}^{K} g_{i}(x, y)\)，噪声图像定义为\(g(x, y)=f(x, y)+\eta(x, y)\)，假设噪声与零均值不相关。（下图只显示了一张噪声图片）
实验结果:

source code:

%% 图像平均
img_orig4 = rgb2gray(imread('cat.jpg'));
img_noise1 = imnoise(img_orig4, 'gaussian', 0, 0.01);
img_noise2 = imnoise(img_orig4, 'gaussian', 0, 0.01);
img_noise3 = imnoise(img_orig4, 'gaussian', 0, 0.01);
img_noise4 = imnoise(img_orig4, 'gaussian', 0, 0.01);
img_noise5 = imnoise(img_orig4, 'gaussian', 0, 0.01);
img_average = imlincomb(0.2,img_noise1, 0.2,img_noise2, 0.2,img_noise3, 0.2,img_noise4, 0.2,img_noise5);
figure();
subplot(1, 3, 1);
imshow(img_orig4);
title('Origin');
subplot(1, 3, 2);
imshow(img_noise1);
title('img_noise1');
subplot(1, 3, 3);
imshow(img_average);
title('img_average');

1.7 Histogram

1.7.1 Histogram Equalization

直方图均衡化是一种经常使用的图像加强技术。假设咱们有一个主要是黑色的图像。而后它的直方图会向灰度的下端倾斜，全部的图像细节都被压缩到直方图的暗端。若是能将暗端的灰度“拉伸”，生成一个更均匀分布的直方图，那么图像就会清晰得多。

具体作法：1. 求出原图的直方图分布 2.计算原图直方图的累计几率分布 3. 映射，公式能够表示为\(f\left(D_{A}\right)=\frac{L}{A_{0}} \sum_{u=0}^{D_{A}} H_{A}(u)\)，\(A\)为原图，\(H\)为直方图，\(L\)为灰度级，\(A_0\)为像素点个数。

1.7.2 adaptive histogram equalization

直方图均衡化中，是直接对全局图像进行均衡化，是Global Histogram Equalization，而没有考虑到局部图像区域(Local Region)，自适应直方图均衡化（AHE）就是在均衡化的过程当中只利用局部区域窗口内的直方图分布来构建映射函数首先，最简单而且直接的想法，对A图像每一个像素点进行遍历，用像素点周围W * W的窗口进行计算直方图变换的CDF（累计几率分布），而后对该像素点进行映射。[9]

1.7.3 Contrast Limited Adaptive Hitogram Equalization（CLAHE）

CLAHE相对于AHE，提出了两个改进的地方。

第一，提出一种限制直方图分布的方法。考虑图像A的直方图，设定一个阈值，假定直方图某个灰度级超过了阈值，就对之进行裁剪，而后将超出阈值的部分平均分配到各个灰度级，这个过程能够用下图[10]来进行解释。图中左上图是原来的直方图分布，能够看出有两处峰值，其对应的CDF为左下图，能够看出有两段梯度比较大，变化剧烈。对于以前频率超过了阈值的灰度级，那么就把这些超过阈值的部分裁剪掉平均分配到各个灰度级上，如右上图，那么这会使得CDF变得较为平缓，如右下图。一般阈值的设定能够直接设定灰度级出现频数，也能够设定为占总像素比例，后者更容易使用。因为右下图所示的CDF不会有太大的剧烈变化，因此能够避免过分加强噪声点。[9]

第二，提出了一种插值的方法，加速直方图均衡化。首先，将图像分块，每块计算一个直方图CDF，其次，对于图像的每个像素点，找到其邻近的四个窗口，分别计算四个窗口直方图CDF对该点像素点的映射值，记做\(f_{u l}(D), f_{u r}(D), f_{d l}(D), f_{d r}(D)\)，而后进行双线性插值获得最终该像素点的映射值，双线性插值(BiLinear)公式为\(f(D)=(1-\Delta y)\left((1-\Delta x) f_{u l}(D)+\Delta x f_{b l}(D)\right)+\Delta y\left((1-\Delta x) f_{u r}(D)+\Delta x f_{b r}(D)\right)\)

其中\(\Delta x, \Delta y\)是像素点相对于左上角窗口中心，即左上角黑色像素点的距离与窗口大小的比值。[9]

实验结果(从左往右依次是原图，HE，AHE，CLAHE）：

source code(hw_1):

![a](assets/a.png)//HE,AHE,CLAHE
import numpy as np
import argparse
import cv2

ap = argparse.ArgumentParser()
ap.add_argument('--image', required=True)
args = vars(ap.parse_args())

image = cv2.imread(args["image"])
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

eh_image = cv2.equalizeHist(image)

ahe = cv2.createCLAHE(clipLimit=0.0, tileGridSize=(8,8))
ahe_image = ahe.apply(image)

clahe = cv2.createCLAHE(clipLimit=10.0, tileGridSize=(8,8))
clahe_image = clahe.apply(image)

cv2.imshow("Histogram Equalization", np.hstack([image, eh_image, 
                                                ahe_image, clahe_image]))
cv2.waitKey(0)

2. Mask Operations

掩模算子经过在图像上滑动掩模，将掩模值与落在它们下面的像素值相乘并得到总和，来执行掩模与图像的卷积。

能够表示为\(g(x, y)=\sum_{s=-a}^{a} \sum_{t=-b}^{b} w(s, t) f(x+s, y+t)\)，总和用做图像上掩模中心位置的值，\(w(s, t)\)即为掩模算子。

2.1 Smoothing operations

图像平滑是指用于突出图像的宽大区域、低频成分、主干部分或抑制图像噪声和干扰高频成分的图像处理方法，目的是使图像亮度平缓渐变，减少突变梯度，改善图像质量。\(\begin{array}{|c|c|c|}\hline 1 & {1} & {1} \\ \hline 1 & {1} & {1} \\ \hline 1 & {1} & {1} \\ \hline\end{array}\)x\(\frac{1}{9}\)就是一个用来平滑图像的掩模算子。

实验结果：

source code:

img_orig5 = imread('a.png');
H1 = fspecial('average', 3);
img_smooth1 = imfilter(img_orig5, H1);
H2 = fspecial('average', 7);
img_smooth2 = imfilter(img_orig5, H2);
H3 = fspecial('average', 11);
img_smooth3 = imfilter(img_orig5, H3 );
figure();
subplot(2, 2, 1);
imshow(img_orig5);
title('Origin');
subplot(2, 2, 2);
imshow(img_smooth1);
title('kernel=3');
subplot(2, 2, 3);
imshow(img_smooth2);
title('kernel=7');
subplot(2, 2, 4);
imshow(img_smooth3);
title('kernel=11');

2.2 Median Filtering

中值滤波是基于排序统计理论的一种能有效抑制噪声的非线性信号处理技术，中值滤波的基本原理是把数字图像或数字序列中一点的值用该点的一个邻域中各点值的中值代替，让周围的像素值接近的真实值，从而消除孤立的噪声点。

实验结果：

source code:

%% 中值滤波
img_orig6 = rgb2gray(imread('lena.png'));
img_noise6=imnoise(img_orig6, 'salt & pepper', 0.02);
img_recover = medfilt2(img_noise6);
figure();
subplot(1, 3, 1);
imshow(img_orig6);
title('Origin');
subplot(1, 3, 2);
imshow(img_noise6);
title('img_salt & pepper');
subplot(1, 3, 3);
imshow(img_recover);
title('img_recover');

2.3 sharpening operations

图像锐化(image sharpening)是补偿图像的轮廓，加强图像的边缘及灰度跳变的部分，使图像变得清晰，分为空间域处理和频域处理两类。图像锐化是为了突出图像上地物的边缘、轮廓，或某些线性目标要素的特征。这种滤波方法提升了地物边缘与周围像元之间的反差，所以也被称为边缘加强。实验用的sobel算子对图像进行锐化。

实验结果：

source code:

%% 锐化
img_orig7 = imread('t.png');
H5 = fspecial('sobel');
edge = imfilter(img_orig7, H5);
sharpened = img_orig7 + edge;
figure();
subplot(1, 3, 1);
imshow(img_orig7);
title('Origin');
subplot(1, 3, 2);
imshow(edge);
title('Edge');
subplot(1, 3, 3);
imshow(sharpened);
title('img_sharpened');

2.4 Derivative operations

常见的偏导算子有sobel算子和laplacian算子，sobel算子是一阶偏导算子，laplacian算子是二阶偏导算子。

sobel算子有两个方向\(Gx = \left[ \begin{array}{ccc}{-1} & {0} & {+1} \\ {-2} & {0} & {+2} \\ {-1} & {0} & {+1}\end{array}\right]\)，\(Gy = \left[ \begin{array}{ccc}{-1} & {-2} & {-1} \\ {0} & {0} & {0} \\ {+1} & {+2} & {+1}\end{array}\right]\)，laplacian常见的算子有四邻域\(\left[ \begin{array}{ccc}{ 0 } & { 1 } & { 0 } \\ { 1 } & { -4 } & { 1 } \\ { 0 } & { 1 } & { 0 }\end{array}\right]\), 八领域\(\left[ \begin{array}{ccc}{ 1 } & { 1 } & { 1 } \\ { 1 } & { -8 } & { 1 } \\ { 1 } & { 1 } & { 1 }\end{array}\right]\).

实验结果：

%% 偏导
img_orig7 = rgb2gray(imread('lena.png'));
H6 = fspecial('laplacian');
H7 = fspecial('sobel');
laplacian = imfilter(img_orig7, H6);
sobel = imfilter(img_orig7, H7);
figure();
subplot(1, 3, 1);
imshow(img_orig7);
title('Origin');
subplot(1, 3, 2);
imshow(laplacian);
title('laplacian');
subplot(1, 3, 3);
imshow(sobel);
title('sobel');

3. Transform operations

这一节是将图片从空间域转到频域，在频域进行操做，公式表示为\(G(u, v)=F(u, v) H(u, v)\)，\(F(u, v)\)是原图通过快速傅里叶变换转换（Fast Fourier Transform, 简称fft）到频域的频谱，\(H(u, v)\)是在频域执行的操做，\(G(u,v)\)是在频域处理后的频谱结果，最后\(G(u, v)\)能够经过快速傅里叶反变换(ifft)获得滤波的图像。

3.1 Low pass filtering

流程：1）原始正常的图像，加噪处理，获得img_noise； 2） img_noise图像进行傅里叶变换，获得频谱； 3）对获得的频谱进行理想低通滤波，低于截止频率\(d_0\)的经过，高于的抑制； 4）对滤波后的频谱进行反傅里叶变换，获得滤波后图像。[2]

实验结果：

source code:

%% 低通滤波, ref[2]
img_origin=rgb2gray(imread('lena.png'));
d0=50;  %阈值
img_noise=imnoise(img_origin,'salt'); % 加椒盐噪声
img_f=fftshift(fft2(double(img_noise)));  %傅里叶变换获得频谱
[m n]=size(img_f);
m_mid=fix(m/2);  
n_mid=fix(n/2);  
img_lpf=zeros(m,n);

for i=1:m
    for j=1:n
        d=sqrt((i-m_mid)^2+(j-n_mid)^2);   %理想低通滤波，求距离
        if d<=d0
            h(i,j)=1;
        else
            h(i,j)=0;
        end
        img_lpf(i,j)=h(i,j)*img_f(i,j);  
    end
end

img_lpf=ifftshift(img_lpf);    %反傅里叶变换
img_lpf=uint8(real(ifft2(img_lpf)));  %取实数部分

subplot(1,3,1);imshow(img_origin);title('原图');
subplot(1,3,2);imshow(img_noise);title('噪声图');
subplot(1,3,3);imshow(img_lpf);title('理想低通滤波');

3.2 High pass filtering

高通滤波与低通滤波相似，区别在于低于截止频率的抑制，高于截止频率的经过。

实验结果：

source code:

%% 高通滤波
img_origin8 = rgb2gray(imread('lena.png'));
g= fftshift(fft2(double(img_origin8)));
[N1,N2]=size(g);
n=2;
d0=30; 
%d0是终止频率
n1=fix(N1/2);
n2=fix(N2/2);
%n1，n2指中心点的坐标，fix()函数是往0取整
for i=1:N1
  for j=1:N2
      d=sqrt((i-n1)^2+(j-n2)^2);  
    if d>=d0
        h=1;  
    else
        h=0;  
    end  
    result(i,j)=h*g(i,j); 
  end
end
final=ifft2(ifftshift(result));
final=uint8(real(final));
figure();
subplot(2,2,1); imshow(img_origin8); title('原图');
subplot(2,2,2); imshow(abs(g),[]); title('原图频谱');
subplot(2,2,3); imshow(final); title('高通滤波后的图像');
subplot(2,2,4); imshow(abs(result), []); title('高通滤波后的频谱');

3.3 Band pass filtering

带通滤波有两个截止频率\(d_0\), \(d_1\)，其中\(d_0\)是较低的频率，\(d_1\)是较高的频率，图像频谱在\([d_0, d_1]\)之间的经过，在区间以外的抑制。

实验结果：

source code:

%% 带通滤波
img_origin9=rgb2gray(imread('lena.png')); 
g= fftshift(fft2(double(img_origin9))); 
[N1,N2]=size(g);  
n=2;  
d0=0;  
d1=200;  
n1=floor(N1/2);  
n2=floor(N2/2);  
for i=1:N1  
   for j=1:N2  
    d=sqrt((i-n1)^2+(j-n2)^2);  
    if d>=d0 || d<=d1  
        h=1;  
    else
        h=0;  
    end  
    result(i,j)=h*g(i,j);
   end  
end 
final=ifft2(ifftshift(result));
final=uint8(real(final));
figure();
subplot(2,2,1); imshow(img_origin9); title('原图');
subplot(2,2,2); imshow(abs(g),[]); title('原图频谱');
subplot(2,2,3); imshow(final); title('带通滤波后的图像');
subplot(2,2,4); imshow(abs(result), []); title('带通滤波后的频谱');

3.4 Homomorphic filtering

同态滤波是把频率过滤和灰度变换结合起来的一种图像处理方法，它依靠图像的照度/ 反射率模型做为频域处理的基础，利用压缩亮度范围和加强对比度来改善图像的质量。使用这种方法可使图像处理符合人眼对于亮度响应的非线性特性，避免了直接对图像进行傅立叶变换处理的失真。[8]

同态滤波的基本原理是：将像元灰度值看做是照度和反射率两个组份的产物。因为照度相对变化很小，能够看做是图像的低频成份，而反射率则是高频成份。经过分别处理照度和反射率对灰度值的影响，达到揭示阴影区细节特征的目的。[8]

流程：\(S(x, y) ---> Log ---> FFT--->filter--->IFFT--->Exp--->T(x, y)\)[8]

实验结果：

%% 同态滤波 ref[8]
img_origin10 = rgb2gray(imread('water.png')); 
[M,N]=size(img_origin10);
rL=0.5;
rH=4.7;
c=2;
d0=10;
log_img=log(double(img_origin10)+1);
FI=fft2(log_img);
n1=floor(M/2);
n2=floor(N/2);
for i=1:M
    for j=1:N
        D(i,j)=((i-n1).^2+(j-n2).^2);
        H(i,j)=(rH-rL).*(exp(c*(-D(i,j)./(d0^2))))+rL;%高斯同态滤波
    end
end
G = H.*FI;
final=ifft2(G);
final=real(exp(final));
figure();
subplot(2,2,1); imshow(img_origin10); title('原图');
subplot(2,2,2); imshow(abs(FI),[]); title('原图频谱');
subplot(2,2,3); imshow(final, []); title('同态滤波后的图像');
subplot(2,2,4); imshow(abs(G), []); title('同态滤波后的频谱');

4. Coloring Operations

4.1 False coloring

将彩色图像转换为灰度图像是一个不可逆的过程，灰度图像也不可能变换为原来的彩色图像。而某些场合须要将灰度图像转变为彩色图像；伪彩色处理主要是把黑白的灰度图像或者多波段图像转换为彩色图像的技术过程。其目的是提升图像内容的可辨识度。

伪彩色图像的含义是，每一个像素的颜色不是由每一个基色份量的数值直接决定，而是把像素值看成彩色查找表(事先作好的)的表项入口地址，去查找一个显示图像时使用的R，G，B强度值，用查找出的R，G，B强度值产生的彩色称为伪彩色。

实验结果：

source code:

%% 伪彩色
img_origin11 = rgb2gray(imread('lena.png')); 
FalseRGB = label2rgb(gray2ind(img_origin11, 255),jet(255));
figure();
subplot(1,2,1); imshow(img_origin11); title('原图');
subplot(1,2,2); imshow(FalseRGB); title('伪彩色');

4.2 Full color processing

上面处理的都是灰度图像，若是要处理全彩色图像，则须要对彩色的每一个通道分别处理，而后叠加在一块儿。下面以中值滤波为例，对彩色图像进行处理。

实验结果:

source code:

%% 全彩色处理，中值滤波为例
img_orig6 = imread('lena.png');
for i = 1:3
img_noise6(:, :, i) = imnoise(img_orig6(:, :, i), 'salt & pepper', 0.02);
img_recover(:, :, i) = medfilt2(img_noise6(:, :, i));
end
figure();
subplot(1, 3, 1);
imshow(img_orig6);
title('Origin');
subplot(1, 3, 2);
imshow(img_noise6);
title('img_salt & pepper');
subplot(1, 3, 3);
imshow(img_recover);
title('img_recover');

5. Retinex

视网膜-大脑皮层（Retinex）理论认为世界是无色的，人眼看到的世界是光与物质相互做用的结果，也就是说，映射到人眼中的图像和光的长波（R）、中波（G）、短波（B）以及物体的反射性质有关。基于这个理论，能够抽象下图中的公式\(I(x, y)=R(x, y) \bullet L(x, y)\)，\(I(x, y)\)表明观察到的图像，\(R(x, y)\)表明物体的反射属性，\(L(x, y)\)表明物体表面的光照。

5.1 Single-Scale Retinex

在Retinex理论中，一个假定是光照\(I(x, y)\)是缓慢变化的，也就是低频的，要从\(I(x, y)\)中还原出\(R(x, y)\)，因此能够经过低通滤波器获得光照份量。

具体作法：

对源公式两边同时取对数，获得\(\log (R)=\log (I)-\log (L)\)
经过高斯模糊低通滤波器对原图进行滤波，公式为\(L = F * I\)
获得\(L\)后，利用第一步的公式获得\(log(R)\)，而后经过\(exp\)函数获得R

5.2 Multi-Scale Retinex

上面是单个尺度下的Retinex算法，固然也存在多尺度的Retinex算法，最为经典的就是3尺度的，大、中、小，既能实现图像动态范围的压缩，又能保持色感的一致性较好。[11]

具体作法: 对每一个尺度分别进行单尺度的Retinex算法，将每一个尺度的结果相加取平均获得最终结果(这里是简单地取权重相同，固然还能够设立权重不等）.

5.3 Multi-Scale Retinex with Color Restoration

在前面的加强过程当中，图像可能会由于增长了噪声，而使得图像的局部细节色彩失真，不能显现出物体的真正颜色，总体视觉效果变差。带色彩恢复因子C的多尺度算法是在多个固定尺度的基础上考虑色彩不失真恢复的结果，在多尺度Retinex算法过程当中，经过引入一个色彩因子C来弥补因为图像局部区域对比度加强而致使的图像颜色失真的缺陷，一般状况下所引入的色彩恢复因子C的表达式为：
\[ R_{M S R C R_{i}}(x, y)=C_{i}(x, y) R_{M S R_{i}}(x, y) \]

\[ C_{i}(x, y)=f\left[\frac{I_{i}(x, y)}{\sum_{j=1}^{N} I_{j}(x, y)}\right] \]

其中，\(C_i\)表示为第\(i\)个颜色通道地色彩恢复系数，它的做用是调节3个通道颜色的比例，\(f(\cdot)\)表示颜色空间的映射函数，一般能够采用下面的形式：
\[ C_{i}(x, y)=\beta \log \frac{\alpha I_{i}(x, y)}{\sum_{j=1}^{N} I_{j}(x, y)}=\beta\left\{\log \left[\alpha I_{i}\right]-\log \left[\sum_{j=1}^{N} I_{j}(x, y)\right]\right\} \]
其中，\(\beta\)是增益常数，\(\alpha\)是受控制的非线性强度，带色彩恢复的多尺度Retinex算法经过色彩恢复因子C这个系数来调整原始图像中3个颜色通道之间的比例关系，从而把相对有点暗的区域的信息凸显出来，以达到消除图像色彩失真缺陷的目的。处理后的图像局域对比度得以提升，并且其亮度与真实的场景很类似，图像在人们的视觉感知下显得更为逼真。所以，MSRCR算法会具备比较好的颜色再现性、亮度恒常性与动态范围压缩等特性。[12]

5.4 Experiment

source code:

#ref [11]
import argparse
import numpy as np
import cv2


def singleScaleRetinexProcess(img, sigma):
    temp = cv2.GaussianBlur(img, (0, 0), sigma)
    gaussian = np.where(temp == 0, 0.01, temp)
    retinex = np.log10(img + 0.01) - np.log10(gaussian)
    return retinex

def multiScaleRetinexProcess(img, sigma_list):
    retinex = np.zeros_like(img * 1.0)
    for sigma in sigma_list:
        retinex = singleScaleRetinexProcess(img, sigma)
    retinex = retinex / len(sigma_list)
    return retinex

def colorRestoration(img, alpha, beta):
    img_sum = np.sum(img, axis=2, keepdims=True)
    color_restoration = beta * (np.log10(alpha * img) - np.log10(img_sum))
    return color_restoration

def multiScaleRetinexWithColorRestorationProcess(img, sigma_list, G, b, alpha, beta):
    img = np.float64(img) + 1.0
    img_retinex = multiScaleRetinexProcess(img, sigma_list)
    img_color = colorRestoration(img, alpha, beta)
    img_msrcr = G * (img_retinex * img_color + b)
    return img_msrcr

def simplestColorBalance(img, low_clip, high_clip):
    total = img.shape[0] * img.shape[1]
    for i in range(img.shape[2]):
        unique, counts = np.unique(img[:, :, i], return_counts=True)
        current = 0
        for u, c in zip(unique, counts):
            if float(current) / total < low_clip:
                low_val = u
            if float(current) / total < high_clip:
                high_val = u
            current += c
        img[:, :, i] = np.maximum(np.minimum(img[:, :, i], high_val), low_val)
    return img

def touint8(img):
    for i in range(img.shape[2]):
        img[:, :, i] = (img[:, :, i] - np.min(img[:, :, i])) / \
                       (np.max(img[:, :, i]) - np.min(img[:, :, i])) * 255
    img = np.uint8(np.minimum(np.maximum(img, 0), 255))
    return img

def SSR(img, sigma=300):
    ssr = singleScaleRetinexProcess(img, sigma)
    ssr = touint8(ssr)
    return ssr

def MSR(img, sigma_list=[15, 80, 250]):
    msr = multiScaleRetinexProcess(img, sigma_list)
    msr = touint8(msr)
    return msr

def MSRCR(img, sigma_list=[15, 80, 250], G=5, b=25, alpha=125, beta=46, low_clip=0.01, high_clip=0.99):
    msrcr = multiScaleRetinexWithColorRestorationProcess(img, sigma_list, G, b, alpha, beta)
    msrcr = touint8(msrcr)
    msrcr = simplestColorBalance(msrcr, low_clip, high_clip)
    return msrcr

def main():
    ap = argparse.ArgumentParser()
    ap.add_argument('--image', required=True)
    args = vars(ap.parse_args())

    image = cv2.imread(args["image"])

    ssr = SSR(image)
    msr = MSR(image)
    msrcr = MSRCR(image)

    cv2.imshow("Retinex", np.hstack([image, ssr, msr, msrcr]))
    cv2.waitKey(0)


if __name__ == "__main__":
    main()

实验结果：

6. Dark Channel Prior

暗通道先验（Dark Channel Prior）是说在绝大多数非天空的局部区域内，某一些像素至少一个颜色通道具备很低的值，这是何凯明等人基于5000多张天然图像的统计获得的定理。根据这必定理，何凯明等人提出了暗通道先验的去雾算法[13]，

以天然图像和雾气图像为例[14]（左边是原图，右边是暗通道）：

暗通道先验公式以下所示：
\[ J^{d a r k}(\mathbf{x})=\min _{c \in\{r, g, b\}}\left(\min _{\mathbf{y} \in \Omega(\mathbf{x})}\left(J^{c}(\mathbf{y})\right)\right) \]
上述公式的意义用代码表达也很简单，首先求出每一个像素RGB份量中的最小值，存入一副和原始图像大小相同的灰度图中，而后再对这幅灰度图进行最小值滤波，滤波的半径由窗口大小决定。暗通道先验的理论指出：\(J^{\mathrm{dark}} \rightarrow 0\).

去雾公认的模型为：
\[ \mathbf{I}(\mathbf{x})=\mathbf{J}(\mathbf{x}) \mathrm{t}(\mathbf{x})+\mathbf{A}(1-\mathbf{t}(\mathbf{x})) \]
其中\(I\)为观测强度，\(J\)为场景亮度，\(A\)为全球大气光，\(t\)为描述非散射到达相机的光部分的介质透射率。去雾的目的是从\(I\)中恢复无雾的\(J\).

对上述公式进行变形，获得以下公式：
\[ \frac{I^{c}(\mathbf{x})}{A^{c}}=t(\mathbf{x}) \frac{J^{c}(\mathbf{x})}{A^{c}}+1-t(\mathbf{x}) \]
上标c表示RGB三个通道的意思，假设t在一个窗口下为常数，对上述公式两边同时取两次最小值，获得：
\[ \min _{\mathbf{y} \in \Omega(\mathbf{x})}\left(\min _{c} \frac{I^{c}(\mathbf{y})}{A^{c}}\right)=\tilde{t}(\mathbf{x}) \min _{\mathbf{y} \in \Omega(\mathbf{x})}\left(\min _{c} \frac{J^{c}(\mathbf{y})}{A^{c}}\right)+1-\tilde{t}(\mathbf{x}) \]
根据暗原色先验理论有：
\[ J^{\operatorname{dark}}(\mathbf{x})=\min _{\mathbf{y} \in \Omega(\mathbf{x})}\left(\min _{c} J^{c}(\mathbf{y})\right)=0 \]
因此前述公式能够化简为：
\[ \tilde{t}(\mathbf{x})=1-\min _{\mathbf{y} \in \Omega(\mathbf{x})}\left(\min _{c} \frac{I^{c}(\mathbf{y})}{A^{c}}\right) \]
在现实生活中，即便是晴天白云，空气中也存在着一些颗粒，所以，看远处的物体仍是能感受到雾的影响，另外，雾的存在让人类感到景深的存在，所以，有必要在去雾的时候保留必定程度的雾，这能够经过在上述式子中引入一个在[0,1] 之间的因子，则上式修正为：
\[ \tilde{t}(\mathbf{x})=1-\omega \min _{\mathbf{y} \in \Omega(\mathbf{x})}\left(\min _{c} \frac{I^{c}(\mathbf{y})}{A^{c}}\right) \]
论文中给出的\(\omega\)等于0.95，对A，论文中给出了一个计算方法：从暗通道图中按照亮度的大小取前0.1%的像素，在这些位置中，在原始有雾图像I中寻找对应的具备最高亮度的点的值，做为A值。这样A知道了，利用上述公式t也就知道了，在根据原始去雾公式，J也就能够计算了。公式为\(\mathrm{J}=(\mathrm{I}-\mathrm{A}) / \mathrm{t}+\mathrm{A}\)。

当t 的值很小时，会致使J的值偏大，从而使得图像总体向白场过分，所以通常可设置一阈值t0，当t值小于t0时，令t=t0，论文中取t0=0.1。
\[ \mathbf{J}(\mathbf{x})=\frac{\mathbf{I}(\mathbf{x})-\mathbf{A}}{\max \left(t(\mathbf{x}), t_{0}\right)}+\mathbf{A} \]

下面实现了一个最简单的版本，简化了不少，没有取窗口，没有用导向滤波等等，更多复杂的操做参考原始论文[13]。

实验结果：

source code:

import cv2
import argparse
import numpy as np


def hazeRemoval(img, w=0.7, t0=0.1):
    #求每一个像素的暗通道
    darkChannel = img.min(axis=2)
    #取暗通道的最大值最为全球大气光
    A = darkChannel.max()
    darkChannel = darkChannel.astype(np.double)
    #利用公式求得透射率
    t = 1 - w * (darkChannel / A)
    #设定透射率的最小值
    t[t < t0] = t0

    J = img
    #对每一个通道分别进行去雾
    J[:, :, 0] = (img[:, :, 0] - (1 - t) * A) / t
    J[:, :, 1] = (img[:, :, 1] - (1 - t) * A) / t
    J[:, :, 2] = (img[:, :, 2] - (1 - t) * A) / t
    return J

def main():
    ap = argparse.ArgumentParser()
    ap.add_argument('--image', required=True)
    args = vars(ap.parse_args())

    hazeImage = cv2.imread(args["image"])

    result = hazeRemoval(hazeImage.copy())

    cv2.imshow("HazeRemoval", np.hstack([hazeImage, result]))
    cv2.waitKey(0)


if __name__ == '__main__':
    main()

Reference

[1] http://homepages.inf.ed.ac.uk/rbf/HIPR2/hipr_top.htm

[2] https://blog.csdn.net/ytang2_/article/details/75451934

[3] https://baike.baidu.com/item/Sobel%E7%AE%97%E5%AD%90/11000092?fr=aladdin

[4] http://www.eie.polyu.edu.hk/~enyhchan/imagee.pdf

[5] http://www.eletel.p.lodz.pl/mstrzel/imageproc/enhancement1.PDF

[6] https://arxiv.org/ftp/arxiv/papers/1003/1003.4053.pdf

[7] https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=8076993

[8] https://blog.csdn.net/scottly1/article/details/42705271#commentBox

[9] https://zhuanlan.zhihu.com/p/44918476

[10] http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.294.8423&rep=rep1&type=pdf

[11] http://www.javashuo.com/article/p-gjretqkm-dr.html

[12] https://blog.csdn.net/baimafujinji/article/details/73824787#commentBox

[13] http://kaiminghe.com/publications/cvpr09.pdf

[14] http://www.cnblogs.com/Imageshop/p/3281703.html