Estoy tratando de hacer una coincidencia de histograma con Python para mejorar el proceso de mosaico de múltiples rásteres superpuestos. Estoy basando mi código en el que se encuentra en:
http://www.idlcoyote.com/ip_tips/histomatch.html
Hasta la fecha, he logrado recortar el área superpuesta de dos rásteres adyacentes y aplanar la matriz.
entonces tengo dos matrices de 1 dimensión de la misma longitud.
Luego escribí el siguiente código basado en el que se encuentra en el sitio web anterior. En el código que se muestra, he sustituido dos conjuntos de datos muy pequeños para las imágenes gd y bd.
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()
Mi programa traza los histogramas y las distribuciones de frecuencia acumulativa con éxito ... y pensé que tenía la parte de obtener la función de transformación 'z' correcta ... pero luego cuando uso la función de distribución 'fi' en el 'bd_hist' para intentar hacer coincidirlo con el conjunto de datos de gd, todo tiene forma de pera.
No soy matemático y es muy probable que haya pasado por alto algo bastante obvio.
cdf_bd = np.cumsum(bd_hist) / float(np.sum(bd_hist))Respuestas:
Si bien no puedo comentar sobre la implementación sugerida, es posible que desee verificar una implementación existente de coincidencia de histograma realizada para GRASS GIS 7 (aquí un complemento):
https://trac.osgeo.org/grass/browser/grass-addons/grass7/imagery/i.histo.match
Para el manual y un ejemplo, vea
http://grass.osgeo.org/grass70/manuals/addons/i.histo.match.html
El código se publica bajo la licencia GPL2 +.
fuente
Como un dulce de azúcar; No estoy seguro de que necesite un PDF si tiene datos de conteo en categorías ...
¿Podría convertir los conteos de cada valor para cada histograma diferente en valores XY y luego usar algún tipo de indicador de regresión para verificar esa coincidencia? Es decir, para dos histogramas perfectamente idénticos, un análisis de correlación proporcionaría y R al cuadrado de 1.0.
fuente
algunos datos de muestra serían buenos, ya que pueden variar de sat a sat. Aquí hay un script simple que hice en un intento de ecualizar histogramas:
https://github.com/rupestre-campos/histogram_equalize
Tal vez puedas obtener alguna idea.
También calcula el cdf como lo hace, pero como he intentado, se volverá loco si calcula banda por banda, por lo que debe considerar todo el ráster.
Parece que pierde el balance de referencia de color y el perfil espectral. También es necesario no contar ningún píxel de datos, debe eliminarse del número total de píxeles de la imagen para calcular el pdf correcto.
Después de algunas pruebas, me gustaron los resultados visuales utilizando el enfoque de trama completo para las bandas Landsat8 de 3-4 y la conversión de 16 bits a 8 bits de rango 0-255.
fuente