# Equal probability quantizer in Python

I am trying to apply Haralick textures to a SAR image (float32). As far as I know, the image first needs to be quantized to a reasonable bit depth prior to calculate the co-ocurrence matrix. In the original paper, Haralick proposes using a "equal probability quantizer". As far as I know is the algorithm employed by Mathlab "imquantize". I have been looking for a pre-existing python/numpy implementation (no success).

I also have reading pre-existing java code from ESA SNAP. As far I know histogram calculations would be straightforward with numpy, but if I understand correctly, every pixel is assigned a "level/histogram bin" using bisection search, so I do not see any straightforward way to properly vectorize the code.

EDITED after this point ===================

I have made simple implementation in Python + Numpy based on the pre-existing java code from ESA SNAP, and the kind advice of Marcus Müller.

import numpy as np

def eqProbQuant(image, levels=32):
# Sort the pixels by value
sorted_image = np.sort(image)
# Get the pixel count
pixel_count = sorted_image.shape[0]
# Get the number of pixels per bin
samples_per_bin = int(pixel_count/levels)
# Get locations where the bin would change
edge_samples = np.arange(levels+1) * samples_per_bin
# Get the values at those locations (bin edges)
bin_edges = sorted_image[edge_samples]
# Use the values to apply quantization
quantized = np.digitize(image, bin_edges)

return quantized


If possible, I would like to know if this code would apply the algorithm correctly (I have very little experience with signal processing)

• Hi! Welcome here. Asking for code written to your specification is off-topic here, but I'm 100% certain you can rephrase your question to describe the signal processing problem you're having implementing that quantizer yourself (it really shouldn't be hard at all). – Marcus Müller Aug 1 at 9:48
• Thanks a lot for your kind advice, I tried to rephrase the question in a more polite way. I want to apologise for the original phrasing. I come from biological sciences, and the topic escapes me a bit. Do you see any way it could be further improved? – opengisapprendice Aug 1 at 10:33
• I shouldn't have used "rephrasing". I meant to say "you can pose a different problem than asking for code that someone else wrote" :) Anyway, the question is fine, let me try to answer it as is – Marcus Müller Aug 1 at 11:02
• Re-focused again (this time without nagging for code!), just asking if the code would properly apply the algorithm. I would greatly appereciate it if you could take a look to know if this would work as intended. If not, I understand, you already helped enough – opengisapprendice Aug 2 at 10:22

There is any way to assign a a "level/histogram bin" to each pixel without iterating through the entire image?

(emphasis mine)

No, because what you describe is an operation that needs to be applied to each pixel.

Also, a histogram needs to be calculated taking all pixels into account (that's the definition of a histogram!).

Anyway, that's not a problem at all – in SAR processing, this operation is definitely of "a lesser concern" when it comes to computational intensity.

1. Calculate a list of pixel values sorted by intensity
2. assign value ranges based on that

# Sort the pixels

That's inherently an operation that can't be fully parallelized – in the end, there is one list for the whole image. However, you can parallelize it for the most part: Let your multiple processors (GPU threads, …) sort disjunct parts of the image, then merge the sorted lists in order (and that can be done in cascaded ways, too). Pseudocode:

N := numbers of pixels
I := image considered as linear array of pixels
O := output list of sorted pixel values, same length as I
T := temporary list of list of pixels
P := numbers of processing cores

FOR p IN P PARALLEL:
T[p][:] := sort(I[p*(N/P) … (p+1)*(N/P)])

merge(T[0…(P-1)])


# Quantize

Then, just take the sorted list of values, divide into the desired number of bins of equal size, note down the "starting" value of each bin (i.e. the values at O[0], O[N/n_bins], O[2*N/n_bins]...).

Then, (again, in parallel if desired – there's no data interdependence between pixels) move through your image and replace the pixel values by the bin number.