SmoothingRecursiveGaussian vs DiscreteGaussian

Beginner Questions
,
1

Here is a small example using both SmoothingRecursiveGaussian and DiscreteGaussian to smooth an image. Both filter should use the same σ of 0.8 - as far as I understood the documentation, the Recursive filter has the standard deviation as an argument, while the Discrete filter wants the variance, thus 0.8^2:

import numpy as np
import SimpleITK as sitk
print(sitk.Version_VersionString())

arr = np.zeros((100, 100), dtype=np.float32)
arr[25:75, 25:] = 1400.0
arr[40:60, 40:60] = 0.0
img = sitk.GetImageFromArray(arr)

smoothing = sitk.SmoothingRecursiveGaussian(img, [0.8, 0.8])
discrete = sitk.DiscreteGaussian(img, [0.8**2, 0.8**2])

v_s = sitk.GetArrayViewFromImage(smoothing)
v_d = sitk.GetArrayViewFromImage(discrete)
print("original", arr.min(), arr.max())
print("smoothing", v_s.min(), v_s.max())
print("discrete", v_d.min(), v_d.max())

What I’m wondering is, that the output is the following:

2.3.1
original 0.0 1400.0
smoothing -2.6740482 1402.6741
discrete 0.0 1400.0

When using the Recursive Variant, I get values smaller and larger than the original image - which I find a bit weird, as I thought that the Gaussian should not produce values lower or larger than the original image. Is this maybe an implementation detail I missed?
Can I safely clamp the result of the filter to the original image range or is there another way to have the filter produce images with no voxels out of the original range?

The reason is, for large images the Recursive variant works much much faster for me, but I need to ensure that voxel values are not outside the original range.

1 Like
2

I would dive into some of the papers related to the details of the algorithm listed here:
https://itk.org/Doxygen/html/classitk_1_1RecursiveGaussianImageFilter.html

It may be interesting to inspect at which pixel locations the outliers and discrepancies occur. It may be related to boundary conditions.

P.S.

The following “flat static” usage has been deprecated:

In favor of:

sitk.Version.VersionString()
1 Like
3

That’s why I deliberately added voxels with non-zero value at the boundary. However, placing the square in the middle of the image does not change the result.

Here is a plot of the Recursive variant and marked all pixels that are <0 or >1400:
myplot

Interesting is the pattern of <0 pixels…

Thanks, I’ll start reading!

edit; Well okay, the first thing I notice is, that in the approximations of the Gaussian there are some constants. Could it just be that due to rounding you get these results?

4

Are you using float32 or a double? Have you compared?

5

I usually use float32. I tested it with the example above and it does not make a difference with double.

Another observation is, that the value depends on the value range. I.e., the maximal value will be value+1.337 and the minimal value-1.337 when a range of 700 is used, independent if I use a pixel value range of 0…700 or 700…1400.
When the range is twice as large, this error will be twice as large.

Also interesting, but might just be due to the image quality: If you look at Fig 2 and 3 here https://inria.hal.science/inria-00074778/document - it appears, that the value of the approximation at x=0 is larger than 1. I guess that would totally explain the behavior. However, I would have to look also in the code now what is actually implemented :smiley: