【圖像處理】-012 同态濾波
在上一篇中，在實作底帽變換用于校正不均勻光照引起的變化時，發現使用同态濾波也可以達到同樣的效果，是以，對同态濾波進行了一些調研。
文章目錄
- 1 理論依據
-
- 1.1 圖像形成模型
- 1.2 同态濾波
- 2 實作
-
- 2.1 Matlab實作
- 2.2 OpenCV實作
- 3 效果
-
- 3.1 Matlab
- 3.2 OpenCV
1 理論依據
1.1 圖像形成模型
形如 f ( x , y ) f(x,y) f(x,y)的二維函數用于表示圖像，在空間坐标 ( x , y ) (x,y) (x,y)處， f f f的值或幅度是一個正的标量，其實體意義由圖像源決定。當一幅圖像由實體過程産生時，其亮度值正比于實體源（例如電磁波）所輻射的能量。是以， f ( x , y ) f(x,y) f(x,y)一定是非零的和有限的，即
(1) 0 < f ( x , y ) < ∞ 0 < f(x,y) < \infty \tag{1} 0<f(x,y)<∞(1)
函數 f ( x , y ) f(x,y) f(x,y)可以由兩個分量來表征：
- 入射分量：入射到被觀察場景的光源照射總量，用 i ( x , y ) i(x,y) i(x,y)表示；
- 反射分量：場景中物體所反射的光照總量，用 r ( x , y ) r(x,y) r(x,y)表示。
這兩個分量的乘積合并形成 f ( x , y ) f(x,y) f(x,y)，
(2) f ( x , y ) = i ( x , y ) r ( x , y ) f(x,y)=i(x,y)r(x,y) \tag{2} f(x,y)=i(x,y)r(x,y)(2)
其中，
(3) 0 < i ( x , y ) < ∞ 0 < i(x,y) < \infty \tag{3} 0<i(x,y)<∞(3)
(4) 0 < r ( x , y ) < 1 0 < r(x,y) < 1 \tag{4} 0<r(x,y)<1(4)
式 ( 4 ) (4) (4)指出反射分量顯示在0（全吸收）和1（全反射）之間。 i ( x , y ) i(x,y) i(x,y)的性質取決于反射源，而 r ( x , y ) r(x,y) r(x,y)的性質則取決于成像物體的特性。這種表示還可以用于照射光通過一個媒體形成圖像的情況。
1.2 同态濾波
上面談到的圖像形成模型可用于開發一種頻率域處理過程，該過程通過同時壓縮灰階範圍和增強對比度來改善一幅圖像的表觀。一幅圖像 f ( x , y ) f(x,y) f(x,y)可以表示為其照射分量和反射分量的乘積，即
(5) f ( x , y ) = i ( x , y ) r ( x , y ) f(x,y)=i(x,y)r(x,y) \tag{5} f(x,y)=i(x,y)r(x,y)(5)
上式不能直接用于對照射和反射分量的操作，因為乘積的傅立葉變換不是傅立葉變換的乘積。
(6) F [ f ( x , y ) ] ̸ = F [ i ( x , y ) ] F [ r ( x , y ) ] F[f(x,y)] \not= F[i(x,y)]F[r(x,y)] \tag{6} F[f(x,y)]̸=F[i(x,y)]F[r(x,y)](6)
由數學可以，乘積的對數等于對數的和，而和的傅立葉積分等于傅立葉積分的和，是以，可以定義：
(7) z ( x , y ) = l n ( f ( x , y ) ) = l n ( i ( x , y ) ) + l n ( r ( x , y ) ) z(x,y)=ln(f(x,y))=ln (i(x,y)) + ln (r(x,y)) \tag{7} z(x,y)=ln(f(x,y))=ln(i(x,y))+ln(r(x,y))(7)
那麼，
(8) F [ z ( x , y ) ] = F [ l n ( f ( x , y ) ) ] = F [ l n ( i ( x , y ) ) ] + F [ l n ( r ( x , y ) ) ] F[z(x,y)]=F[ln(f(x,y))]=F[ln(i(x,y))]+F[ln(r(x,y))] \tag{8} F[z(x,y)]=F[ln(f(x,y))]=F[ln(i(x,y))]+F[ln(r(x,y))](8)
或
(9) Z ( u , v ) = F i ( u , v ) + F r ( u , v ) Z(u,v)=F_{i}(u,v)+F_{r}(u,v) \tag{9} Z(u,v)=Fi(u,v)+Fr(u,v)(9)
其中， F i ( u , v ) F_{i}(u,v) Fi(u,v)和 F r ( u , v ) F_{r}(u,v) Fr(u,v)分别表示 l n ( i ( x , y ) ) ln(i(x,y)) ln(i(x,y))和 l n ( r ( x , y ) ) ln(r(x,y)) ln(r(x,y))的傅立葉變換。
用濾波器 H ( u , v ) H(u,v) H(u,v)對 Z ( u , v ) Z(u,v) Z(u,v)進行濾波，有：
(10) S ( u , v ) = H ( u , v ) Z ( u , v ) = H ( u , v ) F i ( u , v ) + H ( u , v ) F r ( u , v ) S(u,v)=H(u,v)Z(u,v)=H(u,v)F_{i}(u,v)+H(u,v)F_{r}(u,v) \tag{10} S(u,v)=H(u,v)Z(u,v)=H(u,v)Fi(u,v)+H(u,v)Fr(u,v)(10)
在空間域中，濾波後的圖像是：
(11) s ( x , y ) = F − 1 { S ( u , v ) } = F − 1 { H ( u , v ) F i ( u , v ) } + F − 1 { H ( u , v ) F r ( u , v ) } s(x,y)=F^{-1}\{ S(u,v)\}=F^{-1}\{ H(u,v)F_{i}(u,v)\}+F^{-1}\{ H(u,v)F_{r}(u,v)\} \tag{11} s(x,y)=F−1{S(u,v)}=F−1{H(u,v)Fi(u,v)}+F−1{H(u,v)Fr(u,v)}(11)
由定義
(12) i ′ ( x , y ) = F − 1 { H ( u , v ) F i ( u , v ) } i^{'}(x,y)=F^{-1}\{ H(u,v)F_{i}(u,v)\} \tag{12} i′(x,y)=F−1{H(u,v)Fi(u,v)}(12)
(13) r ′ ( x , y ) = F − 1 { H ( u , v ) F r ( u , v ) } r^{'}(x,y)=F^{-1}\{ H(u,v)F_{r}(u,v)\} \tag{13} r′(x,y)=F−1{H(u,v)Fr(u,v)}(13)
是以：
(14) s ( x , y ) = i ′ ( x , y ) + r ′ ( x , y ) s(x,y)=i^{'}(x,y)+r^{'}(x,y) \tag{14} s(x,y)=i′(x,y)+r′(x,y)(14)
因為 z ( x , y ) z(x,y) z(x,y)是通過輸入圖像的自然對數形成的，是以可以通過對濾波後的圖像的指數這一反處理過程來形成輸出圖像。
(15) g ( x , y ) = e s ( x , y ) = e i ′ ( x , y ) e r ′ ( x , y ) = i 0 ( x , y ) r 0 ( x , y ) g(x,y) = e^{s(x,y)}=e^{i^{'}(x,y)}e^{r^{'}(x,y)}=i_{0}(x,y)r_{0}(x,y) \tag{15} g(x,y)=es(x,y)=ei′(x,y)er′(x,y)=i0(x,y)r0(x,y)(15)
其中，
(16) i 0 ( x , y ) = e i ′ ( x , y ) i_{0}(x,y)=e^{i^{'}(x,y)} \tag{16} i0(x,y)=ei′(x,y)(16)
(17) r 0 ( x , y ) = e r ′ ( x , y ) r_{0}(x,y)=e^{r^{'}(x,y)} \tag{17} r0(x,y)=er′(x,y)(17)
是輸出圖像的照射和反射分量。
該方法是以稱之為同态系統的一類系統的特殊情況為基礎的，方法的關鍵是照射分量和反射分量的分離，然後，同态濾波函數 H ( u , v ) H(u,v) H(u,v)可以對各分量進行操作。
由上圖可以看出同态濾波的過程如下：
- (1) 對輸入圖像進行對數變換，原理上對應将照射分量和反射分量分離；
- (2) 進行傅立葉變換；
- (3) 在頻域對圖像進行濾波處理；
- (4) 對濾波處理的結果進行反傅立葉變換；
- (5) 對反傅立葉變換的結果進行指數變換，得到輸出結果圖像。
圖像的照射分量通常由慢的空間變化來表征，而反射分量往往引起突變，特别是在不同物體連接配接不部分。這些特征導緻圖像取對數之後的傅立葉變換的低頻成分與照射相關聯，而高頻成分與反射相聯系。
使用同态濾波器可以更好地控制照射分量和反射分量。這種控制需要指定一個濾波器函數 H ( u , v ) H(u,v) H(u,v),它可以對傅立葉變換的低頻和高頻部分用不同的方法處理。
上圖顯示的是一個可用的濾波器剖面圖。如果 γ H \gamma _{H} γH和 γ L \gamma _{L} γL標明，而 γ H > 1 \gamma _{H}>1 γH>1且 γ L < 1 \gamma _{L}<1 γL<1，那麼該濾波器函數取象于衰減低頻（照射）分量，而增強高頻（反射）分量。最終結果是同僚進行動态範圍的壓縮和對比度的增強。
(18) H ( u , v ) = ( γ H − γ L ) [ 1 − e − c [ D 2 ( u , v ) / D 0 2 ] ] + γ L H(u,v)=(\gamma_{H}-\gamma_{L})[1-e^{-c[D^{2}(u,v)/D_{0}^2]}]+\gamma_{L} \tag{18} H(u,v)=(γH−γL)[1−e−c[D2(u,v)/D02]]+γL(18)
2 實作
2.1 Matlab實作
function [image_out] = myHF(image_in, rh, rl, c, D0)
% 同态濾波器
% 輸入為需要進行濾波的灰階圖像，同态濾波器的參數rh, rl,c, D0
% 輸出為進行濾波之後的灰階圖像
[m, n] = size(image_in);
P = 2*m;
Q = 2*n;
% 取對數
image1 = log(double(image_in) + 1);
fp = zeros(P, Q);
%對圖像填充0,并且乘以(-1)^(x+y) 以移到變換中心
for i = 1 : m
for j = 1 : n
fp(i, j) = double(image1(i, j)) * (-1)^(i+j);
end
end
% 對填充後的圖像進行傅裡葉變換
F1 = fft2(fp);
% 生成同态濾波函數，中心在(m+1,n+1)
Homo = zeros(P, Q);
a = D0^2;
% 計算一些不變的中間參數
r = rh-rl;
for u = 1 : P
for v = 1 : Q
temp = (u-(m+1.0))^2 + (v-(n+1.0))^2;
Homo(u, v) = r * (1-exp((-c)*(temp/a))) + rl;
end
end
%進行濾波
G = F1 .* Homo;
% 反傅裡葉變換
gp = ifft2(G);
% 處理得到的圖像
image_out = zeros(m, n, 'uint8');
gp = real(gp);
g = zeros(m, n);
for i = 1 : m
for j = 1 : n
g(i, j) = gp(i, j) * (-1)^(i+j);
end
end
% 指數處理
ge = exp(g)-1;
% 歸一化到[0, L-1]
mmax = max(ge(:));
mmin = min(ge(:));
range = mmax-mmin;
for i = 1 : m
for j = 1 : n
image_out(i,j) = uint8(255 * (ge(i, j)-mmin) / range);
end
end
end
調用方法：
clear all;
close all;
clc;
image1 = imread('../images/61.jpg');
image1 = rgb2gray(image1);
[m, n] = size(image1);
image2 = myHF(image1, 0.6, 0.55, 1, 80);
% 顯示圖像
subplot(1,2,1), imshow(image1), title('原圖像');
subplot(1,2,2), imshow(image2), title('D0 = 80');
2.2 OpenCV實作
#include "../include/importOpenCV.h"
#include "../include/baseOps.h"
#include "../include/opencv400/opencv2/core.hpp"
#include <iostream>
/*
//High-Frequency-Emphasis Filters
H(u,v)=(\gamma_{H}-\gamma_{L})[1-e^{-c[D^{2}(u,v)/D_{0}^2]}]+\gamma_{L}
根據輸入的濾波器尺寸，半徑D0，調節系數c，上界rh,下界rl，計算H(u,v)表示的高通濾波器，
傳回複數形式的濾波器，通過realIm傳回該濾波器的實部。
注意，傳回的複數濾波器的虛部為0。
*/
cv::Mat Butterworth_Homomorphic_Filter(cv::Size sz, int D, int c, float rh, float rl, cv::Mat& realIm)
{
cv::Mat single(sz.height, sz.width, CV_32F);
cv::Point centre = cv::Point(sz.height/2, sz.width/2);
double radius;
float upper = rh ;
float lower = rl ;
long dpow = D*D;
float W = (upper - lower);
for(int i = 0; i < sz.height; i++)
{
for(int j = 0; j < sz.width; j++)
{
radius = pow((float)(i - centre.x), 2) + pow((float) (j - centre.y), 2);
float r = exp(-c*radius/dpow);
if(radius < 0)
single.at<float>(i,j) = upper;
else
single.at<float>(i,j) =W*(1 - r) + lower;
}
}
single.copyTo(realIm);
//将實數濾波器擴充成虛部為0的複數濾波器
cv::Mat butterworth_complex;
//make two channels to match complex
cv::Mat butterworth_channels[] = { cv::Mat_<float>(single), cv::Mat::zeros(sz, CV_32F)};
cv::merge(butterworth_channels, 2, butterworth_complex);
return butterworth_complex;
}
/*
//DFT 傳回功率譜Power
計算輸入圖像的傅立葉變換，
傳回該圖像的傅立葉變換的功率譜；
通過image_comples傳回該圖像的複數形式的傅立葉變換值；
通過image_phase傳回該圖像傅立葉變換對應的相位；
通過image_mag傳回該圖像傅立葉變換的幅度譜
注意：這裡計算的不是标準的傅立葉變換，因為為了進行同态濾波，對輸入圖像進行了src = log(1+src)的操作。
*/
cv::Mat FFT2(cv::Mat frame_bw, cv::Mat& image_complex, cv::Mat &image_phase, cv::Mat &image_mag)
{
cv::Mat frame_log;
//将8UC1轉換成32FC1；
frame_bw.convertTo(frame_log, CV_32F);
frame_log = frame_log/255;//歸一化
//對輸入圖像進行自然對數變換，ln(src+1)
frame_log += 1;
cv::log( frame_log, frame_log);
/*2.Expand the image to an optimal size The performance of the DFT depends of the image size.
It tends to be the fastest for image sizes that are multiple of 2, 3 or 5.
We can use the copyMakeBorder() function to expand the borders of an image.*/
cv::Mat padded;
int M = cv::getOptimalDFTSize(frame_log.rows);
int N = cv::getOptimalDFTSize(frame_log.cols);
cv::copyMakeBorder(frame_log, padded, 0, M - frame_log.rows, 0, N - frame_log.cols, cv::BORDER_CONSTANT, cv::Scalar::all(0));
/*Make place for both the complex and real values The result of the DFT is a complex.
Then the result is 2 images (Imaginary + Real), and the frequency domains range is much
larger than the spatial one. Therefore we need to store in float ! That's why we will
convert our input image "padded" to float and expand it to another channel to hold
the complex values. Planes is an arrow of 2 matrix (planes[0] = Real part, planes[1] = Imaginary part)*/
cv::Mat image_planes[] = { cv::Mat_<float>(padded), cv::Mat::zeros(padded.size(), CV_32F)};
/*Creates one multichannel array out of several single-channel ones.*/
cv::merge(image_planes, 2, image_complex);
/*Make the DFT The result of thee DFT is a complex image : "image_complex"*/
cv::dft(image_complex, image_complex);
/***************************/
//Create spectrum magnitude//
/***************************/
/*Transform the real and complex values to magnitude NB: We separe Real part to Imaginary part*/
cv::split(image_complex, image_planes);
//Starting with this part we have the real part of the image in planes[0] and the imaginary in planes[1]
cv::phase(image_planes[0], image_planes[1], image_phase);
cv::magnitude(image_planes[0], image_planes[1], image_mag);
//Power
cv::pow(image_planes[0],2,image_planes[0]);
cv::pow(image_planes[1],2,image_planes[1]);
cv::Mat Power = image_planes[0] + image_planes[1];
return Power;
}
/*
傅立葉逆變換
注意，這裡為了進行同态濾波，對逆變換的結果進行了e為底的指數變換之後，
對結果進行了-1之後歸一化的操作，對應在傅立葉正變換的過程中的ln(src+1)操作。
*/
void IFFT2(cv::Mat input, cv::Mat& inverseTransform)
{
/*Make the IDFT*/
cv::Mat result;
cv::idft(input, result, cv::DFT_SCALE);
/*Take the exponential*/
cv::exp(result, result);
std::vector<cv::Mat> planes;
cv::split(result, planes);
cv::magnitude(planes[0], planes[1], planes[0]);
planes[0] = planes[0] - 1.0;
cv::normalize(planes[0], planes[0], 0, 255, cv::NORM_MINMAX);
planes[0].convertTo(inverseTransform, CV_8U);
}
//傅立葉變換時交換四個象限
void ShiftFFT2(cv::Mat &fImage)
{
//For visualization purposes we may also rearrange the quadrants of the result, so that the origin (0,0), corresponds to the image center.
cv::Mat tmp, q0, q1, q2, q3;
/*First crop the image, if it has an odd number of rows or columns.
Operator & bit to bit by -2 (two's complement : -2 = 111111111....10)
to eliminate the first bit 2^0 (In case of odd number on row or col,
we take the even number in below)*/
fImage = fImage(cv::Rect(0, 0, fImage.cols & -2, fImage.rows & -2));
int cx = fImage.cols/2;
int cy = fImage.rows/2;
/*Rearrange the quadrants of Fourier image so that the origin is at the image center*/
q0 = fImage(cv::Rect(0, 0, cx, cy));
q1 = fImage(cv::Rect(cx, 0, cx, cy));
q2 = fImage(cv::Rect(0, cy, cx, cy));
q3 = fImage(cv::Rect(cx, cy, cx, cy));
/*We reverse each quadrant of the frame with
its other quadrant diagonally opposite*/
/*We reverse q0 and q3*/
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
/*We reverse q1 and q2*/
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
}
int main()
{
//将工作目錄設定到EXE所在的目錄。
SetCurrentDirectoryToExePath();
cv::Mat src = cv::imread("../images/61.jpg");
cv::cvtColor(src, src, cv::COLOR_BGR2GRAY);
cv::imshow("原圖", src);
cv::Mat img;
cv::Mat imgHls;
std::vector<cv::Mat> vHls;
if (src.channels() == 3)
{
cvtColor(src, imgHls, cv::COLOR_BGR2HSV);
split(imgHls, vHls);
vHls[2].copyTo(img);
}
else
src.copyTo(img);
cv::Mat img_complex,img_mag,img_phase;
cv::Mat fpower = FFT2(img,img_complex,img_phase,img_mag);
ShiftFFT2(img_complex);
ShiftFFT2(fpower);
int D0 = 80;
int n = 1;
int w = img_complex.cols;
int h = img_complex.rows;
float rH = 0.6;
float rL =0.55;
cv::Mat filter,filter_complex;
filter_complex = Butterworth_Homomorphic_Filter(cv::Size(w,h),D0,n,rH,rL,filter);
cv::mulSpectrums(img_complex, filter_complex, filter_complex, 0);
ShiftFFT2(filter_complex);
cv::Mat result;
IFFT2(filter_complex,result);
cv::Mat dst;
if(src.channels() == 3)
{
vHls.at(2)= result(cv::Rect(0,0,src.cols,src.rows));
merge(vHls,imgHls);
cvtColor(imgHls, dst, cv::COLOR_HSV2BGR);
}
else
result.copyTo(dst);
cv::imshow("dst", dst);
cv::waitKey();
cv::destroyAllWindows();
return 0;
}
3 效果
為了測試算法的效果，我在網上找了些圖檔進行測試。測試圖檔的版權歸原作者所有，本人不主張對測試圖檔及結果的任何版權。圖檔來源于百度圖檔搜尋。
3.1 Matlab
3.2 OpenCV
可以看出，處理結果的尺寸相比原圖發生了變化，因為進行DFT時對圖像進行了填充，我暫時沒有注意去将結果的邊緣取消。