1

EDIT: After some more testing and a response form the scipy mailing list, the issue appears to be with fspecial(). To get the same output I need to generate the same kind of kernel in Python as the Matlab fspecial command is producing. For now I will try to export the kernel from matlab and work from there. Added as a edit since question has been "closed"

I am trying to port the following MATLAB code to Python. It seems to work but the output is different form MATLAB. I think the problem is with apply a "mean" filter to the log(amplituide). Any help appreciated.

The MATLAB code is from: http://www.klab.caltech.edu/~xhou/projects/spectralResidual/spectralresidual.html

%% Read image from file
inImg = im2double(rgb2gray(imread('1.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);

Here is my attempt in python:

from skimage import img_as_float                                                       
from skimage.io import imread
from skimage.color import rgb2gray                                                    
from scipy import fftpack, ndimage, misc                                        
from scipy.ndimage import uniform_filter
from matplotlib.pyplot as plt

# Read image from file
image = img_as_float(rgb2gray(imread('1.jpg')))
image = misc.imresize(image, 64.0 / image.shape[0])

# Spectral Residual
fft = fftpack.fft2(image)                                                   
logAmplitude = np.log(np.abs(fft))                                          
phase = np.angle(fft)                                                       
avgLogAmp = uniform_filter(logAmplitude, size=3, mode="nearest") #Is this same a applying "mean" filter            
spectralResidual = logAmplitude - avgLogAmp                                 
saliencyMap = np.abs(fftpack.ifft2(np.exp(spectralResidual + 1j * phase))) ** 2

# After Effect
saliencyMap = ndimage.gaussian_filter(sm, sigma=2.5)
plt.imshow(sm)
plt.show()

For completness here is a input image and the output from MATLAB and python.

Input Image, Island Saliency Map (MATLAB) Saliency Map (PYTHON)

Michael
  • 559
  • 6
  • 15
  • Ummm... not sure how to improve the question. I am off to investigate differences between FFT implementations (Thanks Colin). Can't see much chance of it reopening. Annoying that I have to wait 2 days to delete? – Michael May 15 '13 at 22:35
  • 1
    After some more testing and a response form the scipy mailing list, the issue appears to be with fspecial(). To get the same output I need to generate the same kind of kernel in Python as the Matlab fspecial command is producing. For now I will try to export the kernel from Matlab and work from there. I will try to add this as a edit to the original post. – Michael May 16 '13 at 05:24

1 Answers1

2

I doubt anyone will be able to give you a firm answer on this. It could be any number of things... Could be that one FFT is 0-centered while the other isn't, could be a float vs double somewhere, could be mishandling of absolute value, could be a filter setting, ...

If I were you, I'd write out some intermediate values for both computations and find a way to compare them. Start in the middle, if they compare well then move down, if they don't compare well then move up. Maybe write an intermediate value from the python script out to a file, import into matlab, take the element-wise difference, and graph. If they're not the same dimensions, that's clue #1.

Colin Nichols
  • 714
  • 4
  • 12
  • Thanks, very helpful. I had assumed there was an error in the translation to python, didn't think there might be difference in the way the FFT is generated. I will do more research. – Michael May 15 '13 at 22:37
  • Found this SO post: http://stackoverflow.com/questions/8680909/fft-in-matlab-and-numpy-scipy-give-different-results which may explain the differences. – Michael May 15 '13 at 22:55
  • Okay. The transpose didn't help in my case. If I remove the imresize(), then the images appear the same, well almost the same. I now suspect the "real" problem/difference is being cause by the final gaussian filter applied, which has more noticble effect on a smaller image. Just need to find python equivalent to imfilter(saliencyMap, fspecial('gaussian', [10, 10], 2.5))" – Michael May 16 '13 at 01:50