圖像增強綜述

本文轉載自查看原文 2019-04-19 11:57 11101 直方圖均衡/ 圖像相減/ 暗通道先驗/ Retinex/ 圖像反轉/ 圖像濾波/ 圖像增強

作者：方陽, 轉載請注明地址。
文件和代碼可以在Github下載, markdown推薦用typora打開。
這篇文章是DIP的第二次作業，對圖像增強技術進行綜述，目錄如下：

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）在許多應用中都很有用，例如顯示醫學圖像和用單色正片拍攝屏幕，其想法是將產生的負片用作投影片。

轉換方程：\(T: G(x, y) = L - F(x,y)\)，其中\(L\)是最大強度值，灰度圖像\(L\)為255。

實驗結果：

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）背后的想法是增加圖像處理中灰度的動態范圍。

轉換公式：\(G(x, y) = g_1 + (g_2 - g_1)/ (f_2 - f_1)[F(x, y) - 1]\)，其中這里\([f1, f2]\)為灰度在新范圍\([g1, g2]\)上的映射，這里f1為圖像的最小強度值，f2為圖像的最大強度值。該函數增強了圖像的對比度，顯示了均勻的強度分布。

實驗結果：

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）。

轉換公式：\(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算子/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] https://www.cnblogs.com/wangyong/p/8665434.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

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 X光圖像增強圖像增強 PyTorch【9】-圖像增強圖像增強之各種算子圖像增強處理-1 圖像增強之圖像去霧圖像增強-圖像銳化 python,opencv圖像增強低光圖像增強學習用python實現對圖像增強