Histogramm-Matching mit Python zur Verbesserung des Mosaikprozesses mehrerer überlappender Raster?

11

Ich versuche, mit Python einen Histogrammabgleich durchzuführen, um den Mosaikprozess mehrerer überlappender Raster zu verbessern. Ich stütze meinen Code auf den folgenden:

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

Bisher habe ich es geschafft, den überlappenden Bereich zweier benachbarter Raster zu beschneiden und das Array zu reduzieren.

Ich habe also zwei eindimensionale Arrays gleicher Länge.

Ich habe dann den folgenden Code geschrieben, der auf dem auf der obigen Website gefundenen basiert. In dem gezeigten Code habe ich die gd- und bd-Bilder durch zwei sehr kleine Datensätze ersetzt. 

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()

Mein Programm zeichnet die Histogramme und kumulativen Häufigkeitsverteilungen erfolgreich auf ... und ich dachte, ich hätte den Teil, die Transformationsfunktion 'z' korrekt zu machen ... aber dann, wenn ich die Verteilungsfunktion 'fi' auf der 'bd_hist' verwende. Um zu versuchen, es mit dem GD-Datensatz abzugleichen, wird alles birnenförmig.

Ich bin kein Mathematiker und es ist sehr wahrscheinlich, dass ich etwas ziemlich Offensichtliches übersehen habe.

python raster histogram Becky
quelle
Ich weiß nicht viel über Histogramm-Matching, aber müssen Ihre CDFs (per Definition) zu 1 summieren? cdf_bd = np.cumsum(bd_hist) / float(np.sum(bd_hist))
Jeff G

Antworten:

1

Als wilder Fudge; Ich bin mir nicht sicher, ob Sie ein PDF benötigen, wenn Sie Zähldaten in Kategorien haben ...
Könnten Sie die Zählungen jedes Werts für jedes unterschiedliche Histogramm in XY-Werte konvertieren und dann eine Art Regressionsindikator verwenden, um diese Übereinstimmung zu überprüfen? Das heißt, für zwei vollkommen identische Histogramme würde eine Korrelationsanalyse ein R-Quadrat von 1,0 liefern.

Mox
quelle
0

Einige Beispieldaten wären nett, da sie von Sat zu Sat variieren können. Hier ist ein einfaches Skript, das ich erstellt habe, um Histogramme auszugleichen:

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

Vielleicht können Sie einen Einblick bekommen.

Es berechnet auch cdf wie Sie, aber wie ich versucht habe, wird es verrückt, wenn Sie Band für Band berechnen, also haben Sie das ganze Raster in Betracht gezogen.

Sieht so aus, als ob Sie die Farbreferenzbalance und das Spektralprofil verlieren. Es besteht auch die Notwendigkeit, keine Datenpixel zu zählen, sondern diese dann aus der Gesamtbildpixelzahl zu entfernen, um das PDF korrekt zu berechnen.

Nach einigen Tests gefielen mir die visuellen Ergebnisse unter Verwendung des gesamten Raster-Ansatzes für Landsat8 mit 3-4 Bändern und der Konvertierung von 16 Bit auf 8 Bit 0-255.

ckc
quelle