Histogram Matching using Python to improve mosaicking process of multiple overlapping rasters?

Question

I am trying to do histogram matching using Python to improve the mosaicking process of multiple overlapping rasters. I am basing my code on that found at:

http://www.idlcoyote.com/ip_tips/histomatch.html

To date, I have managed to clip the overlapping area of two adjacent rasters and flatten the array.

so I have two 1 dimensional arrays of the same length.

I have then written the following code based on that found at the above website. In the code shown I have substituted two very small datasets for the gd and bd images.

import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

bins = range(0,100, 10)

gd_hist = [1,2,3,4,5,4,3,2,1]

bd_hist = [2,4,6,8,10,8,6,4,2]

nPixels = len(gd_hist)

# here we are creating the cumulative distribution frequency for the bad image
cdf_bd = []
for k in range(0, len(bins)-1):
    b = sum(bd_hist[:k]) 
    cdf_bd.append(float(b)/nPixels)

# here we are creating the cumulative distribution frequency for the good image
cdf_gd = []
for l in range(0, len(bins)-1):
    g = sum(gd_hist[:l])
    cdf_gd.append(float(g)/nPixels) 


# we plot a histogram of the number of 
plt.plot(bins[1:], gd_hist, 'g')
plt.plot(bins[1:], bd_hist, 'r--')
plt.show()        

# we plot the cumulative distribution frequencies of both images
plt.plot(bins[1:], cdf_gd, 'g')
plt.plot(bins[1:], cdf_bd, 'r--')
plt.show()

z = []
# loop through the bins
for m in range(0, len(bins)-1):

    p = [cdf_bd.index(b) for b in cdf_bd if b < cdf_gd[m]] 
    if len(p) == 0:
        z.append(0)
    else:
        # if p is not empty, find the last value in the list p
        lastval = p[len(p)-1]

        # find the bin value at index 'lastval'
        z.append(bins[lastval])

plt.plot(bins[1:], z, 'g')
plt.show()

# look into the 'bounds_error'
fi = interp1d(bins[1:], z, bounds_error=False, kind='cubic')  
plt.plot(bins[1:], gd_hist, 'g')
plt.show
plt.plot(bins[1:], fi(bd_hist), 'r--')
plt.show()

My program plots the histograms and cumulative frequency distributions successfully...and I thought that I had the part of getting the transformation function 'z' correct....but then when I use the distribution function 'fi' on the 'bd_hist' to try to match it to the gd dataset it all goes pear-shaped.

I am not a mathematician and it is highly likely I have overlooked something fairly obvious.

Answer 1

While I cannot comment on the suggested implementation, you may want to check an existing implementation of histogram matching done for GRASS GIS 7 (here an addon):

https://trac.osgeo.org/grass/browser/grass-addons/grass7/imagery/i.histo.match

For the manual and an example, see

http://grass.osgeo.org/grass70/manuals/addons/i.histo.match.html

The code is published under the GPL2+ licence.

Answer 2

As a wild fudge; I'm not sure you need a PDF if you've got count data in categories...
Could you convert the counts of each value for each different histogram into XY values, and then use some sort of regression indicator to check that match? Ie, for two perfectly identical histograms, a correlation analysis would provide and R squared of 1.0.

Answer 3

some sample data would be nice as it may vary from sat to sat. here is one simple script I've made in a try to equalize histograms:

https://github.com/rupestre-campos/histogram_equalize

Maybe you can get some insight.

It also computes cdf as you do, but as I've tried it will go crazy if you compute band-per-band, so you got consider the whole raster.

Looks like you loose color reference balance and spectral profile. Also there is the need to not count no data pixels, gotta remove then from the total image pixel number in order to compute pdf correct.

After some testing I liked the visual results using the whole raster approach to 3-4 bands Landsat8 and converting from 16bit to 8bit 0-255 range.