我正在尝试实现本研究文章所实现的同态过滤 - 步骤 4 (pdf)
我试图翻译成 Python 的原始 Matlab 代码如下（警告：只有开关的第一种情况是相关的）：
%%%%%%%%%%%%%%%%%%%%%%%
% Homomorphic filtering
% Stéphane BAZEILLE
% 15/10/07
%%%%%%%%%%
function result=homomorphic_filtering(I,type,RH,RL,cutoff)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% I = image (en niveau de gris et double entre 0 et 1) %%
%% type= 'HighPassLiao' 'HighPassCufi' 'HighPassFilter' %%
%% 'LowPassFilter' 'BandPassFilter' 'OnesFilter' %%
%% !!!!!!!!! afficher avec imshow sans contrast_stretching !!!!!!!!! %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[rows,cols,dim]=size(I);
%figure;imshow(I);title('homo in');
% figure;hist(I,0.01:0.01:1.49);
disp(['Homomorphic filtering (',type,',',num2str(RH),',',num2str(RL),')']);
%Modélisation du filtre
if mod(cols,2) xrange = [-(cols-1)/2:(cols-1)/2]/(cols-1);
else xrange = [-cols/2:(cols/2-1)]/cols;
end
if mod(rows,2) yrange = [-(rows-1)/2:(rows-1)/2]/(rows-1);
else yrange = [-rows/2:(rows/2-1)]/rows;
end
[x,y] = meshgrid(xrange, yrange);
radius = sqrt(x.^2 + y.^2);
switch type
case 'HighPassLiao'
DX=cols/cutoff;
G = ones(rows,cols);
for m = 1:rows
for n = 1:cols
G (m,n) = ((RH-RL)*(1-exp(-(((m-rows/2))*((m-rows/2))+((n-cols/2))*((n-cols/2)))/(2*DX*DX))))+RL;
%G (m,n) = ((RH-RL)*(1-exp(-(((m-rows/8))*((m-rows/8))+((n-cols/8))*((n-cols/8)))/(2*DX*DX))))+RL;
end
end
case 'HighPassCufi'
cutoff = 0.4; % 0 - 0.5
offset = 0.5;
boost = 2;
scale = 50;
G = boost./(1.0 + exp(-scale*(radius-cutoff))) + offset;
case 'HighPassFilter'
cutoff = 0.3; % 0 - 0.5
n = 3; % n >= 1
boost = 0.2; % 0<= boost=<1
G = (1-boost)*(1 - (1 ./ (1.0 + (radius ./ cutoff).^(2*n))))+boost;
case 'LowPassFilter'
cutoff = 0.1; % 0 - 0.5
n = 3; % n >= 1
boost = 0.2; % 0<= boost <1
G = (1-boost) ./ (1.0 + (radius ./ cutoff).^(2*n))+boost;
case 'BandPassFilter'
cutoff = 0.1; % 0 - 0.5
cutin = 0.01; % 0 - 0.5
n = 3; % n >= 1
G = (1 ./ (1.0 + (radius ./ cutoff).^(2*n)))-(1 ./ (1.0 + (radius ./ cutin).^(2*n)));
case 'OnesFilter'
G=ones(rows,cols);
end
%Filtrage homomorphique
if(dim==1)
Idft=fft2(log(I+0.01));
filtree=G.*fftshift(Idft);
result=exp(real(ifft2(ifftshift(filtree))));
else
Idft=fft2(log(I(:,:,1)+0.01));
filtree=G.*fftshift(Idft);
result(:,:,1)=exp(real(ifft2(ifftshift(filtree))));
Idft=fft2(log(I(:,:,2)+0.01));
filtree=G.*fftshift(Idft);
result(:,:,2)=exp(real(ifft2(ifftshift(filtree))));
Idft=fft2(log(I(:,:,3)+0.01));
filtree=G.*fftshift(Idft);
result(:,:,3)=exp(real(ifft2(ifftshift(filtree))));
end
end
它被称为如下：
%param: full image name
function prefilter_demo(name)
ori = imread(name);
ori = double(ori)/255;
disp('Image Loaded');
figure;subplot(1,2,1);imshow(ori);axis off
tic;
disp('BEGIN FILTERING');
% EXTEND AND RESIZE IMAGE
[nl nc I ori]=picture_power_of_2(ori,5);
% CONVERT RGB TO YCBCR
ycbcr = rgb2ycbcr(I);
% HOMOMORPHIC FILTERING
ycbcr(:,:,1)=homomorphic_filtering(ycbcr(:,:,1),'HighPassLiao',2.5,0.5,32);
我的问题是在最终的指数计算中出现溢出错误：
test_CV.py:64: RuntimeWarning: overflow encountered in exp
我不明白为什么，因为我认为我的 Python 代码等同于 Matlab 原始代码。我没有预处理阶段来将图像大小调整为 2 像素平方的幂，但我仍然在 512x512 图像上执行我的脚本。
这是我的 Python 脚本：
#!/usr/bin/env python3
# coding: utf-8
import cv2
import numpy as np
img = cv2.imread('bleu_carree_512.jpg',-1)
rows,cols,dim=img.shape
rh, rl, cutoff = 2.5,0.5,32
imgYCrCb = cv2.cvtColor(img, cv2.COLOR_BGR2YCrCb)
y,cr,cb = cv2.split(imgYCrCb)
y_log = np.log(y+0.01)
y_fft = cv2.dft(np.float32(y_log),flags = cv2.DFT_COMPLEX_OUTPUT)
y_fft_shift = np.fft.fftshift(y_fft)
DX = cols/cutoff
G = np.ones((rows,cols,2))
for i in range(rows):
for j in range(cols):
G[i][j][0]=((rh-rl)*(1-np.exp(-((i-rows/2)**2+(j-cols/2)**2)/(2*DX**2))))+rl
G[i][j][1]=((rh-rl)*(1-np.exp(-((i-rows/2)**2+(j-cols/2)**2)/(2*DX**2))))+rl
result_filter = np.multiply(y_fft_shift,G)
result_interm = np.real(cv2.idft(np.fft.ifftshift(result_filter)))
#################################
### OVERFLOW SOURCE BELOW #######
#################################
result = np.exp(result_interm)
print(result.shape)
cv2.imshow('résultat',result)
cv2.waitKey(0)
cv2.destroyAllWindows()
与文章相关的所有代码都可以在这里找到
谢谢 ！
编辑和修复：
错误确实是缺乏归一化，并且使用归一化cv2.dft()会返回一个 3 维数组。为了等效于 的计算matlab fft2()，我切换到numpy.fft()which 返回一个二维数组，就像 Matlab 所做的那样。
调整大小和对称化似乎不是绝对必要的（更正的脚本在没有它的情况下工作）。
尽管如此，它似乎matlab rgb2ycbcr()没有给出与python cv2.cvtColor(img, cv2.COLOR_BGR2YCrCb)（结果数组略有不同）相同的 Y 分量。这会影响所有后续处理和最终结果。
此外，Matlab 不需要对生成的图像进行反规范化以显示后者，而我需要result*255在我的 Python 脚本中将其可视化。这有一种趋势，以及 RGB->RCrCb 转换的细微差异，将过滤器转换为 Python 的“质量”改变。
话虽如此，如果有人感兴趣，这里是功能性 Python 脚本：
#!/usr/bin/python2.7
# coding: utf-8
import cv2
import numpy as np
import matplotlib.pyplot as plt
img = cv2.imread('U352.jpg',-1)
img = np.float32(img)
img = img/255
rows,cols,dim=img.shape
rh, rl, cutoff = 2.5,0.5,32
imgYCrCb = cv2.cvtColor(img, cv2.COLOR_BGR2YCrCb)
y,cr,cb = cv2.split(imgYCrCb)
y_log = np.log(y+0.01)
y_fft = np.fft.fft2(y_log)
y_fft_shift = np.fft.fftshift(y_fft)
DX = cols/cutoff
G = np.ones((rows,cols))
for i in range(rows):
for j in range(cols):
G[i][j]=((rh-rl)*(1-np.exp(-((i-rows/2)**2+(j-cols/2)**2)/(2*DX**2))))+rl
result_filter = G * y_fft_shift
result_interm = np.real(np.fft.ifft2(np.fft.ifftshift(result_filter)))
result = np.exp(result_interm)