Computing a neighborhood window efficiently

[Serafino]

I'm trying to implement a version of the Spatial Fuzzy C-Means as a plugin for image segmentation. I need to calculate a spatial function over a neighborhood window for each pixel of an image and this step slows thing down a bit if I simply implement the windowing as in the mean shift filter (ie by simply cycling through each each and creating a window ex nihilo).
I thought up some optimizations, but doing them in Java leads to worse results than the unoptimized one!
Here is a link to the question I posted on StackOverflow (it shows the code for the algorithm optimization).
http://stackoverflow.com/questions/4699396/java-code-optimization-on-matrix-windowing-computes-in-more-time
Is there another way to compute that? What am I doing wrong in Java ?(in C++ it is optimized)

[Gabriel Landini]

On Sunday 16 Jan 2011, Serafino <[hidden email]> wrote:
> I'm trying to implement a version of the Spatial Fuzzy C-Means as a plugin
> for image segmentation. I need to calculate a spatial function over a
> neighborhood window for each pixel of an image and this step slows thing
> down a

In stackoverfow you posted:

> I have a matrix which represents an image and I need to cycle over each
> pixel  and for each one of those I have to compute the sum of all its
> neighbors, ie  the pixels that belong to a window of radius rad centered on
> the pixel.

If it is the "sum" you want, I would take a look at the implementation of the
built-in mean filter. Just do not divide by the number of pixels and you've
got the sum computed.

Cheers

Gabriel

[dscho]

Hi,

On Sun, 16 Jan 2011, Serafino wrote:

> I'm trying to implement a version of the Spatial Fuzzy C-Means as a
> plugin for image segmentation. I need to calculate a spatial function
> over a neighborhood window for each pixel of an image and this step
> slows thing down a bit if I simply implement the windowing as in the
> mean shift filter (ie by simply cycling through each each and creating a
> window ex nihilo). I thought up some optimizations, but doing them in
> Java leads to worse results than the unoptimized one!

In Java, creating a lot of objects or indirection leads to suboptimal
results.

For example, using a 1-dimensional array instead of a 2-dimensional one
helps.

For example, using a fixed-size array always beats using a collection.

For example, using int instead of Integer accelerates the code. Note that
this is not possible with Java Collections; your int would always be boxed
into Integer instances.

In general, there are quite a few pitfalls for C++ programmers who think
they can program Java just the same way. Have a look at
http://pacific.mpi-cbg.de/wiki/index.php/Tips_for_C%2B%2B_developers for
a list of common mistakes (not including your problems, though).

You specific question, how sums of rectangular[*1*] neighborhoods can be
calculated efficiently, has a totally different answer, though: construct
a second matrix S containing the sum of the rectangles between the upper
left corner and the current coordinate. Then, the sum of arbitrary
rectangles (x1 - x2, y1 - y2) can be calculated as

         S(x1, y1) + S(x2, y2) - S(x1, y2) - S(x2, y1)

A sketch using ImageJ's API for 32-bit images:

        int w = image.getWidth(), h = image.getHeight();
        float[] pixels = (float[])image.getProcessor().getPixels();
        float[] sum = new int[w * h];

        sum[0] = pixels[0];
        for (int i = 1; i < w; i++)
                sum[i] = sum[i - 1] + pixels[i];
        for (int j = 1; j < h; j++) {
                sum[j * w] = sum[(j - 1) * w] + pixels[j * w];
                for (int i = 1; i < w; i++)
                        sum[j * w + i] = sum[(j - 1) * w + i]
                                + pixels[j * w + i - 1] + pixels[j * w + i];
        }

Note: I wrote this in the mail client, so it is not even compile-tested.

Note2: there might be another 1% performance gain from using
        for (int j = w; j < w * h; j += w)
together with the obvious changes in the rest of the code.

Note3: the performance of this algorithm does not depend on the size of
the neighborhood, but the memory complexity is larger than your
solutions'.

Ciao,
Johannes

Footnote *1*: you almost mislead me by using the term "radius" suggesting
circular neighborhoods :-)

[Stephan Saalfeld]

Hi,

getting the sum or mean of arbitrarily large rectangles can be
accelerated by preprocessing the image into an integral image:

http://en.wikipedia.org/wiki/Summed_area_table

That way, accessing the sum is constant in time, namely 4 add
operations:

I(x_max,y_max) + I(x_min,y_min) - I(x_min,y_max) - I(x_max,y_min)

Note that you have to make sure that this will get you quite large
numbers which requires storing them in a larger type (e.g. summing up
bytes into byte wouldn't work, you would have to go for short, int or
long).  No insertion or replacement required, that is you can use []
instead of Collections which, in Java, is almost always faster and
allows you to not use Objects but basic types that do not acquire space
in the heap.

Do that first and later check how to optimize the code.

Good luck,
Stephan


On Sun, 2011-01-16 at 03:36 -0800, Serafino wrote:

> I'm trying to implement a version of the Spatial Fuzzy C-Means as a plugin
> for image segmentation. I need to calculate a spatial function over a
> neighborhood window for each pixel of an image and this step slows thing
> down a bit if I simply implement the windowing as in the mean shift filter
> (ie by simply cycling through each each and creating a window ex nihilo).
> I thought up some optimizations, but doing them in Java leads to worse
> results than the unoptimized one!
> Here is a link to the question I posted on StackOverflow (it shows the code
> for the algorithm optimization).
> http://stackoverflow.com/questions/4699396/java-code-optimization-on-matrix-windowing-computes-in-more-time
> http://stackoverflow.com/questions/4699396/java-code-optimization-on-matrix-windowing-computes-in-more-time 
> Is there another way to compute that? What am I doing wrong in Java ?(in C++
> it is optimized)
>

[Stephan Saalfeld]

Hi Johannes and all,

ah---sorry for repeating your reply---haven't read it before posting...

Note that the given example suffers from writing the sums of many many
floats into float.  That might be o.k. for images of average size but
could generate problems due to float having limited accuracy for larger
numbers than for smaller ones.  This may not only affect the accuracy of
the sum itself but can accumulate errors during summing up.  Adding a
very small B to a very large A may result in A.  One would have to
carefully check if, in a specific situation, the results are fine.

To prevent error accumulation during summing up, I once wrote a helper
class:

http://pacific.mpi-cbg.de/cgi-bin/gitweb.cgi?p=mpicbg.git;a=blob;f=mpicbg/util/RealSum.java

which stores intermediate sums and always adds up sums of equally sized
sets of numbers.  By that, it reduces the effect of the accumulated
mistakes.

Best,
Stephan



On Sun, 2011-01-16 at 13:55 +0100, Johannes Schindelin wrote:

> Hi,
>
> On Sun, 16 Jan 2011, Serafino wrote:
>
> > I'm trying to implement a version of the Spatial Fuzzy C-Means as a
> > plugin for image segmentation. I need to calculate a spatial function
> > over a neighborhood window for each pixel of an image and this step
> > slows thing down a bit if I simply implement the windowing as in the
> > mean shift filter (ie by simply cycling through each each and creating a
> > window ex nihilo). I thought up some optimizations, but doing them in
> > Java leads to worse results than the unoptimized one!
>
> In Java, creating a lot of objects or indirection leads to suboptimal
> results.
>
> For example, using a 1-dimensional array instead of a 2-dimensional one
> helps.
>
> For example, using a fixed-size array always beats using a collection.
>
> For example, using int instead of Integer accelerates the code. Note that
> this is not possible with Java Collections; your int would always be boxed
> into Integer instances.
>
> In general, there are quite a few pitfalls for C++ programmers who think
> they can program Java just the same way. Have a look at
> http://pacific.mpi-cbg.de/wiki/index.php/Tips_for_C%2B%2B_developers for
> a list of common mistakes (not including your problems, though).
>
> You specific question, how sums of rectangular[*1*] neighborhoods can be
> calculated efficiently, has a totally different answer, though: construct
> a second matrix S containing the sum of the rectangles between the upper
> left corner and the current coordinate. Then, the sum of arbitrary
> rectangles (x1 - x2, y1 - y2) can be calculated as
>
> S(x1, y1) + S(x2, y2) - S(x1, y2) - S(x2, y1)
>
> A sketch using ImageJ's API for 32-bit images:
>
> int w = image.getWidth(), h = image.getHeight();
> float[] pixels = (float[])image.getProcessor().getPixels();
> float[] sum = new int[w * h];
>
> sum[0] = pixels[0];
> for (int i = 1; i < w; i++)
> sum[i] = sum[i - 1] + pixels[i];
> for (int j = 1; j < h; j++) {
> sum[j * w] = sum[(j - 1) * w] + pixels[j * w];
> for (int i = 1; i < w; i++)
> sum[j * w + i] = sum[(j - 1) * w + i]
> + pixels[j * w + i - 1] + pixels[j * w + i];
> }
>
> Note: I wrote this in the mail client, so it is not even compile-tested.
>
> Note2: there might be another 1% performance gain from using
> for (int j = w; j < w * h; j += w)
> together with the obvious changes in the rest of the code.
>
> Note3: the performance of this algorithm does not depend on the size of
> the neighborhood, but the memory complexity is larger than your
> solutions'.
>
> Ciao,
> Johannes
>
> Footnote *1*: you almost mislead me by using the term "radius" suggesting
> circular neighborhoods :-)

[Serafino]

Thank you for pointing this out, I'll have a look at it

[Serafino]

Hi Johannes,
Thank you for replying and posting the C++-to-java link (by the way I'm aware about what is written there). I guess the main disadvantage with collection, like someone pointed out on stackoverflow.com is due to the conversions between primitive types and wrappers.
I'll have now a look at the integral image technique (thanks to Stephan Saalfeld also ^^), and I'll let you know what I can accomplish.
The code I wrote till now for the plugin is located at https://github.com/arranger1044/SFCM if anyone is interested.

PS
we say "window radius" (well, the italian equivalent) and maybe we are misusing the word

[Serafino]

It seems that the mean filter does something like my second optimization, for each row it calculates the sum and for each row it subtract the values gone out the window and sums the new ones that entered the window (it uses arrays and not collections). I did a new version of it using one array for the row and now it performs better but still worse than the simple way sometimes.
I will do more work on it to have a deeper understanding.
Concerning the integral image I guess that the formula to compute the sum vector shall be:
        int[] sum = new int[w * h];

        sum[0] = mat[0];
        for (int i = 1; i < w; i++)
                sum[i] = sum[i - 1] + mat[i];
        for (int j = 1; j < h; j++)
        {
            sum[j * w] = sum[(j - 1) * w] + mat[j * w];
            for (int i = 1; i < w; i++)
                sum[j * w + i] = sum[(j - 1) * w + i] + sum[j * w + i - 1] +
                       mat[j * w + i] - sum[(j - 1) * w + i -1];
        }
am I wrong?

[Stephan Saalfeld]

Hi,

> [...]
> Concerning the integral image I guess that the formula to compute the sum
> vector shall be:
> int[] sum = new int[w * h];
>
> sum[0] = mat[0];
> for (int i = 1; i < w; i++)
> sum[i] = sum[i - 1] + mat[i];
> for (int j = 1; j < h; j++)
>         {
>             sum[j * w] = sum[(j - 1) * w] + mat[j * w];
>             for (int i = 1; i < w; i++)
> sum[j * w + i] = sum[(j - 1) * w + i] + sum[j * w + i - 1] +
>                        mat[j * w + i] - sum[(j - 1) * w + i -1];
> }
> am I wrong?

No---looks correct.  I would just save a few operations and array
accesses and, probably, put it into a long:

long[] sum = new long[w * h];
int s = 0;
for (int i = 0; i < w; ++i) {
        s += mat[i]
        sum[i] = s;
}
for (int j = 1; j < h; ++j)
{
        int jw = j * w;
        sum[jw] = sum[jw - w] + mat[jw];
        for (int i = 1; i < w; ++i) {
                int jwi = jw + i;
                sum[jwi] =
                        sum[jwi - w] + sum[jwi - 1] +
                        mat[jwi] - sum[jwi - w - 1];
}

Best,
Stephan

[Stephan Saalfeld]

Hi,

infected by the discussion on integral images I wrote a small plugin
that demonstrates the speed gain with an interactive mean-smooth.  With
an image open, call the plugin `Integral Smooth' and click and drag with
the mouse in the image canvas to change the size of the smooth-kernel.
ENTER leaves the smoothed image, ESC cancels.

http://fly.mpi-cbg.de/saalfeld/download/integral_smooth.jar

The jar contains both binaries and sources.

Best,
Stephan




On Mon, 2011-01-17 at 00:55 +0100, Stephan Saalfeld wrote:

> Hi,
>
> > [...]
> > Concerning the integral image I guess that the formula to compute the sum
> > vector shall be:
> > int[] sum = new int[w * h];
> >
> > sum[0] = mat[0];
> > for (int i = 1; i < w; i++)
> > sum[i] = sum[i - 1] + mat[i];
> > for (int j = 1; j < h; j++)
> >         {
> >             sum[j * w] = sum[(j - 1) * w] + mat[j * w];
> >             for (int i = 1; i < w; i++)
> > sum[j * w + i] = sum[(j - 1) * w + i] + sum[j * w + i - 1] +
> >                        mat[j * w + i] - sum[(j - 1) * w + i -1];
> > }
> > am I wrong?
>
> No---looks correct.  I would just save a few operations and array
> accesses and, probably, put it into a long:
>
> long[] sum = new long[w * h];
> int s = 0;
> for (int i = 0; i < w; ++i) {
> s += mat[i]
> sum[i] = s;
> }
> for (int j = 1; j < h; ++j)
> {
> int jw = j * w;
> sum[jw] = sum[jw - w] + mat[jw];
> for (int i = 1; i < w; ++i) {
> int jwi = jw + i;
> sum[jwi] =
> sum[jwi - w] + sum[jwi - 1] +
> mat[jwi] - sum[jwi - w - 1];
> }
>
> Best,
> Stephan

[Serafino]

That plugin is supernice ^^ thank you for sharing it.
I had to tweak your formula for my purposes since it did not computed the precise square window sum since I have odd square sides.
Once I implemented it , it results the fastest approach when the window 'radius' is growing, but at small level the column cached array approach is better (it is always better for memory saving purposes).
I guess that I'm not that good at comparing algorithm in Java (I still need to find a good profiler), in the meantime I will attach the java code showing all the methods I wrote as a proof of concept
Main.java