I'm trying to implement the homomorphic filtering as implemented by this research article - step 4 (pdf)
The original Matlab code, which I'm trying to translate into Python, is as follows (WARNING: only the first case of the switch is relevant):
%%%%%%%%%%%%%%%%%%%%%%%
% 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
It is called as follows:
%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);
My problem is that I get an overflow error at the final exponential computation:
test_CV.py:64: RuntimeWarning: overflow encountered in exp
I don't understand why, since I think my Python code is equivalent to the Matlab original one. I don't have a pre-processing stage to adapt the image size as a power of 2 pixels square, still, I execute my script on a 512x512 image.
Here is my Python script:
#!/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()
All the code associated with the article is available here
Thanks !
EDIT & FIX:
Errors were indeed the lack of normalization and the usage of cv2.dft()
which returns a 3 dimensional array. To be equivalent to the computation of matlab fft2()
, I switched to numpy.fft()
which returns a 2 dimensional array as Matlab does.
The resizing and symmetrization don't seem absolutely necessary (the corrected script worked without it).
Still, it appears matlab rgb2ycbcr()
doesn't give the same Y component as python cv2.cvtColor(img, cv2.COLOR_BGR2YCrCb)
(slight difference in the resulting arrays). That impacts all the following processing and final result.
Furthermore, Matlab doesn't need to de-normalize the resulting image in order to show the latter, whereas I needed to result*255
to visualize it in my Python script. That has a tendency, along with the slight difference in the RGB->RCrCb transformation, to alter the "quality" of the filter translation into Python.
That being said, here is the functional python script, if it interests anyone:
#!/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)
result_interm
about the same size as the number causing problems in that question (1234.1
) ? if so, it seems like you may need to scale things to fit. $\endgroup$ – Peter K.♦ Jul 17 '17 at 14:48ori = double(ori)/255;
. I can't see something similar in the python version. $\endgroup$ – Peter K.♦ Jul 17 '17 at 15:39