Fast Way to Remove Isolated Pixels in a Highly Quantized Image

Question

I have an image that I've quantized like so:

And I would like to get rid of pixels that aren't much like their neighbors (basically do a low-pass filter). The goal would be to get rid of isolated pixels of a particular value or narrow strips of such pixels.

I've tried doing a Gaussian blur and then re-quantizing, but it still leaves behind narrow strips, e.g.

One idea I had was for each pixel to compute the most common pixel among it's neighbors in region

this mostly works

But my naive implementation takes O(num_pixels*radius**2)

def knn(img,radius=4):

    imgOut=img.copy()

    for i in trange(img.size[0]):
        for j in range(img.size[1]):
          counts={}
          for k in range(-radius,radius+1):
                for l in range(-radius,radius+1):
                    if i+k>=0 and i+k<img.size[0] and j+l>=0 and j+l<img.size[1]:
                        v=img.getpixel((i+k,j+l))
                        counts[v]=counts.get(v,0)+1
          l=list(counts.items())
          l.sort(key=lambda x:x[1])
          imgOut.putpixel((i,j),l[-1][0])
    
    return imgOut

And it doesn't even get rid of all the isolated pixels, since small regions shrink, creating more isolated pixels.

Is there a faster/better way to achieve the same effect?

I feel like there ought to be some kind of discrete Fourier transform + low pass method, but I'm not sure what it would be.

=============

Update:

scipy.ndimage.median_filter

Is almost what I want, but it still leaves behind narrow edges (like guassian blur did). This is because if you have (for example) an area that looks like [0,0,0,1,2,2,2] the median will be 1, leaving behind a narrow strip.

I don't think this is the most elegant way, but I get around this by doing multiple passes and randomizing the labels at each pass

def smooth_image(img,r=2,k=16):
    a= np.asarray(img5)
    aa=np.copy(a)
    for i in trange(k):
        v=np.random.permutation(NUM_COLORS)
        ip=np.argsort(v)
        #shuffle labels
        aa=v[aa]
        aa=scipy.signal.medfilt2d(aa,2*r+1)
        #unshuffle
        aa=ip[aa]    
        bmp = Image.fromarray(aa.astype(np.uint8))
        bmp.putpalette(img5.palette) 
    bmp = Image.fromarray(aa.astype(np.uint8))
    bmp.putpalette(img.palette)     
    return bmp

which produces the follow result

Answer 1

This is something you'd typically accomplish with an opening and/or a closing filter. Those filters however are biased (towards dark and light regions, respectively). To avoid the bias, you could apply an alternating sequential filter, which is a sequence of openings and closings of increasing size.

Still, those filters will also leave narrow regions of the intermediate color.

I came up with this "multi-region" version of the closing. It basically applies an erosion to each color independently, then fills in the area in between the remaining regions by growing those regions until they touch each other. I'm using DIPlib (disclaimer: I'm an author).

The first part is just converting the PNG attached to the question into a labeled image, one that has pixels with values 1 through 4:

import diplib as dip
import numpy as np

original = dip.ImageRead('/Users/cris/Downloads/akVCv.png')
values = np.array([45,135,180,216])
input = dip.MultipleThresholds(original(1), (values[:3] + values[1:]) / 2) + 1

values are the four intensities of the green channel.

We can now apply the suggested operation:

ero = 0
for ii in range(1,5):
   ero = ero + dip.Erosion(input == ii, 11) * ii
ero.Convert('UINT8')
out = dip.GrowRegions(ero)

ero is the result of the erosion, we erode the binary image obtained for each of the colors. 11 is the diameter of the circular structuring element. This is a large value, which I picked on purpose to make the simplification of the image more obvious. A smaller value might be better in practice. In the display below, the darkest color represents "no color" (the image has values of 0 here). Note how the highly textured area of the top-left quadrant is mostly "no color" now, which is not ideal.

out is the result, where the blank areas in ero are filled in with the nearest color. The top-left quadrant is the region that looks least like the original image. One way to avoid this effect could be to apply this filter two or three times, with increasing structuring element size.