[SciPy-User] Matlab imfilter/fspecial equivalent

[Brickle Macho]

I am porting some Matlab code to python. When I run the ported version
the output differs. It appears the difference may be due to how I have
translated Matlab's imfilter() and/or fspecial(). Below are my
translation, are there any Scipy/Matlab gurus that can let me know if
there are correct or if I have made some errors, or could suggest better
translations.

Matlab: imfilter(saliencyMap, fspecial('gaussian', [10, 10], 2.5))
Scipy: ndimage.gaussian_filter(saliencyMap, sigma=2.5)

Also the following. At the rsik of highlighting my lack of
understanding of the temrinology, Is uniform_filter the same as 'average'?

Matlab: imfilter(myLogAmplitude, fspecial('average', 3), 'replicate');
Scipy: scipy.ndimage.uniform_filter(mylogAmplitude, size=3, mode="nearest")


Thanks,

Michael.
--


If curious here are the 8 lines of matlab code I am trying to port:

%% Read image from file
inImg = im2double(rgb2gray(imread('yourImage.jpg')));
inImg = imresize(inImg, 64/size(inImg, 2));

%% Spectral Residual
myFFT = fft2(inImg);
myLogAmplitude = log(abs(myFFT));
myPhase = angle(myFFT);
mySpectralResidual = myLogAmplitude - imfilter(myLogAmplitude,
fspecial('average', 3), 'replicate');
saliencyMap = abs(ifft2(exp(mySpectralResidual + i*myPhase))).^2;

%% After Effect
saliencyMap = mat2gray(imfilter(saliencyMap, fspecial('gaussian', [10,
10], 2.5)));
imshow(saliencyMap);

[Travis Oliphant]

Hi Michael,

Thanks for sharing your code and example. The imfilter command is
equivalent to scipy.signal.correlate and scipy.ndimage.correlate (the one
in scipy.ndimage is faster I believe). These implement simple
correlation-based filtering given a finite kernel.

To get the same output you would need to generate the same kind of kernel
in Python as the Matlab fspecial command is producing. one approach to get
you started and to help separate the problem into 2 stages (reproducing the
imfilter and reproducing the fspecial filter) is to export the results of
the Matlab fspecial command and use that kernel in Python code (save it as
a .mat file and read that .mat file with Python).

SciPy does not have the equivalent to the fspecial command but you can
generate all kinds of 1-d special filters with scipy.signal.get_window.
You can also "generate" the filters you need directly from code.

Here is some untested code for generating something close to what
fspecial('gaussian", [10,10], 2.5) would be producing

import numpy as np

def fgaussian(size, sigma):
m,n = size
h, k = m//2, n//2
x, y = np.mgrid[-h:h, -k:k]

Post by Brickle Macho
I am porting some Matlab code to python. When I run the ported version
the output differs. It appears the difference may be due to how I have
translated Matlab's imfilter() and/or fspecial(). Below are my
translation, are there any Scipy/Matlab gurus that can let me know if there
are correct or if I have made some errors, or could suggest better
translations.
Matlab: imfilter(saliencyMap, fspecial('gaussian', [10, 10], 2.5))
Scipy: ndimage.gaussian_filter(saliencyMap, sigma=2.5)
Also the following. At the rsik of highlighting my lack of understanding
of the temrinology, Is uniform_filter the same as 'average'?
Matlab: imfilter(myLogAmplitude, fspecial('average', 3), 'replicate');
Scipy: scipy.ndimage.uniform_filter(mylogAmplitude, size=3,
mode="nearest")
Thanks,
Michael.
--
%% Read image from file
inImg = im2double(rgb2gray(imread('yourImage.jpg')));
inImg = imresize(inImg, 64/size(inImg, 2));
%% Spectral Residual
myFFT = fft2(inImg);
myLogAmplitude = log(abs(myFFT));
myPhase = angle(myFFT);
mySpectralResidual = myLogAmplitude - imfilter(myLogAmplitude,
fspecial('average', 3), 'replicate');
saliencyMap = abs(ifft2(exp(mySpectralResidual + i*myPhase))).^2;
%% After Effect
saliencyMap = mat2gray(imfilter(saliencyMap, fspecial('gaussian', [10,
10], 2.5)));
imshow(saliencyMap);
_______________________________________________
SciPy-User mailing list
http://mail.scipy.org/mailman/listinfo/scipy-user

[Travis Oliphant]

I accidentally hit send before I finished the original email (Gmail
interface quirk strikes again...) I'll including the entire original
email....

Hi Michael,

Thanks for sharing your code and example. The imfilter command is
equivalent to scipy.signal.correlate and scipy.ndimage.correlate (the one
in scipy.ndimage is faster I believe). These implement simple
correlation-based filtering given a finite kernel.

To get the same output you would need to generate the same kind of kernel
in Python as the Matlab fspecial command is producing. one approach to get
you started and to help separate the problem into 2 stages (reproducing the
imfilter and reproducing the fspecial filter) is to export the results of
the Matlab fspecial command and use that kernel in Python code (save it as
a .mat file and read that .mat file with Python).

SciPy does not have the equivalent to the fspecial command but you can
generate all kinds of 1-d special filters with scipy.signal.get_window.
You can also "generate" the filters you need directly from code.

Here is some untested code for generating something close to what
fspecial('gaussian", [10,10], 2.5) would be producing

import numpy as np

def fgaussian(size, sigma):
m,n = size
h, k = m//2, n//2
x, y = np.mgrid[-h:h, -k:k]
return np.exp(-(x**2 + y**2)/(2*sigma**2))

Up to a scaling constant (and possibly shifted a bit) this should return a
similar filter to fspecial('gaussian', size, sigma).

I'm not completely sure if the ndimage.gaussian_filter does a finite-window
convolution like imfilter (and correlate) or if it does something else.

The uniform filter is the same as an averaging filter (up to a scaling
constant). To other approach is to just use scipy.ndimage.correlate with
np.ones(size) / np.product(size) where size is the size of the kernel.

Perhaps you would be willing to post your code when you get it to work.

Best,

-Travis

Post by Brickle Macho
I am porting some Matlab code to python. When I run the ported version
the output differs. It appears the difference may be due to how I have
translated Matlab's imfilter() and/or fspecial(). Below are my
translation, are there any Scipy/Matlab gurus that can let me know if there
are correct or if I have made some errors, or could suggest better
translations.
Matlab: imfilter(saliencyMap, fspecial('gaussian', [10, 10], 2.5))
Scipy: ndimage.gaussian_filter(saliencyMap, sigma=2.5)
Also the following. At the rsik of highlighting my lack of understanding
of the temrinology, Is uniform_filter the same as 'average'?
Matlab: imfilter(myLogAmplitude, fspecial('average', 3), 'replicate');
Scipy: scipy.ndimage.uniform_filter(mylogAmplitude, size=3,
mode="nearest")
Thanks,
Michael.
--
%% Read image from file
inImg = im2double(rgb2gray(imread('yourImage.jpg')));
inImg = imresize(inImg, 64/size(inImg, 2));
%% Spectral Residual
myFFT = fft2(inImg);
myLogAmplitude = log(abs(myFFT));
myPhase = angle(myFFT);
mySpectralResidual = myLogAmplitude - imfilter(myLogAmplitude,
fspecial('average', 3), 'replicate');
saliencyMap = abs(ifft2(exp(mySpectralResidual + i*myPhase))).^2;
%% After Effect
saliencyMap = mat2gray(imfilter(saliencyMap, fspecial('gaussian', [10,
10], 2.5)));
imshow(saliencyMap);
_______________________________________________
SciPy-User mailing list
http://mail.scipy.org/mailman/listinfo/scipy-user

[Brickle Macho]

Hi Travis,

Thanks for the response, it has given a better understanding of the
functions. Once working, I will post the new version.

The code comes from [1], but [2] questions the use/need of log-amplitude
and experimentally show that just using the phase produces the same
result. This new approach reduces the computation by about 1/3. I
have been able to confirm that using just the phase component produces
the same result within either environment, Matlab or Scipy. Obviously
I still have the problem that the output from Matlab differs from Scipy.

I also discovered that if I don't resize the image (it is currently
being reduced to approx 64x64) the output from Matlab is almost the same
as Python. This suggests to me that, since keeping the kernels the same
and using a smaller/larger images implies the differences in output is
due to the kernel being used (I think). So for now I will focus my
effort on implementing the the fspecial('gaussian', [10,10], 2.5) using
the one provided below as a starting point.

Thanks again for the pointers.

Michael.
--
[1] Xiaodi Hou and Liqing Zhang , Saliency Detection: A Spectral
Residual Approach, Journal Computer Vision and Pattern Recognition,
IEEE Computer Society Conference 2007

[2] Guo, Chenlei and Ma, Qi and Zhang, Liming, Spatio-temporal saliency
detection using phase spectrum of quaternion fourier transform, Computer
Vision and Pattern Recognition, 2008. CVPR 2008.

Post by Travis Oliphant
I accidentally hit send before I finished the original email (Gmail
interface quirk strikes again...) I'll including the entire original
email....
Hi Michael,
Thanks for sharing your code and example. The imfilter command is
equivalent to scipy.signal.correlate and scipy.ndimage.correlate (the
one in scipy.ndimage is faster I believe). These implement simple
correlation-based filtering given a finite kernel.
To get the same output you would need to generate the same kind of
kernel in Python as the Matlab fspecial command is producing. one
approach to get you started and to help separate the problem into 2
stages (reproducing the imfilter and reproducing the fspecial filter)
is to export the results of the Matlab fspecial command and use that
kernel in Python code (save it as a .mat file and read that .mat file
with Python).
SciPy does not have the equivalent to the fspecial command but you can
generate all kinds of 1-d special filters with
scipy.signal.get_window. You can also "generate" the filters you
need directly from code.
Here is some untested code for generating something close to what
fspecial('gaussian", [10,10], 2.5) would be producing
import numpy as np
m,n = size
h, k = m//2, n//2
x, y = np.mgrid[-h:h, -k:k]
return np.exp(-(x**2 + y**2)/(2*sigma**2))
Up to a scaling constant (and possibly shifted a bit) this should
return a similar filter to fspecial('gaussian', size, sigma).
I'm not completely sure if the ndimage.gaussian_filter does a
finite-window convolution like imfilter (and correlate) or if it does
something else.
The uniform filter is the same as an averaging filter (up to a scaling
constant). To other approach is to just use scipy.ndimage.correlate
with np.ones(size) / np.product(size) where size is the size of the
kernel.
Perhaps you would be willing to post your code when you get it to work.
Best,
-Travis
I am porting some Matlab code to python. When I run the ported
version the output differs. It appears the difference may be due
to how I have translated Matlab's imfilter() and/or fspecial().
Below are my translation, are there any Scipy/Matlab gurus that
can let me know if there are correct or if I have made some
errors, or could suggest better translations.
Matlab: imfilter(saliencyMap, fspecial('gaussian', [10, 10], 2.5))
Scipy: ndimage.gaussian_filter(saliencyMap, sigma=2.5)
Also the following. At the rsik of highlighting my lack of
understanding of the temrinology, Is uniform_filter the same as 'average'?
Matlab: imfilter(myLogAmplitude, fspecial('average', 3),
'replicate');
Scipy: scipy.ndimage.uniform_filter(mylogAmplitude, size=3, mode="nearest")
Thanks,
Michael.
--
%% Read image from file
inImg = im2double(rgb2gray(imread('yourImage.jpg')));
inImg = imresize(inImg, 64/size(inImg, 2));
%% Spectral Residual
myFFT = fft2(inImg);
myLogAmplitude = log(abs(myFFT));
myPhase = angle(myFFT);
mySpectralResidual = myLogAmplitude - imfilter(myLogAmplitude,
fspecial('average', 3), 'replicate');
saliencyMap = abs(ifft2(exp(mySpectralResidual + i*myPhase))).^2;
%% After Effect
saliencyMap = mat2gray(imfilter(saliencyMap, fspecial('gaussian',
[10, 10], 2.5)));
imshow(saliencyMap);
_______________________________________________
SciPy-User mailing list
http://mail.scipy.org/mailman/listinfo/scipy-user
--
---
Travis Oliphant
Continuum Analytics, Inc.
http://www.continuum.io

[Brickle Macho]

I am getting closer. Here are some observations:

* If I don't imresize() the original image, then visually I can't tell
the difference in the output image from Matlab or Scipy.
* Empirically the results using different kernels are the same. That is
the difference between two images is essentially zero. That is doesn't
seem to matter if I use the imported the "fspecial()" kernel from
matlab; generate my own kernel using a fgaussion(), or use
scipy.ndimage.gaussian_filter().

So look like the problem is with ndimage.resize().

Michael.
--

Post by Travis Oliphant
I accidentally hit send before I finished the original email (Gmail
interface quirk strikes again...) I'll including the entire original
email....
Hi Michael,
Thanks for sharing your code and example. The imfilter command is
equivalent to scipy.signal.correlate and scipy.ndimage.correlate (the
one in scipy.ndimage is faster I believe). These implement simple
correlation-based filtering given a finite kernel.
To get the same output you would need to generate the same kind of
kernel in Python as the Matlab fspecial command is producing. one
approach to get you started and to help separate the problem into 2
stages (reproducing the imfilter and reproducing the fspecial filter)
is to export the results of the Matlab fspecial command and use that
kernel in Python code (save it as a .mat file and read that .mat file
with Python).
SciPy does not have the equivalent to the fspecial command but you can
generate all kinds of 1-d special filters with
scipy.signal.get_window. You can also "generate" the filters you
need directly from code.
Here is some untested code for generating something close to what
fspecial('gaussian", [10,10], 2.5) would be producing
import numpy as np
m,n = size
h, k = m//2, n//2
x, y = np.mgrid[-h:h, -k:k]
return np.exp(-(x**2 + y**2)/(2*sigma**2))
Up to a scaling constant (and possibly shifted a bit) this should
return a similar filter to fspecial('gaussian', size, sigma).
I'm not completely sure if the ndimage.gaussian_filter does a
finite-window convolution like imfilter (and correlate) or if it does
something else.
The uniform filter is the same as an averaging filter (up to a scaling
constant). To other approach is to just use scipy.ndimage.correlate
with np.ones(size) / np.product(size) where size is the size of the
kernel.
Perhaps you would be willing to post your code when you get it to work.
Best,
-Travis
I am porting some Matlab code to python. When I run the ported
version the output differs. It appears the difference may be due
to how I have translated Matlab's imfilter() and/or fspecial().
Below are my translation, are there any Scipy/Matlab gurus that
can let me know if there are correct or if I have made some
errors, or could suggest better translations.
Matlab: imfilter(saliencyMap, fspecial('gaussian', [10, 10], 2.5))
Scipy: ndimage.gaussian_filter(saliencyMap, sigma=2.5)
Also the following. At the rsik of highlighting my lack of
understanding of the temrinology, Is uniform_filter the same as 'average'?
Matlab: imfilter(myLogAmplitude, fspecial('average', 3),
'replicate');
Scipy: scipy.ndimage.uniform_filter(mylogAmplitude, size=3, mode="nearest")
Thanks,
Michael.
--
%% Read image from file
inImg = im2double(rgb2gray(imread('yourImage.jpg')));
inImg = imresize(inImg, 64/size(inImg, 2));
%% Spectral Residual
myFFT = fft2(inImg);
myLogAmplitude = log(abs(myFFT));
myPhase = angle(myFFT);
mySpectralResidual = myLogAmplitude - imfilter(myLogAmplitude,
fspecial('average', 3), 'replicate');
saliencyMap = abs(ifft2(exp(mySpectralResidual + i*myPhase))).^2;
%% After Effect
saliencyMap = mat2gray(imfilter(saliencyMap, fspecial('gaussian',
[10, 10], 2.5)));
imshow(saliencyMap);
_______________________________________________
SciPy-User mailing list
http://mail.scipy.org/mailman/listinfo/scipy-user
--
---
Travis Oliphant
Continuum Analytics, Inc.
http://www.continuum.io

[Michael Borck]

An update. On further testing, the problem appeared to be caused by
converting the RGB image to grayscale, not imresize(). I was getting
different image values using a rgb2grayscale convertor form skimage.
If I use ndimage.imread(, flatten=True) then I get the same values in
the image array/matrix. Since I was using skimage to do the
conversion, I posted a query on the skimage mail list, and got back the
following response:

They're just different color conversion factors. Based on
http://www.mathworks.com/help/images/ref/rgb2gray.html, Matlab uses:
0.2989 R + 0.5870 G + 0.1140 B

Based on the docstring for `color.rgb2gray`:
0.2125 R + 0.7154 G + 0.0721 B

So I think this will explain why the output is similar in structure but
different in detail. I will write a convertor using the Matlab's
conversion factors to convince myself that the output is the same as in
the paper.

I will post the final ported version later.

Michael.
--

Post by Brickle Macho
* If I don't imresize() the original image, then visually I can't tell
the difference in the output image from Matlab or Scipy.
* Empirically the results using different kernels are the same. That
is the difference between two images is essentially zero. That is
doesn't seem to matter if I use the imported the "fspecial()" kernel
from matlab; generate my own kernel using a fgaussion(), or use
scipy.ndimage.gaussian_filter().
So look like the problem is with ndimage.resize().
Michael.
--