homomorphic filter - python overflow

Question

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)

Answer 1

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)