Questions aboutGaussianBlur

[José]

Hi,

I'm studying the code in
https://github.com/imagej/imagej1/blob/master/ij/plugin/filter/GaussianBlur.javaand
I have some questions about some code and the maths behind it:

1) On lines 229-235 there are some adjustments to sigma based on
downscaling the image, applying the filter and upscaling again that i don't
fully understand. Could you explain more clearly that adjustments or
suggest any references i can check about applying the Gaussian blur with
that downscale-convolve-upscale approach?

2) On lines 179-186 one can see that, when bluring, extra lines are added
to the X component but not to the Y component. Why is that?

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Michael Schmid]

Hi José,

unfortunately I am not aware of any reference for the Gaussian Blur
downscaling algorithm.  When I added this to ImageJ (because the Gaussien
Blur with large sigma was too slow for me) I did search the literature for
ways to improve it and tried various methods described in the literature,
but none was satisfying.
So I did the math for downscaling and upscaling myself; it is pretty simple:
- Downscaling already slightly smooths the data.
- Then the line (or column) of the downscaled image is blurred the
conventional way.
- Finally, a kind of 'soft' cubic interpolation algorithm is used to
upscale the line (or column) to the original size.
  I just had a look at the kernels used for downscaling and upscaling and
calculate the standard deviation sigma for each of them.  Since we have
successive convolutions, one can simply add the squares of the sigma
values of them.
  For the upscaling and downscaling kernels, see the makeDownscaleKernel
and makeUpscaleKernel functions.

If I find time I might write a short paper on it, but this won't be in the
immediate future.

---

Concerning point (2), extra lines when blurring in x direction, this can
be explained more easily:
Convolving with a 2D Gaussian is separable, one can first convolve with a
1D Gaussian in x and then with a 1D Gaussian in y.
This is straightforward when convolving the full image. If one wants to
treat only a part of the image, one must make sure that everything that is
an input for blurring in y direction has been blurred in x direction
before.  E.g., if we need the output for lines 50-100, and blurring in y
has an overall kernel size of +/- 10 lines, blurring in y needs lines
40-110 as an input.  Therefore we first need to blur lines 40-110 as an
input for blurring in y direction.

I hope that this is clear...

Michael
_________________________________________________________________



On Thu, February 20, 2014 13:40, José  wrote:

> Hi,
>
> I'm studying the code in
> https://github.com/imagej/imagej1/blob/master/ij/plugin/filter/GaussianBlur.javaand
> I have some questions about some code and the maths behind it:
>
> 1) On lines 229-235 there are some adjustments to sigma based on
> downscaling the image, applying the filter and upscaling again that i
> don't
> fully understand. Could you explain more clearly that adjustments or
> suggest any references i can check about applying the Gaussian blur with
> that downscale-convolve-upscale approach?
>
> 2) On lines 179-186 one can see that, when bluring, extra lines are added
> to the X component but not to the Y component. Why is that?
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Michael Schmid]

Hi José,

the 'edge correction' for the Gaussian kernel is shortly described in the description of the method (line 504):

"To avoid a step due to the cutoff at a finite value, the near-edge values are replaced by a 2nd-order polynomial with its minimum=0 at the first out-of-kernel  pixel. Thus, the kernel function has a smooth 1st derivative in spite of finite length."

In other words:
In principle, a Gaussian kernel would have infinite size.  Although the Gaussian function decreases rather quickly, the full kernel size required for 32-bit floats with the full precision would be rather large and slow down the calculation.
Simply truncating the kernel would lead to a step, which is undesirable.
E.g., when doing a derivative or Laplace filter after the Gaussian Blur, a kernel with a sharp cutoff will create artifacts.

So there is a 'transition zone' between the exact Gaussian and the edge of the kernel, where the function decreases with a 2nd-order polynomial. Thus, the first derivative of the kernel is smooth everywhere.

You can see how well it works with the following macro:

newImage("Untitled", "32-bit black", 100, 100, 1);
makeRectangle(50, 50, 1, 1);
run("Set...", "value=100");
run("Select None");
run("Gaussian Blur...", "sigma=7");
run("Convolve...", "text1=[0 -1 0\n-1 4 -1\n0 -1 0\n]");

The last command is a Laplacian filter, i.e., the 2nd derivative, so it is highly sensitive to any deviation from an ideal Gaussian.  You need to crank up the contrast a lot (-0.0001 to +0.0001) to see four dark lines around the center annulus.  These lines are artifacts due to the finite size of the kernel.

For comparison, if you filter with a constant ('rectangular') kernel three times, which is the usual algorithm often used as a 'low-quality Gaussian blur', the Laplacian will still look awful.  Nevertheless, the three rectangular blurs will be hardly any faster.  Unless one carefully selects the kernel sizes, one needs five iterations of convolution with the constant kernel for better performance that the current implementation of the Gaussian blur.  In addition, in the usual simple implementation, multiple blurring with rectangular kernels gives you only some selected values of the 'blur length' sigma.

By the way, the 'edge correction' needs some minimum kernel size to work.  That's why the 'accuracy' parameter, which defines the kernel size, should not be above 0.02.

So far a bit about the background of the GaussianBlur code...
---

By the way, I start to wonder why you are analyzing the GaussianBlur in such detail.   Is there a problem with it?


Michael
________________________________________________________________


On Feb 24, 2014, at 23:15, José wrote:

> Thanks. Also, could you give me some information about how the 'edge correction' works (what does it do and how?)? (line 532)
>
>
> 2014-02-22 17:38 GMT+01:00 Michael Schmid <[hidden email]>:
> Hi José,
>
> unfortunately I am not aware of any reference for the Gaussian Blur
> downscaling algorithm.  When I added this to ImageJ (because the Gaussien
> Blur with large sigma was too slow for me) I did search the literature for
> ways to improve it and tried various methods described in the literature,
> but none was satisfying.
> So I did the math for downscaling and upscaling myself; it is pretty simple:
> - Downscaling already slightly smooths the data.
> - Then the line (or column) of the downscaled image is blurred the
> conventional way.
> - Finally, a kind of 'soft' cubic interpolation algorithm is used to
> upscale the line (or column) to the original size.
>   I just had a look at the kernels used for downscaling and upscaling and
> calculate the standard deviation sigma for each of them.  Since we have
> successive convolutions, one can simply add the squares of the sigma
> values of them.
>   For the upscaling and downscaling kernels, see the makeDownscaleKernel
> and makeUpscaleKernel functions.
>
> If I find time I might write a short paper on it, but this won't be in the
> immediate future.
>
> ---
>
> Concerning point (2), extra lines when blurring in x direction, this can
> be explained more easily:
> Convolving with a 2D Gaussian is separable, one can first convolve with a
> 1D Gaussian in x and then with a 1D Gaussian in y.
> This is straightforward when convolving the full image. If one wants to
> treat only a part of the image, one must make sure that everything that is
> an input for blurring in y direction has been blurred in x direction
> before.  E.g., if we need the output for lines 50-100, and blurring in y
> has an overall kernel size of +/- 10 lines, blurring in y needs lines
> 40-110 as an input.  Therefore we first need to blur lines 40-110 as an
> input for blurring in y direction.
>
> I hope that this is clear...
>
> Michael
> _________________________________________________________________
>
>
>
> On Thu, February 20, 2014 13:40, José  wrote:
> > Hi,
> >
> > I'm studying the code in
> > https://github.com/imagej/imagej1/blob/master/ij/plugin/filter/GaussianBlur.javaand
> > I have some questions about some code and the maths behind it:
> >
> > 1) On lines 229-235 there are some adjustments to sigma based on
> > downscaling the image, applying the filter and upscaling again that i
> > don't
> > fully understand. Could you explain more clearly that adjustments or
> > suggest any references i can check about applying the Gaussian blur with
> > that downscale-convolve-upscale approach?
> >
> > 2) On lines 179-186 one can see that, when bluring, extra lines are added
> > to the X component but not to the Y component. Why is that?
> >
> > --
> > ImageJ mailing list: http://imagej.nih.gov/ij/list.html
> >
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Michael Schmid]

Hi Jose,


[I'll copy this to the mailing list; in case someone else is interested in the inner workings of the GaussianBlur...]


The idea for the 'edge correction' code, doing a smooth transition from the Gaussian to a polynomial that becomes zero at the kernel boundary, is as follows:

Assume a second-order polynomial (parabola) which is zero at the first point outside the 'support' of the kernel:
  y(x) = a^2*(x-kRadius)^2

This polynomial ensures that there is no step where the kernel becomes zero and also the first derivative is smooth.

At some point r, its value should be equal to that of the Gaussian, which is kernel[0][r]:
  y(r) = kernel[0][r]               yields
  a^2 = kernel[0][r]/(x-kRadius)^2
  a = sqrt(kernel[0][r]) / (x-kRadius)

The polynomial should touch the Gaussian from the bottom at the transition point, then also the 1st derivative is smooth.  (The Gaussian is nonzero at x=kRadius, so it must be above the polynomial there)
So we search for the polynomial with the smallest value 'a', where the equation is fulfilled in some point.  To understand this, it is best to draw the two functions, the tail of the Gaussian and the parabola.
When stepping in with r from the outermost point, the calculated 'a' value will become smaller and smaller, until a minimum is found.  That's where the Gaussian and the polynomial touch in just one point, and we keep that value 'a' as variable 'sqrtSlope' (the name is not very nice, is should better be aMin).
If the 'a' value gets larger again, we have stepped too far in r.  At this r value + 1, the polynomial and Gaussian have the same value.  We have to calculate the values of the polynomial starting at this r value + 2, and this is finally done in line 544.

By the way, this algorithm does not work if the support size of the Gaussian is too small, e.g. 1 sigma each direction.  That's why it is specified that 'accuracy' should not be > 0.02.  If I remember it correctly (I wrote that code in 2007), it would even work up to somewhat higher values, but I wanted to keep it on the safe side.
The accuracy values actually used by ImageJ are 0.002 for 8-bit data and 0.0002 for 16-bit and float data.


Hope this helps,

Michael
________________________________________________________________
On Mar 14, 2014, at 17:58, José wrote:

> Thanks for the explanations!
>
> I am having some troubles trying to understand how the polynomial is found (lines 532-542 on https://github.com/imagej/imagej1/blob/master/ij/plugin/filter/GaussianBlur.java )
>
> Could you clarify where does that expresion (Math.sqrt(kernel[0][r])/(kRadius-r)) come from and how the two degree polynomial is calculated?
>
> José
>
>
> 2014-02-25 17:12 GMT+01:00 Michael Schmid <[hidden email]>:
> Hi José,
>
> the 'edge correction' for the Gaussian kernel is shortly described in the description of the method (line 504):
>
> "To avoid a step due to the cutoff at a finite value, the near-edge values are replaced by a 2nd-order polynomial with its minimum=0 at the first out-of-kernel  pixel. Thus, the kernel function has a smooth 1st derivative in spite of finite length."
>
> In other words:
> In principle, a Gaussian kernel would have infinite size.  Although the Gaussian function decreases rather quickly, the full kernel size required for 32-bit floats with the full precision would be rather large and slow down the calculation.
> Simply truncating the kernel would lead to a step, which is undesirable.
> E.g., when doing a derivative or Laplace filter after the Gaussian Blur, a kernel with a sharp cutoff will create artifacts.
>
> So there is a 'transition zone' between the exact Gaussian and the edge of the kernel, where the function decreases with a 2nd-order polynomial. Thus, the first derivative of the kernel is smooth everywhere.
>
> You can see how well it works with the following macro:
>
> newImage("Untitled", "32-bit black", 100, 100, 1);
> makeRectangle(50, 50, 1, 1);
> run("Set...", "value=100");
> run("Select None");
> run("Gaussian Blur...", "sigma=7");
> run("Convolve...", "text1=[0 -1 0\n-1 4 -1\n0 -1 0\n]");
>
> The last command is a Laplacian filter, i.e., the 2nd derivative, so it is highly sensitive to any deviation from an ideal Gaussian.  You need to crank up the contrast a lot (-0.0001 to +0.0001) to see four dark lines around the center annulus.  These lines are artifacts due to the finite size of the kernel.
>
> For comparison, if you filter with a constant ('rectangular') kernel three times, which is the usual algorithm often used as a 'low-quality Gaussian blur', the Laplacian will still look awful.  Nevertheless, the three rectangular blurs will be hardly any faster.  Unless one carefully selects the kernel sizes, one needs five iterations of convolution with the constant kernel for better performance that the current implementation of the Gaussian blur.  In addition, in the usual simple implementation, multiple blurring with rectangular kernels gives you only some selected values of the 'blur length' sigma.
>
> By the way, the 'edge correction' needs some minimum kernel size to work.  That's why the 'accuracy' parameter, which defines the kernel size, should not be above 0.02.
>
> So far a bit about the background of the GaussianBlur code...
> ---
>
> By the way, I start to wonder why you are analyzing the GaussianBlur in such detail.   Is there a problem with it?
>
>
> Michael
> ________________________________________________________________
>
>
> On Feb 24, 2014, at 23:15, José wrote:
>
> > Thanks. Also, could you give me some information about how the 'edge correction' works (what does it do and how?)? (line 532)
> >
> >
> > 2014-02-22 17:38 GMT+01:00 Michael Schmid <[hidden email]>:
> > Hi José,
> >
> > unfortunately I am not aware of any reference for the Gaussian Blur
> > downscaling algorithm.  When I added this to ImageJ (because the Gaussien
> > Blur with large sigma was too slow for me) I did search the literature for
> > ways to improve it and tried various methods described in the literature,
> > but none was satisfying.
> > So I did the math for downscaling and upscaling myself; it is pretty simple:
> > - Downscaling already slightly smooths the data.
> > - Then the line (or column) of the downscaled image is blurred the
> > conventional way.
> > - Finally, a kind of 'soft' cubic interpolation algorithm is used to
> > upscale the line (or column) to the original size.
> >   I just had a look at the kernels used for downscaling and upscaling and
> > calculate the standard deviation sigma for each of them.  Since we have
> > successive convolutions, one can simply add the squares of the sigma
> > values of them.
> >   For the upscaling and downscaling kernels, see the makeDownscaleKernel
> > and makeUpscaleKernel functions.
> >
> > If I find time I might write a short paper on it, but this won't be in the
> > immediate future.
> >
> > ---
> >
> > Concerning point (2), extra lines when blurring in x direction, this can
> > be explained more easily:
> > Convolving with a 2D Gaussian is separable, one can first convolve with a
> > 1D Gaussian in x and then with a 1D Gaussian in y.
> > This is straightforward when convolving the full image. If one wants to
> > treat only a part of the image, one must make sure that everything that is
> > an input for blurring in y direction has been blurred in x direction
> > before.  E.g., if we need the output for lines 50-100, and blurring in y
> > has an overall kernel size of +/- 10 lines, blurring in y needs lines
> > 40-110 as an input.  Therefore we first need to blur lines 40-110 as an
> > input for blurring in y direction.
> >
> > I hope that this is clear...
> >
> > Michael
> > _________________________________________________________________
> >
> >
> >
> > On Thu, February 20, 2014 13:40, José  wrote:
> > > Hi,
> > >
> > > I'm studying the code in
> > > https://github.com/imagej/imagej1/blob/master/ij/plugin/filter/GaussianBlur.javaand
> > > I have some questions about some code and the maths behind it:
> > >
> > > 1) On lines 229-235 there are some adjustments to sigma based on
> > > downscaling the image, applying the filter and upscaling again that i
> > > don't
> > > fully understand. Could you explain more clearly that adjustments or
> > > suggest any references i can check about applying the Gaussian blur with
> > > that downscale-convolve-upscale approach?
> > >
> > > 2) On lines 179-186 one can see that, when bluring, extra lines are added
> > > to the X component but not to the Y component. Why is that?
> > >
> > > --
> > > ImageJ mailing list: http://imagej.nih.gov/ij/list.html
> > >
> >
> > --
> > ImageJ mailing list: http://imagej.nih.gov/ij/list.html
> >
>
>

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html