7

Numpy を使用した画像処理、具体的には実行中の標準偏差ストレッチの実行に取り組んでいます。これは X 列を読み取り、Std を見つけます。パーセンテージの線形ストレッチを実行します。次に、列の次の「グループ」に反復し、同じ操作を実行します。入力画像は 1 GB、32 ビット、シングル バンド ラスターであり、処理に非常に長い時間 (数時間) かかります。以下はコードです。

おそらくボトルネックが発生している場所である 3 つの入れ子になった for ループがあることに気付きました。「ボックス」で画像を処理する場合、つまり、[500,500] の配列をロードし、画像処理時間を繰り返すのは非常に短いです。残念ながら、カメラのエラーにより、バンディングを避けるために非常に長いストリップ (52,000 x 4) (y,x) で反復する必要があります。

これをスピードアップするための提案は大歓迎です：

def box(dataset, outdataset, sampleSize, n):

    quiet = 0
    sample = sampleSize
    #iterate over all of the bands
    for j in xrange(1, dataset.RasterCount + 1): #1 based counter

        band = dataset.GetRasterBand(j)
        NDV = band.GetNoDataValue()

        print "Processing band: " + str(j)       

        #define the interval at which blocks are created
        intervalY = int(band.YSize/1)    
        intervalX = int(band.XSize/2000) #to be changed to sampleSize when working

        #iterate through the rows
        scanBlockCounter = 0

        for i in xrange(0,band.YSize,intervalY):

            #If the next i is going to fail due to the edge of the image/array
            if i + (intervalY*2) < band.YSize:
                numberRows = intervalY
            else:
                numberRows = band.YSize - i

            for h in xrange(0,band.XSize, intervalX):

                if h + (intervalX*2) < band.XSize:
                    numberColumns = intervalX
                else:
                    numberColumns = band.XSize - h

                scanBlock = band.ReadAsArray(h,i,numberColumns, numberRows).astype(numpy.float)

                standardDeviation = numpy.std(scanBlock)
                mean = numpy.mean(scanBlock)

                newMin = mean - (standardDeviation * n)
                newMax = mean + (standardDeviation * n)

                outputBlock = ((scanBlock - newMin)/(newMax-newMin))*255
                outRaster = outdataset.GetRasterBand(j).WriteArray(outputBlock,h,i)#array, xOffset, yOffset


                scanBlockCounter = scanBlockCounter + 1
                #print str(scanBlockCounter) + ": " + str(scanBlock.shape) + str(h)+ ", " + str(intervalX)
                if numberColumns == band.XSize - h:
                    break

                #update progress line
                if not quiet:
                    gdal.TermProgress_nocb( (float(h+1) / band.YSize) )

ここに更新があります:プロファイルモジュールを使用せずに、コードの小さなセクションを関数にラップしたくなかったので、print ステートメントと exit ステートメントを組み合わせて使用​​して、どの行が最も時間がかかっているかについて大まかなアイデアを得ました。幸いなことに（そして、私がどれほど幸運だったかを理解しています）、1行ですべてが引き下げられました。

    outRaster = outdataset.GetRasterBand(j).WriteArray(outputBlock,h,i)#array, xOffset, yOffset

出力ファイルを開いて配列を書き出すとき、GDALは非常に非効率的であるようです。これを念頭に置いて、変更した配列「outBlock」を Python リストに追加し、チャンクを書き出すことにしました。変更したセグメントは次のとおりです。

outputBlock が変更されました...

         #Add the array to a list (tuple)
            outputArrayList.append(outputBlock)

            #Check the interval counter and if it is "time" write out the array
            if len(outputArrayList) >= (intervalX * writeSize) or finisher == 1:

                #Convert the tuple to a numpy array.  Here we horizontally stack the tuple of arrays.
                stacked = numpy.hstack(outputArrayList)

                #Write out the array
                outRaster = outdataset.GetRasterBand(j).WriteArray(stacked,xOffset,i)#array, xOffset, yOffset
                xOffset = xOffset + (intervalX*(intervalX * writeSize))

                #Cleanup to conserve memory
                outputArrayList = list()
                stacked = None
                finisher=0

フィニッシャーは、エッジを処理する単純なフラグです。リストから配列を作成する方法を理解するのに少し時間がかかりました。その中で、numpy.array を使用すると 3 次元配列が作成され (理由を説明したい人はいますか?)、write array には 2 次元配列が必要です。合計処理時間は、2 分弱から 5 分までさまざまです。時間の範囲が存在する理由は何ですか?

投稿してくださった皆様、本当にありがとうございました！次のステップは、実際に Numpy に取り掛かり、追加の最適化のためのベクトル化について学習することです。

4

3 に答える 3

7

numpyデータの操作を高速化する 1 つの方法は、 を使用することvectorizeです。基本的に、ベクトル化は関数を受け取り、配列にマップfする新しい関数を作成します。次に、次のように呼び出されます。gfagg(a)

>>> sqrt_vec = numpy.vectorize(lambda x: x ** 0.5)
>>> sqrt_vec(numpy.arange(10))
array([ 0.        ,  1.        ,  1.41421356,  1.73205081,  2.        ,
        2.23606798,  2.44948974,  2.64575131,  2.82842712,  3.        ])

使用しているデータが利用可能でなければ、これが役立つかどうかははっきりとは言えませんが、上記をvectorized. おそらくこの場合、インデックスの配列を にベクトル化できますReadAsArray(h,i,numberColumns, numberRows)。潜在的なメリットの例を次に示します。

>>> print setup1
import numpy
sqrt_vec = numpy.vectorize(lambda x: x ** 0.5)
>>> print setup2
import numpy
def sqrt_vec(a):
    r = numpy.zeros(len(a))
    for i in xrange(len(a)):
        r[i] = a[i] ** 0.5
    return r
>>> timeit.timeit(stmt='a = sqrt_vec(numpy.arange(1000000))', setup=setup1, number=1)
0.30318188667297363
>>> timeit.timeit(stmt='a = sqrt_vec(numpy.arange(1000000))', setup=setup2, number=1)
4.5400981903076172

15倍のスピードアップ！ndarraynumpy slicing はs のエッジをエレガントに処理することにも注意してください。

>>> a = numpy.arange(25).reshape((5, 5))
>>> a[3:7, 3:7]
array([[18, 19],
       [23, 24]])

ReadAsArrayしたがって、データを に取り込むことができればndarray、エッジ チェックの悪ふざけを行う必要はありません。

再形成に関するあなたの質問について - 再形成は基本的にデータをまったく変更しません。numpyデータにインデックスを付ける「ストライド」を変更するだけです。メソッドを呼び出すと、reshape返される値はデータの新しいビューです。データはまったくコピーまたは変更されず、古いストライド情報を含む古いビューもありません。

>>> a = numpy.arange(25)
>>> b = a.reshape((5, 5))
>>> a
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24])
>>> b
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19],
       [20, 21, 22, 23, 24]])
>>> a[5]
5
>>> b[1][0]
5
>>> a[5] = 4792
>>> b[1][0]
4792
>>> a.strides
(8,)
>>> b.strides
(40, 8)
于 2011-07-11T19:44:44.620 に答える
5

要求通りに答えた。

IO バウンドの場合は、読み取り/書き込みをチャンクする必要があります。約 500 MB のデータを ndarray にダンプし、すべてを処理して書き出し、次の約 500 MB を取得してみてください。必ず ndarray を再利用してください。

于 2011-07-13T19:24:28.600 に答える
2

あなたが何をしているのかを完全に理解しようとせずに、numpy スライスや配列ブロードキャストを使用していないことに気付きました。どちらもコードを高速化するか、少なくとも読みやすくする可能性があります。これらがあなたの問題と密接に関係していない場合は、申し訳ありません。

于 2011-07-11T19:06:12.780 に答える