Image smoothing with own smoothing function

I have thousands images similar to the following one: This image was produced with the following color table:

colTable = {{Black},
Table[{Blend[{Blue, Green, Yellow, Red}, x]}, {x, 1/255, 1,
1/255}]};
colTable = Flatten[colTable];
...
image =
Colorize[Image[image], ColorFunction -> (Blend[colTable, #] &)]

Now I would like to smooth this image pixel by pixel following a defined function.

The idea behind the smoothing is: I want to calculate around each pixel the number of black pixels plus the number of red pixels in a defined sourrounding rectangle (e.g. rectangle size: dx=dy=40 pixels). This number value is devided by the number of rectangle pixels and the final value determines a certain smoothed color for the corresponding pixel.

The code which I wrote is the following:

dim = ImageDimensions[image];

sx = dim[];
sy = dim[];

dx = 40;
dy = 40;

smoothImage = ConstantArray[0, {sx, sy}];

Table[

Table[

subImage =
ImageTake[
image, {iy - dy/2, iy + dy/2 - 1}, {ix - dx/2, ix + dx/2 - 1}];
data = ImageData[subImage];
sumRB =
Count[data, {0., 0., 0.}, Infinity] +
Count[data, {1., 0., 0.}, Infinity];
smoothValues[[ix, iy]] = sumRB/(dx*dy);

, {ix, dx/2 + 1, sx - dx/2 + 1}

];

, {iy, dy/2 + 1, sy - dy/2 + 1}

]

max = Max[smoothValues];

smoothImage = Image[smoothValue/max];
smoothImage =
Colorize[Image[smoothImage], ColorFunction -> (Blend[colTable, #] &)]

The problem is: this code is extremely slow and lasts for hours (only for one image) ... it might also be I made somewhere a mistake ...

Is there another way (e.g. by compilation, parallel processing, cuda or any other fast routines) to solve this problem?

My expectation is to get such a smoothed image depending on the color table and the used color range (this image I have produced with IDL which takes 1 sec). • Why don't you process the source image (before the Colorize step) - wouldn't that make the counting easier? – Niki Estner Jul 7 '15 at 17:29
• could you not use ImageFilter ? also ParallelTable would be a quick way to parallelize. – amr Jul 7 '15 at 17:37

As long as you're only counting/summing pixel values in a sliding window, you don't have to recalculate the whole sum for every pixel. You can just "count" red/black pixels in a "1x1 window" and then use a moving average or "box" filter to sum over neighborhoods.

First calculate two images where every red/black pixel is 1, other pixels 0:

image = Import["http://i.stack.imgur.com/Pgppv.png"]
dx = 40;
dy = 40;

red = {1, 0, 0};
black = {0, 0, 0};
redPixels =
Binarize[ColorNegate@ImageApply[Norm[# - red] &, image], 10^-5];
blackPixels =
Binarize[ColorNegate@ImageApply[Norm[# - black] &, image], 10^-5];

Then apply a box filter to those images and add them:

box = BoxMatrix[Floor[{dx, dy}/2]];
ImageConvolve[blackPixels, box]] // ImageAdjust • Norm[# - {0, 0, 0}] & is the same as plain Norm, no? – J. M. will be back soon Jul 7 '15 at 18:38
• @Guesswhoitis. I wanted to make the step from Count[data, {0., 0., 0.}, Infinity] to Norm[# - {0,0,0}] more obvious. – Niki Estner Jul 8 '15 at 6:18
• Ah, of course; more illustrative that way. :) – J. M. will be back soon Jul 8 '15 at 6:26
• @nikie: great and fast. thanks a lot. How would the code look like for a circular filter with a certain radius r? – mrz Jul 8 '15 at 9:45
• @mrz: You would simply use DiskMatrix instead of BoxMatrix – Niki Estner Jul 8 '15 at 9:59

You can use ColorReplace:

image = Import["http://i.stack.imgur.com/Pgppv.png"];
dx = 40;
dy = 40;
redblack =
ColorReplace[image, {Red -> White, Black -> White, _ -> Black}];
box = BoxMatrix[Floor[{dx, dy}/2]]; 