0

I am trying to replicate Matlab's imgaussfilter behaviour in Python but I have been unable to reproduce the results. The docs don't help much since there is not explanation of what exactly is being done (for instance how is that function different from fspecial or how it differs from Octave's imsmooth (with Gaussian as an argument).

Matlab Code

imgaussfilt(image,sigma)

with output

 [[-0.02936392 -0.03168419 -0.0343706  ...  0.03136455  0.02864487
   0.02585145]
 [-0.03212093 -0.0347433  -0.03775943 ...  0.03507484  0.03209807
   0.02906134]
 [-0.03512981 -0.03808864 -0.04147075 ...  0.03873834  0.03549163
   0.03219311]
 ...
 [-0.07713804 -0.08337475 -0.08975262 ... -0.04314206 -0.03945251
  -0.03606256]
 [-0.07145457 -0.07714807 -0.08297654 ... -0.03986635 -0.03641579
  -0.03323718]
 [-0.06605107 -0.07122191 -0.07651892 ... -0.03684535 -0.03362917
  -0.03065128]]

The closest I have managed to get with Python has been:

Python

skimage.filters.gaussian(image, sigma=s,mode = 'nearest',truncate=2.0)

with output

[[-0.02936397 -0.03168425 -0.03437067 ...  0.03136461  0.02864492
   0.0258515 ]
 [-0.03212099 -0.03474336 -0.03775951 ...  0.03507491  0.03209814
   0.0290614 ]
 [-0.03512988 -0.03808872 -0.04147083 ...  0.03873842  0.03549169
   0.03219317]
 ...
 [-0.07713819 -0.08337491 -0.08975279 ... -0.04314214 -0.03945259
  -0.03606263]
 [-0.07145471 -0.07714821 -0.0829767  ... -0.03986643 -0.03641586
  -0.03323725]
 [-0.06605119 -0.07122205 -0.07651907 ... -0.03684542 -0.03362923
  -0.03065134]]

which is similar but not exactly the same result. Is there a method that approximates it better? Do I need to change one of the parameters?

EDIT: If you are looking for a close approximation using OpenCV @avisionx has provided a good starting point in the solutions.

Jorge Diaz
  • 491
  • 5
  • 11

3 Answers3

2

There's multiple possible causes for this. Since Matlab is not open source software, it's impossible to know the exact cause, but we can take some educated guesses.

Regarding what Gaussian filtering is doing: it replaces every pixel with a weighted sum of its neighbors. The weighting is determined by the Gaussian probability density function, and the exact weight is determined by sigma but also by the truncate parameter that you can see in the scipy.ndimage.gaussian_filter documentation. This parameter is needed because, technically, the Gaussian function never decays to zero, it just approaches zero as the distance goes to infinity, so a "perfect" Gaussian blur would need to sum an infinite number of neighboring pixels, which is impossible. Therefore, the function needs to decide how far away you are close enough to zero to stop summing things.

An additional complicating factor is that Gaussian blur is separable, which means that you can do a 2D (or nD) blur by blurring separately along each axis in sequence, which can reduce computing costs. Separability affects the order of computation, and floating point computation is not exact, so if Matlab and Python are doing the operations in a different order, you would expect slightly different results.

Even if you filter the axes in the exact same order (you might achieve this by e.g. transposing the image in Python, filtering, then transposing it back), the underlying, low-level array computation libraries might be summing arrays in a different order, which again would change the exact results.

In short, some possible causes for the differences you're seeing:

  • different truncation
  • different implementations (separating the filter or not)
  • different order of processing of the axes
  • different ordering of the low-level sums

There may be no way of disentangling these, and in the end my advice mirrors Cris Luengo's comment: your results should probably not rely on values more precise than six significant digits. An exact match between Matlab and Python here will be meaningless, because neither is able to guarantee perfect accuracy with respect to the theoretical value of a Gaussian filter. Both are approximations and an exact match only means that you are making the same approximation errors in both.

Juan
  • 5,433
  • 21
  • 23
  • 1
    Nice summary. One more thing to add: MATLAB uses either a separable or a Fourier-domain implementation depending on `sigma`. There are also recursive filters, but I don’t know if either of OP’s methods use those. – Cris Luengo Oct 18 '19 at 02:10
  • Oh yeah, I totally forgot about Fourier convolution! ‍♂️ – Juan Oct 18 '19 at 02:34
  • Thank you very much for your answer Juan. I understand that trying to produce the same result from two languages is not something that is usually performed or is desirable. Sticking to approximations is good. However, numpy and other image processing libraries try in many cases to emulate the Matlab results so I was hoping to find a way to reproduce the exact result. This might be a thing, afterall. I am still going to try and keep the question open for a bit to see if anyone has a way to solve it, otherwise, I will choose yours as the answer since it is a very good answer to the question. – Jorge Diaz Oct 20 '19 at 12:33
1

I was able to get a relatively good match (better than 10e-10 percent) between imgaussfilt called with the circular padding option and scipy.ndimage.gaussian_filter called with mode='wrap' parameter, provided that imgaussfilt uses the default filter size of 2*ceil(2*sigma)+1 and gaussian_filter specifies truncate=np.ceil(2*sigma)/sigma. These settings result in the same radius and overall size for the underlying Gaussian kernel, for the circular padding case. I have not tested other padding types.

Ekin
  • 111
  • 2
1

Recently I was working on the conversion of a Matlab code to python and looking about the same on the internet and docs for the functions found the exact same results with OpenCV implementation.

% In Matlab
imgaussfilt(image, 100)

is equivalent to

import cv2
# in python
cv2.GaussianBlur(image, ksize=(0, 0), sigmaX=100, borderType=cv2.BORDER_REPLICATE)

enter image description here

avisionx
  • 220
  • 1
  • 8
  • Exact as in identical floating-point values? That is highly unlikely. Even different versions of MATLAB might not produce numerically identical output. Also, OpenCV in that code likely produces integer output, not floating-point output as OP’s case. – Cris Luengo Nov 15 '20 at 14:50
  • I think this attempt is worthy of praise as it is an approximation to a potentially exact solution. I have upvoted -and have highlighted it onmy original question- it because I think that people looking to solve this problem in the future might take into consideration this approach. However, as @Cris Luengo says the correct answer still stands: there is no perfect replication of the Matlab result into Python. – Jorge Diaz Nov 15 '20 at 18:42
  • @CrisLuengo By exact here I could get the similar/closest result in Matlab and python obviously how rounding and algorithms are written in each library will affect the output but my main point is borderType=cv2.BORDER_REPLICATE is something one needs to change as by default output won't be same – avisionx Nov 17 '20 at 15:39
  • @avisionx: OP already had results identical up to 6 decimal digits, and an answer showed equality up to 10 decimal digits. They obviously figured out how to make the boundary condition match up. This answer just shows visually similar results, I don't see how it adds any value. – Cris Luengo Nov 17 '20 at 15:57