圖像加強綜述

時間 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()