Locate Non Homogeneous Areas in an Image

Question

I need to choose some points on an image. These points should be chosen more where there is lots of color changes, transitions and variations. Which techniques can I use for determining where most color changes and transitions occur in an image?

Answer 1

In general, the approach to take, is to have a local feature which has high value for such areas in the image.

There are many approaches to shape such a feature.
Probably the easiest one would be by local variance.

I tried 3 different approaches to this:

1. Local Variance by a Filter.
2. Local Variance of a Super Pixel.
3. Using the Weak Texture from Noise Level Estimation from a Single Image (By Masayuki Tanaka).

I applied it on the Lenna image and got:

Higher values shows non homogenous area.
It seems that for high SNR images you can work with the local Variance. But the super pixel approach seems to be more robust.

In my opinion the Super Pixel result is the best of all 3.
It was achieved by:

1. Apply Super Pixel based Segmentation (SLIC Based).
2. Per Label Mean
   Calculate the mean value of each Super Pixel by the indices of each super pixel.
3. Per Label Variance
   Calculate the variance of each super pixel using only its pixels.

In a more general form if you look for features for homogeneity and use them to find homogenous regions and then select the inverse. Those are very popular in segmentation.

Another approach would be using more advanced features such as:

1. BRISK Feature.
2. FAST Feature.
3. HOG Feature.
4. MSER Feature.

Then count how many of those are found within each Super Pixel. A Super Pixel with more features will be less homogenous.

The full code is available on my StackExchange Signal Processing Q75536 GitHub Repository (Look at the SignalProcessing\Q75536 folder)..

Update: today I encountered Robust Segmentation Free Algorithm for Homogeneity Quantification in Images.

Answer 2

2D wavelet transform is well suited. It's an extension of 1D CWT where we correlate wavelets of different center frequencies and "scales" (widths in time domain).

Wavelets can be calibrated to detect

• fast or slow variations
• over small, localized or large, spread out parts of image

The output is a 3D array indexed as:

• x: x-coordinate of wavelet centering over image
• y: y-coordinate of wavelet centering over image
• n: wavelet (that produced the output), controlling frequency (rate of variation), width (spatial extent of variation), and angle (orientation of variation)

Then one can find e.g. "greatest variation over a 2cm x 2cm region" by indexing with a proper unit conversion and taking argmax over the 2D slice; each slice is an "intensity heatmap" of variations per n. Outputs are made robust to noise at expense of spatial localization by lowpassing the modulus of output.

Example

The image (source):

Example wavelets (j = scale index (width + frequency), theta = angle index):

Grouped by scale

scale=0 captures fast variations over small intervals, while scale=3 captures slow variations over large regions

Grouped by angle

Orientation information is obfuscated by averaging over all the scales. Suppose we wish to detect this transition:

Then we want a large scale wavelet that's steeply oriented. And indeed:

Higher-order variation

Obtained by taking wavelet transform of wavelet transform, capturing variations of variations:

(scale=0 is missing for theoretical reasons) Comparing against first-order, we see bottom-left is more intense, which follows from colors shifting at a changing rate (faster toward fractal singularity).

Prioritizing the fast

To preserve fastest variations, set J=1, which omits low-freq wavelets and makes lowpass narrow (also reducing subsampling):

Too many outputs?

Just use those you need. Inspect the wavelets used, and keep a subset (or just one); examples in code.

J=3 no subsampling / lowpassing

Can't do with library code, see comment below answer.

---

Lena example

Code; can probably improve

Going further

A conv-net can be trained on top of scattering features with segmentation objective; wavelet scattering attains SOTA on many benchmarks in limited data settings.

More examples -- scattering lecture.

Code

import numpy as np
import PIL
from kymatio.numpy import Scattering2D
from kymatio.visuals import imshow
def load_rgb(img):
    img = np.array(PIL.Image.open(path).convert("RGB")).astype('float64')
    img /= np.abs(img).max(axis=(0, 1))
    return img
def group_by_angle(out, J, L):
    n_S1 = J * L  # number of first-order coeffs
    S1_all = out[1:n_S1 + 1]
theta_all = list(range(L))
S1_slices = {k: [] for k in theta_all}

i = 0
for p1 in S.psi:
    S1_slices[p1['theta']].append(S1_all[i][None])
    i += 1

for theta in S1_slices:
    S1_slices[theta] = np.vstack(S1_slices[theta]).mean(axis=0)[None]
S1_theta = np.vstack(list(S1_slices.values()))
return S1_theta


def group_by_scale(out, J, L):
    n_S1 = J * L  # number of first-order coeffs
    S1_all = out[1:n_S1 + 1]
    S1_scale = S1_all.reshape(L, -1, *out.shape[-2:], order='F').mean(axis=0)
    return S1_scale
def viz_by_angle(S1_theta, L, ticks=0):
    for theta in range(L):
        th = int(L - L/2 - 1) - theta
        angle = 90 * th / L
        imshow(S1_theta[theta], title=f"angle={angle}",
               w=.5, h=.5, abs=1, ticks=ticks)
def viz_by_scale(S1_scale, J, ticks=0, second_order=False):
    J = J - 1 if second_order else J
    for j in range(J):
        j_title = j + 1 if second_order else j
        imshow(S1_scale[j], title=f"scale={j_title}",
               w=.5, h=.5, abs=1, ticks=ticks)
def viz_all_first_order(out, J, L, ticks=0):
    n_S1 = J * L  # number of first-order coeffs
    S1_all = out[1:n_S1 + 1]
for i, o in enumerate(S1_all):
    j, theta = S.psi[i]['j'], S.psi[i]['theta']
    th = int(L - L/2 - 1) - theta
    angle = 90 * th / L
    imshow(o, title=f"angle={angle}, scale={j}",
           abs=1, w=.5, h=.5, ticks=ticks)


def group_by_angle_second_order(out, J, L):
    n_S1 = J * L
    S2_all = out[n_S1 + 1:]
    theta_all = list(range(L))
    S2_slices = {k: [] for k in theta_all}
i = 0
for p1 in S.psi:
    j1 = p1['j']
    for n2, p2 in enumerate(S.psi):
        j2 = p2['j']
        if j2 <= j1:
            continue
        S2_slices[p2['theta']].append(S2_all[i][None])
        i += 1

for theta in S2_slices:
    if S2_slices[theta] != []:
        S2_slices[theta] = np.vstack(S2_slices[theta]).mean(axis=0)[None]

S2 = np.vstack([arr for arr in S2_slices.values() if len(arr) != 0])
return S2


def group_by_scale_second_order(out, J, L):
    n_S1 = J * L
    S2_all = out[n_S1 + 1:]
    j_all = list(range(J))
    S2_slices = {k: [] for k in j_all}
i = 0
for p1 in S.psi:
    j1 = p1['j']
    for n2, p2 in enumerate(S.psi):
        j2 = p2['j']
        if j2 <= j1:
            continue
        S2_slices[p2['j']].append(S2_all[i][None])
        i += 1

for j in S2_slices:
    if S2_slices[j] != []:
        S2_slices[j] = np.vstack(S2_slices[j]).mean(axis=0)[None]

S2 = np.vstack([arr for arr in S2_slices.values() if len(arr) != 0])
return S2


#%%
path = r"C:\Desktop\colors.png"
img = load_rgb(path)
#%%
imshow(img, w=.6, h=.6, title="%s x %s image" % img.shape[:2])
#%%
J, L = 4, 8  # largest scale; number of angles
S = Scattering2D(shape=img.shape[:2], L=L, J=J)
#%%
take 3 wavelet transforms on each channel and average
out = np.mean([S(im) for im in img.transpose(-1, 0, 1)], axis=0)
#%%
S1_scale = group_by_scale(out, J, L)
viz_by_scale(S1_scale, J)
#%%
S1_theta = group_by_angle(out, J, L)
viz_by_angle(S1_theta, L)
#%%
S2_theta = group_by_angle_second_order(out, J, L)
viz_by_angle(S2_theta, L)
#%%
S2_scale = group_by_scale_second_order(out, J, L)
viz_by_scale(S2_scale, J, second_order=True)
#%%
viz_all_first_order(out, J, L)