Image smoothing with own smoothing function

Question

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[[1]];
sy = dim[[2]];

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).

Answer 1

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["https://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]];
ImageAdd[ImageConvolve[redPixels, box], 
  ImageConvolve[blackPixels, box]] // ImageAdjust

Answer 2

You can use ColorReplace:

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