# homomorphic filter - python overflow

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);
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 = double(ori)/255;
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

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 = 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)

• Try using 128 bit precision as suggested here. This seems more like a SO question than a DSP question. – Peter K. Jul 17 '17 at 14:25
• I don't know why, but the numpy.longdouble "data type = 13 is not supported" on my system (Debian Stretch 64). Furthermore, according to your link, my exponential computation would overflow with float128. There must be an algorithmic mistake, I'm going to compare middle results with the original Matlab code. Thanks anyway ! (and sorry if off-topic) – Blupon Jul 17 '17 at 14:41
• Are any of the entries in 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. – Peter K. Jul 17 '17 at 14:48
• there are of magnitude e+05, and I've tested my script with the test image of the original Matlab algorithm (with original preprocessing for resizing and extending symmetrically - cf. article), I still have those abnormal big values for result_interm (abnormal because way bigger than the max of the associated matrix in the original code). That surely means I've made an actual mistake in my translation of the algorithm/forgotten one entire step. – Blupon Jul 17 '17 at 15:10
• The original code appears to normalize the image with the line ori = double(ori)/255;. I can't see something similar in the python version. – Peter K. Jul 17 '17 at 15:39

# 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 = 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)