Histograma Combinando usando Python para melhorar o processo de mosaico de várias rasters sobrepostas?

Eu estou tentando fazer a correspondência de histograma usando Python para melhorar o processo de mosaico de várias rasters sobrepostas. Estou baseando meu código no encontrado em:

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

Até o momento, consegui recortar a área sobreposta de duas rasters adjacentes e achatar a matriz.

então eu tenho duas matrizes unidimensionais do mesmo comprimento.

Escrevi o código a seguir com base no código encontrado no site acima. No código mostrado, substituí dois conjuntos de dados muito pequenos pelas imagens gd e 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()

Meu programa plota os histogramas e distribuições de frequência cumulativa com êxito ... e eu pensei que tinha o papel de corrigir a função de transformação 'z' .... mas depois, quando eu uso a função de distribuição 'fi' no 'bd_hist' para tentar combiná-lo com o conjunto de dados gd, tudo fica em forma de pêra.

Não sou matemático e é muito provável que tenha esquecido algo bastante óbvio.

Embora eu não possa comentar sobre a implementação sugerida, convém verificar uma implementação existente da correspondência de histograma feita para o GRASS GIS 7 (aqui um complemento):

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

Para o manual e um exemplo, consulte

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

O código é publicado sob a licença GPL2 +.

Como um fudge selvagem; Não tenho certeza se você precisa de um PDF se tiver dados de contagem em categorias ...
Você pode converter as contagens de cada valor para cada histograma diferente em valores XY e usar algum tipo de indicador de regressão para verificar essa correspondência? Ou seja, para dois histogramas perfeitamente idênticos, uma análise de correlação forneceria e R ao quadrado de 1,0.

alguns dados de amostra seriam bons, pois podem variar de sab para sab. Aqui está um script simples que eu fiz na tentativa de equalizar histogramas:

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

Talvez você possa ter uma ideia.

Ele também calcula o cdf como você faz, mas como eu tentei, fica louco se você calcular banda por banda, então você deve considerar toda a varredura.

Parece que você perde o equilíbrio de referência de cores e o perfil espectral. Além disso, é necessário não contar pixels de dados, é necessário remover o número total de pixels da imagem para calcular o pdf corretamente.

Após alguns testes, gostei dos resultados visuais, usando toda a abordagem de varredura para 3-4 bandas do Landsat8 e convertendo de 16 bits para 8 bits no intervalo 0-255.