9

私は python/numpy/scipy を使用して、地形の側面と勾配に基づいて 2 つのデジタル標高モデル (DEM) を整列させるためのこのアルゴリズムを実装しています。

「氷河の厚さの変化を定量化するための衛星標高データ セットの共同登録とバイアス補正」、C. Nuth および A. Kääb、doi:10.5194/tc-5-271-2011

私はフレームワークをセットアップしましたが、scipy.optimize.curve_fit によって提供される適合の質が悪いです。

def f(x, a, b, c):
    y = a * numpy.cos(numpy.deg2rad(b-x)) + c
    return y

def compute_offset(dh, slope, aspect):
    import scipy.optimize as optimization

    idx = random.sample(range(dh.compressed().size), 10000)
    xdata = numpy.array(aspect.compressed()[idx], float)
    ydata = numpy.array((dh/numpy.tan(numpy.deg2rad(slope))).compressed()[idx], float)

    #Generate synthetic data to test curve_fit
    #xdata = numpy.arange(0,360,0.01)
    #ydata = f(xdata, 20.0, 130.0, -3.0) + 20*numpy.random.normal(size=len(xdata))

    print xdata
    print ydata

    x0 = numpy.array([0.0, 0.0, 0.0])

    fit = optimization.curve_fit(f, xdata, ydata, x0)[0]
    #optimization.leastsq(f, x0[:], args=(xdata, ydata))
    genplot(xdata, ydata, fit)
    return fit

def genplot(x, y, fit):
    a = (numpy.arange(0,360))
    f_a = f(a, fit[0], fit[1], fit[2])
    idx = random.sample(range(x.size), 10000)
    plt.figure()
    plt.xlabel('Aspect (deg)')
    plt.ylabel('dh/tan(slope) (m)')
    plt.plot(x[idx], y[idx], 'r.')
    plt.axhline(color='k')
    plt.plot(a, f_a, 'b')
    plt.ylim(-80,80)
    plt.show()

#Input DEMs
dem1_fn = sys.argv[1]
dem2_fn = sys.argv[2]

dem1_ds = gdal.Open(dem1_fn, gdal.GA_ReadOnly)
dem2_ds = gdal.Open(dem2_fn, gdal.GA_ReadOnly)

#Extract band 1 from each dataset as masked array using internal nodata value
dem1 = getperc_new.gdal_getma(dem1_ds, 1)
dem2 = getperc_new.gdal_getma(dem2_ds, 1)

#Produce slope and aspect maps using gdaldem and load into masked arrays
dem1_slope = gdaldem_slope(dem1_fn)
dem1_aspect = gdaldem_aspect(dem1_fn)

#Compute common mask and apply to all products
common_mask = dem1.mask + dem2.mask + dem1_slope.mask + dem1_aspect.mask

diff_euler = numpy.ma.array(dem2-dem1, mask=common_mask)
dem1_slope.__setmask__(common_mask)
dem1_aspect.__setmask__(common_mask)

#Compute relationship between elevation difference, slope and aspect
fit = compute_offset(diff_euler, dem1_slope, dem1_aspect)
print fit

これは、最初は約 200 万ポイントで構成されている私のデータに適合しますが、テスト/プロットの目的でランダムにサンプリングしました。

[ -14.9639559 216.01093596 -41.96806735]

plot_one

適切な適合を示すデータはたくさんありますが、curve_fit の結果はよくありません。合成データで実行すると、次のようにうまく適合します。

元の入力パラメータ [20.0、130.0、-3.0]

Curve_fit の結果 [-19.66719631 -49.6673076 -3.12198723]

plot_two

これがマスクされた配列の使用、curve_fit の制限と関係があるのか​​ 、それとも単純なものを見落としているだけなのかはわかりません。提案をありがとう。

==========================

編集 2013 年 9 月 4 日 16:30 PDT

@Evert などで示唆されているように、問題は間違いなく外れ値に関連していました。外れ値を取り除いた後、より良い適合を得ることができました。私の古いコードを見ると、各側面の絶対偏差の中央値を計算し、フィッティングの前に 2*mad 以外のものをすべて削除したようです。

2012 年 11 月にいくつかの追加プロットを生成しました。

ヒストグラム

ここに画像の説明を入力

しかし、これらをもう一度見てみると、異なる入力データに対して生成されたものであることがほぼ確実です。今見つけられるのはこれだけなので、偏ったサンプリングの例としてここに含めます。この DEM アライメントの方法は、このような場合には失敗するはずです。また、scipy のカーブ フィッティング機能とは何の関係もありません。

最終的に、2 つのマスクされた 2D numpy 配列の正規化された相互相関、サブピクセルの調整、および垂直オフセットの除去を含む、位置合わせのための別のアプローチを開発しました。それはより速く、一貫してより良い結果を提供します。そのアプローチでさえ、 NASA Ames Stereo Pipelineの一部として Oleg Alexandrov によって開発された Iterative Closest Point (ICP) ツール (pc_align) に取って代わられました。

すべての回答に感謝します。この質問を放棄したことをお詫びします。

4

1 に答える 1

1

位相オフセットのある正弦波を得ようとしているだけなら、非線形フィットは必要ありません。

a * sin(x - b) + cオフセットのある正弦は、正弦と余弦の適切な組み合わせとして記述できるため、これをで置き換えることがa * sin(x) + b * cos(x) + cできます (フーリエ変換のような「フェーザー加算」)。

それでも同じ結果が得られる場合、問題は「非線形」フィットではありません。

于 2012-11-03T16:52:00.493 に答える