pythonを用いてフラクタル次元解析[ボックスカウンティング法]の高速化を実現したい。

Question

前提・実現したいこと

ボックスカウンティング法のプログラムを自作で作成したのですが、あまりにも処理が遅いため改善したいと考えております。
プログラムは二次元のグラフに対しての処理で、100データを50データずらしながらボックスカウンティング法を用いることを検討しています。現状総データ3000に対して3388秒かかってしまいます。

該当のソースコード

df:データ[x値、y値]
x: xの値における最大値
y: yの絶対値における最大値
b: ボックスの大きさ

def count(df, x, y, b):
    i = 0
    j = 0
    c = 0
    x = math.ceil(x/b) * b
    y = math.ceil(y/b) * b
    flag1 = False
    flag2 = False
    while i <= x and j <= y:
        flag1 = False
        flag2 = False
        if i + b <= x and j + b <= y:
            for m, n in df:
                if (i <= abs(m) < i+b) and (j <= abs(n) < j+b):
                    if (np.sign(n) == 1) and (flag1 == 0):
                        c += 1
                        flag1 = True
                    elif  (np.sign(n) == -1) and (flag2 == 0):
                        c += 1
                        flag2 = True
                if (flag1 == 1) and (flag2 == 1):
                    break
        i += b
        if i >= x:
            i = 0
            j += b
    return c

試したこと

higich法やマルチフラクタル解析について調べてみたのですがよく理解できませんでした。

Comments

[tiitoi]

こういうループが沢山あるコードは Python はかなり遅いです。
Python で利用するにしても、このアルゴリズム部分だけ CかC++で書き直して、Python から呼び出すようにしたらどうでしょうか。
ロジックがまったく同じコードでも C++ で書き直すだけで300倍ぐらい早くなります。

[A_kirisaki]

いまなら C++ よりちょっと遅いけどメモリ安全な Rust もアルヨ！

[grintea]

ありがとうございます！
ほかの言語については詳しくないのですがやってみます！

[ppaul]

フラクタル次元を推定したいのなら、bを6個ぐらいとって収束の様子を調べればいいですね。
pythonでやるにしても、C/C++でやるにしても、局所的にコードをいじるのではなく全体のアルゴリズムを見直せばo(nlog(n))にできるので、まずそれを考えてみたほうが良いと思います。

Answer 1

grinteaさんのコードと端の余りの処理がちょっと違うかもしれませんが、cvを使って高速化してみました。

import cv2

def boxcount(image, b):
    x, y = image.shape
    xrepeat = x // b
    yrepeat = y // b
    x = xrepeat * b
    y = yrepeat * b
    _, img = cv2.threshold(cv2.bitwise_not(image[:x, :y]), 1, 255, cv2.THRESH_BINARY)
    _, img_b = cv2.threshold(cv2.resize(img, (xrepeat,yrepeat)), 1, 255, cv2.THRESH_BINARY)
    return img_b.sum()//255

if __name__ == '__main__':
    import sys
    import time
    imagefile = sys.argv[1]
    image = cv2.imread(imagefile, cv2.IMREAD_GRAYSCALE)
    print(f'ImageFile: {imagefile}  Size: ({image.shape[0]}, {image.shape[1]})')
    start = time.time()
    print('BoxCount:', boxcount(image, 2))
    elapsed_time = time.time() - start
    print (f"elapsed_time:{elapsed_time}sec")

実行結果は以下です。

> python boxcount.py koch.png
ImageFile: koch.png  Size: (16384, 16384)
BoxCount: 64876
elapsed_time:8.804325342178345sec

16384 × 16384 の画像を処理して9秒弱でした。
cvの関数はたぶんC/C++で書いてあるので、それほど遅くはないと思います。

Comments

[grintea]

たくさんのご意見ありがとうございます！
是非参考にしていただきます！