Gaussian Smoothing

[gusciora@…]

When applying the formula for Gaussian smoothing in the ADD, I do not obtain
meaningful results. Basically, the center coefficient c_0 will have some
non-zero value, while all other coefficients are zero or very close to zero.
The reason is that as n becomes large, then e is raised to a very large negative
number, which is of course, zero. The value I use for sigma in that formula is
basically the stdev for the entire vector (when calculating robust stats)
divided by 4. Am I applying the formula correctly?

[Paul Price]

sigma is a measure of the width of the smoothing. My first guess is that you're
taking a sigma in the coordinate, and applying it to the ordinate. Can you post
more information? The following might help:

• The values
• The measured sigma
• How you're applying the smoothing.

[gusciora@…]

I measure the clipped standard deviation for the entire image. I then divide
that by 4 and use that as the "sigma" in the ADD formula for the Gaussian
coefficients.
Some examples of the actual values are:

sigma is 0.511959

For a "Gaussian width" (corresponds to N in the ADD) of 20, the coefficients are:
p_psVectorsmoothHistGaussian(): gaussianCoefs[0] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[1] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[2] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[3] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[4] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[5] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[6] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[7] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[8] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[9] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[10] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[11] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[12] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[13] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[14] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[15] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[16] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[17] is 0.000002
p_psVectorsmoothHistGaussian(): gaussianCoefs[18] is 0.002616
p_psVectorsmoothHistGaussian(): gaussianCoefs[19] is 0.226155
p_psVectorsmoothHistGaussian(): gaussianCoefs[20] is 1.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[21] is 0.226155
p_psVectorsmoothHistGaussian(): gaussianCoefs[22] is 0.002616
p_psVectorsmoothHistGaussian(): gaussianCoefs[23] is 0.000002
p_psVectorsmoothHistGaussian(): gaussianCoefs[24] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[25] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[26] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[27] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[28] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[29] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[30] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[31] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[32] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[33] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[34] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[35] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[36] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[37] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[38] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[39] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[40] is 0.000000

How wide should the Gaussian width be?

[gusciora@…]

I measure the clipped standard deviation for the entire image. I then divide
that by 4 and use that as the "sigma" in the ADD formula for the Gaussian
coefficients.
Some examples of the actual values are:

sigma is 0.511959

For a "Gaussian width" (corresponds to N in the ADD) of 20, the coefficients are:
p_psVectorsmoothHistGaussian(): gaussianCoefs[0] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[1] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[2] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[3] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[4] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[5] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[6] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[7] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[8] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[9] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[10] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[11] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[12] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[13] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[14] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[15] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[16] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[17] is 0.000002
p_psVectorsmoothHistGaussian(): gaussianCoefs[18] is 0.002616
p_psVectorsmoothHistGaussian(): gaussianCoefs[19] is 0.226155
p_psVectorsmoothHistGaussian(): gaussianCoefs[20] is 1.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[21] is 0.226155
p_psVectorsmoothHistGaussian(): gaussianCoefs[22] is 0.002616
p_psVectorsmoothHistGaussian(): gaussianCoefs[23] is 0.000002
p_psVectorsmoothHistGaussian(): gaussianCoefs[24] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[25] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[26] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[27] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[28] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[29] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[30] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[31] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[32] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[33] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[34] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[35] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[36] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[37] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[38] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[39] is 0.000000
p_psVectorsmoothHistGaussian(): gaussianCoefs[40] is 0.000000

How wide should the Gaussian width be?

I did not understand your comment about taking a sigma in the coordinate and
applying it to the ordinate.

[Paul Price]

Looking at your function....

I think psHistogramAlloc and psVectorHistogram use different types (U32, S32).

Aside from the bad Gaussian (why aren't you using psGaussian?), it seems you're
confusing the definition of sigma --- it's the width of the distribution in real
units (whatever the units of the data are), not in number of bins.

Try this algorithm:

Iterate over each bin (i):

i_mid = (bounds[i+1] + bounds[i]) / 2.0
smooth[i] = 0
For 20 bins either side of bin i (j):

j_mid = (bounds[j+1] + bounds[j])/2.0
smooth += num[j]*psGaussian(j_mid, i_mid, sigma, true)

[gusciora@…]

I have corrected the S32/U32 bugs with the histogram.

Regarding the Gaussian code... I did not use psGaussian() because the ADD
instructed me to explicitly calculate the Gaussian coefficients and multiply
them by the histogram bin values. I don't mind using the algorithm you
suggested, however we'll have to put the new algorithm in the ADD. In fact, I
already coded (it's fairly straightforward), but the results do not seem correct.

I'm wondering if you forgot a statement in your algorithm where you divide
smooth by the total number of bins on either side of bin i...

I don't understand what you're saying about using the width of the distribution
in real units, not bins. I used sigma/4 as the sigma for the Gaussian
smoothing, as directed by the ADD.

[Paul Price]

The bin index is not the same as the value of the centre of the bin. In your
implementation, "sigma" is in the same units as the vector on which statistics
are being performed, but "i" is in units of bins. That is to say, the
dimensions in your Gaussian are not consistent. The Gaussian (neglecting the
normalisation) that is implemented in rel_3 is:

exp(-(float)((i - GAUSS_WIDTH) * (i - GAUSS_WIDTH))/(2.0 * sigma * sigma))

Say that the vector I put into psStats is the number of apples that different
people have. Then the dimensions of "sigma" is apples --- it's a measure of the
number of apples. But the dimensions of "i" is bin number, so that the exponent
is *not* dimensionless (as it should be), but is in units of "bin/apples". To
correct this, you need to calculate the mid-point of the bin, and use that
instead of the bin number.

Perhaps "algorithm" wasn't a good choice of word --- the "algorithm" I presented
above is merely an implementation of the algorithm specified in the ADD. Note
that the major difference between my implementation and the one in psLib is that
I am using consistent dimensions by using the value in the middle of the bin
(not the bin index) and sigma.

Using psGaussian will save you from making mistakes in calculating the Gaussian
coefficients --- presumably this code has already been tested, so you should
make use of it.

[gusciora@…]

Thanks for the clarification.

I agree that the algorithm you supplied in this bug report is correct, and it
produces meaningful results. However, I don't think that that algorithm is an
implementation of the algorithm in the ADD. The ADD doesn't specify that the
dimension of "i" should be in apples, or whatever; it specifies bin numbers, or
index numbers for the coefficients. I think I originally implemented what was
in the ADD faithfully.

[Paul Price]

Extract of ADD with the smoothing update

[Paul Price]

Closing subsequent to release of SDRS-08, ADD-07.