# Locate Non Homogenous Areas in an Image

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?

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.

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

• What do you mean by "Variance by Super Pixel"? After I find the superpixels say I have 10 superpixels, then should I calculate a sliding window say 5x5 variance on super pixels too? Is it a two layered process? Isn't a "local variance" only on source image enough? I just calculated a super pixel of lena. Should I calculate a variance on this too? imgur.com/a/3rnkSUe – Lyrk Jun 7 at 20:31
• I mean after calculating super pixels and their average colors, should I apply variance on the resultant image too? (I want to choose the points based on most non-homogenous to least non-homogenous) @Royi – Lyrk Jun 7 at 20:38
• @Lyrk, I think you marked another answer as answering your question. From my experience, using wavelets is the 90's answer to this. We have more modern approaches as I wrote above. The code shows exactly how to calculate the variance per Super Pixel. Template matching of Wavelets is sub optimal to this problem. You can do better things with Wavelets. But a simpler and more elegant would be to chose the correct tool to begin with. – Royi Jun 9 at 3:45
• I mean is this the process? 1- Find superpixels 2- Find mean color of superpixels 3- Apply variance on result of 2. – Lyrk Jun 9 at 7:21
• @Lyrk, Indeed. I will update the answer above to be explicit on the process. – Royi Jun 9 at 16:09

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.

### Library

Currently in dev stage, you can obtain latest version (from here) via:

pip install git+https://github.com/kymatio/kymatio.git@refs/pull/674/head


or clone dev branch that's without kymatio.visuals. Modulus + lowpassing + subsampling is built-in, and center frequency is tied to scale.

### Code

import numpy as np
import PIL
from kymatio.numpy import Scattering2D
from kymatio.visuals import imshow

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"
#%%
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)

• @Lyrk The visuals are only available from latest version which you can install with code I provided under "Library": pip install git+https://github.com/kymatio/kymatio.git@refs/pull/674/head – OverLordGoldDragon Jun 7 at 11:40
• @Lyrk Unfortunately subsampling is built in, but I've made a modified code you can use that's also without lowpassing - paste this in kymatio/scattering2d/core/scattering2d.py - and updated answer with its results for my example. – OverLordGoldDragon Jun 7 at 13:02
• @Royi I don't follow. Wavelets can calibrate to detect variation of any rate, in any orientation, at any scale. They can also extract higher order variations that other methods can't. Everything is captured. – OverLordGoldDragon Jun 7 at 20:33
• @OverLordGoldDragon, Yet the output is local. Namely they don't really have any real advantage for this task over the local mean. Better approaches are based on local segmentation. Something like extracting features of Super Pixels. Wavelets are not an optimized tool for this task. – Royi Jun 8 at 3:22
• @Royi Local mean can't target rate or orientation, or higher order variation. Wavelets enable specifying desired homogeneity template, with controlled tradeoff between robustness and localization. I don't deny Super Pixel etc may do better. Wavelets may not be optimal, but have their advantages. – OverLordGoldDragon Jun 8 at 16:37