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.
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 but if anyone can offer any help that would be absolutely amazing. Thank you for your time.
Becky
أكثر...
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.
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 but if anyone can offer any help that would be absolutely amazing. Thank you for your time.
Becky
أكثر...