UBLAS Matrix Finding Surrounding Values of a Cell?

Question

I am looking for an elegant way to implement this. Basically i have a m x n matrix. Where each cell represents the pixel value, and the rows and columns represent the pixel rows and pixel columns of the image.

Since i basically mapped points from a HDF file, along with their corresponding pixel values. We basically have alot of empty pixels. Which are filled with 0.

Now what i need to do is take the average of the surrounding cell's, to average out of a pixel value for the missing cell.

Now i can brute force this but it becomes ugly fast. Is there any sort of elegant solution for this?

Answer 1

There's a well-known optimization to this filtering problem.

• Integrate the cells in one direction (say horizontally)
• Integrate the cells in the other direction (say vertically)
• Take the difference between each cell and it's N'th neighbor to the left.
• Take the difference between each cell and it's N'th lower neighbor

Like this:

    for (i = 0; i < h; ++i)
    for (j = 0; j < w-1; ++j)
       A[i][j+1] += A[i][j];
    for (i = 0; i < h-1; ++i)
    for (j = 0; j < w; ++j)
       A[i+1][j] += A[i][j]
    for (i = 0; i < h; ++i)
    for (j = 0; j < w-N; ++j)
       A[i][j] -= A[i][j+N];
    for (i = 0; i < h-N; ++i)
    for (j = 0; j < w; ++j)
       A[i][j] -= A[i-N][j];

What this does is:

• The first pass makes each cell the sum of all of the cells on that row to it's left, including itself.
• After the 2nd pass , each cell is the sum of all of the cells in a rectangle above and left of itselt (including it's own row and column)
• After the 3rd pass, each cell is the sum of a rectangle above and to the right of itself, N columns wide.
• After the 4th pass each cell is the sum of an NxN rectangle below and to the right of itself.

This takes 4 operations per cell to compute the sum, as opposed to 8 for brute force (assuming you're doing a 3x3 averaging filter).

The cool thing is that if you use ordinary two's-complement arithmetic, you don't have to worry about any overflows in the first two passes; they cancel out in the last two passes.

Answer 2

The main issues here are utilizing all available cores and cache effeciency.
You might be interested in checking fast implementation of convolution.
However, since you do it with Boost, you can check how this is done in this Boost example
I beleive you have to change only the convolution kernel for your specialized task.