J'essaie de faire une correspondance d'histogramme en utilisant Python pour améliorer le processus de mosaïquage de plusieurs rasters qui se chevauchent. Je fonde mon code sur celui trouvé sur:
http://www.idlcoyote.com/ip_tips/histomatch.html
À ce jour, j'ai réussi à découper la zone de chevauchement de deux rasters adjacents et à aplatir le réseau.
j'ai donc deux tableaux 1 dimension de la même longueur.
J'ai ensuite écrit le code suivant basé sur celui trouvé sur le site Web ci-dessus. Dans le code montré, j'ai substitué deux très petits jeux de données aux images gd et 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()
Mon programme trace les histogrammes et les distributions de fréquences cumulées avec succès ... et je pensais que j'avais la partie d'obtenir la fonction de transformation 'z' correcte .... mais ensuite quand j'utilise la fonction de distribution 'fi' sur le 'bd_hist' pour essayer de le faire correspondre au jeu de données gd, tout se passe en forme de poire.
Je ne suis pas mathématicien et il est fort probable que j'ai oublié quelque chose d'assez évident.
cdf_bd = np.cumsum(bd_hist) / float(np.sum(bd_hist))
Réponses:
Bien que je ne puisse pas commenter l'implémentation suggérée, vous voudrez peut-être vérifier une implémentation existante de la correspondance d'histogramme effectuée pour GRASS GIS 7 (ici un addon):
https://trac.osgeo.org/grass/browser/grass-addons/grass7/imagery/i.histo.match
Pour le manuel et un exemple, voir
http://grass.osgeo.org/grass70/manuals/addons/i.histo.match.html
Le code est publié sous la licence GPL2 +.
la source
Comme un fudge sauvage; Je ne suis pas sûr que vous ayez besoin d'un PDF si vous avez des données de comptage dans les catégories ...
Pourriez-vous convertir le nombre de chaque valeur pour chaque histogramme différent en valeurs XY, puis utiliser une sorte d'indicateur de régression pour vérifier cette correspondance? Autrement dit, pour deux histogrammes parfaitement identiques, une analyse de corrélation fournirait et R au carré de 1,0.
la source
certains exemples de données seraient bien car ils peuvent varier de sat à sat. voici un script simple que j'ai fait pour essayer d'égaliser les histogrammes:
https://github.com/rupestre-campos/histogram_equalize
Vous pouvez peut-être obtenir un aperçu.
Il calcule également le cdf comme vous le faites, mais comme je l'ai essayé, il deviendra fou si vous calculez bande par bande, vous devez donc considérer l'ensemble du raster.
On dirait que vous perdez l'équilibre de référence des couleurs et le profil spectral. Il faut également ne pas compter aucun pixel de données, puis le supprimer du nombre total de pixels de l'image afin de calculer correctement le pdf.
Après quelques tests, j'ai aimé les résultats visuels en utilisant toute l'approche raster de 3-4 bandes Landsat8 et en convertissant de 16 bits à 8 bits 0-255.
la source