numpyの連続した数値列のカウント

2018/04/12 13:28

そうです！写真から最終的にはグラフの作成を考えています

umyu

2018/04/12 13:33

＞itakovichさんへ　質問文にその旨を追加してくださいな。

https://www.geeksforgeeks.org/maximum-consecutive-ones-or-zeros-in-a-binary-array/

2018/04/12 14:33

分かりにくくて、すみません…　追加しておきます

行動規範の内容に同意します

回答8件

普通に考えて全部のピクセルを見ないといけないのでO(N)になります。(worstケース)

最も簡単な方法は並列化です。
今回の場合、単純並列なので、立てることのできるスレッドの数分の1に計算時間を減らすことができます。

あとはアルゴリズムの話になりますが、例えばk番目まで見て、今のマスは0でこれまでの最長の1たちは100個でした。

そうすれば、とりあえず次は、k+100を見てみることにします。
こいつが1であればk+50を見てみることにします。
そして…

なんとなくお察しいただけるかと思いますが、このように探索することで枝切りの可能性が増えます。
例えばk+50の時点で0に出逢えば、k+1~~k+49、k+51~~k+99はもう探索する必要はありません。

このようにすることで常に153620481000個のピクセルを見なければいけなかったところを、より少ないピクセルを探索すれば良いことになります。

後はデータの性質に依存します。
01010101010101010101010101010のように並んでいるのであれば判定が必要な分遅くなります。
11111111110000011100000000011のようであれば、大幅に探索を減らすことができます。

ブロック分けしてdivide-and-conquerでもよいですが。

↓    ↓↓　　↓
0000001 0111011 ...
0000010 1111011 ...
1000010 0000000 ...
1100001 1000001 ...
1000010 0111110 ...
...

比較。（更新：groupby, numpy.where, backward, forward）

python
1from itertools import groupby
2import numpy as np
3import random
4
5from contextlib import contextmanager
6import time
7@contextmanager
8def timer(name):
9    t0 = time.time()
10    yield
11    print(f'[{name}] done in {time.time() - t0:.5f} s')
12
13slen = 1300
14tslen = 1500
15N = 10000
16
17ss = []
18for _i in range(N):
19    random.seed(_i)
20    s_ = np.array([1]*(tslen-slen) + [random.randint(0, 1) for _ in range(slen)])
21    ss.append(s_)
22
23with timer('groupby'):
24    mm0 = []
25    for _i in range(N):
26        s = ss[_i]
27        ans = [len(list(g)) for k, g in groupby(s)]
28        m = max(ans[::2])
29        mm0.append(m)
30        
31with timer('numpy.where'):
32    mm1 = []
33    for _i in range(N):
34        a = ss[_i]
35        offset = 1
36        diff_a = a[offset:] - a[:-offset]
37
38        pos1to0 = np.where(diff_a ==-1)[0] + offset
39        pos0to1 = np.where(diff_a ==+1)[0] + offset
40
41        pos0to1 = np.append([0], pos0to1)
42        m = (pos1to0 - pos0to1[:len(pos1to0)]).max()
43        mm1.append(m)
44
45with timer('backward-search'):
46    def forward_search(s):
47        j = 0
48        while j < len(s) and s[j] == 1:
49            j += 1
50        return j
51    
52    def backward_search(s, m):
53        if m >= len(s):
54            return False, len(s), 0
55        t = s[m::-1]
56        c = forward_search(t)
57        return True, m, c
58
59    mm2 = []
60    for _i in range(N):
61        s = ss[_i]
62        m = 0
63        i = forward_search(s)
64        m = i
65        while i < len(s):
66            fforward, j, c = backward_search(s[i:], m)
67            i += j
68            if fforward:
69                j = forward_search(s[i:])
70                c += j
71                i += j
72                if c > m:
73                    m = c
74        mm2.append(m)
75        
76with timer('forward-search'):
77    mm3 = []
78    for _i in range(N):
79        s = ss[_i]
80        m = 0
81        c = 0
82        for v in s:
83            if v == 1:
84                c += 1
85            else:
86                if c > m:
87                    m = c
88                c = 0
89        if c > m:
90            m = c
91        mm3.append(m)
92
93assert np.all(mm0 == mm1)
94assert np.all(mm0 == mm2)
95assert np.all(mm0 == mm3)

雑談

高速化を行う際には、

アルゴリズムの改良
コーディングの改良

に分けて考えるとスッキリします。

今回の場合、最低でもO(N)で各ピクセルを見ないと答えが出ないような問題です。
すると、探索の打ち切りの可能性を検討した後、コーディングの改良を検討します。

（目的とあっていない可能性がありますが、
最長の1の塊を探索するのであれば、最長になりえない領域をできるだけスキップする作戦を考えます。）

次に書き方を工夫します。
（ブロック化してSIMD命令を利用するのが、C++などのコンパイル言語を使う際に検討する必要があります。）
Pythonの場合、numpy、opencvのライブラリを出来る限り利用できるようにするのが近道です。
すると、numpyのメソッド、opencvのメソッド、pythonのbuilt-inメソッドを使用することを検討すべきです。

一個ずつピクセルを見ていくより、ベクトルとみなしてsumが速いことは予想できます。
ならば、sumだけで済むように前処理を検討するのが良さそうです。

opencvのオープニングを使うことで先にノイズを消してから、sumを行うのが一番速そうです。

他にもバイナリとして扱ってビット演算を用いることで高速化することも考えられます。
Pythonでそれをやるのは少しちぐはぐ感が残ってしまいます。
C++などで書き換えるだけで、もとのアルゴリズムで高速に処理できることが見込まれるためです。

（重複するようなケースが出現する場合、キャッシュをうまく使うことも重要です。
https://pypi.python.org/pypi/fastcache/0.4.3
メモリと演算のバランスをチューニングするがよいコード書くために大切です。）

高速化することも大事ですが、楽して高速化することも大事です。

投稿2018/04/12 13:40

編集2018/04/15 13:41

総合スコア8562

2018/04/12 14:43

分かりやすい説明ありがとうございます！さっそく実装してみます

2018/04/12 19:21

やってみました！意外と難しくて、なかなかうまくいきませんデータの性質では 111111111111111111111110000000000000000000 111111111111111111111111100000100000000000 111111111111111000000000001000000000000000 のように一列目はすべて１から始まっていますノイズをカウントしたくないので連続しているという１の数を求めたいという標記にしましたよろしければ探索を減らすソースコードを乗せていただけないでしょうか

2018/04/15 05:27

今気づいたのですが、連続したブロックのうち最長のものがほしいのか、 1個以上連続で出現している1の総数がほしいのか、どちらでしょう？

2018/04/15 05:43

前者ならデータの性質の応じて最後に追記したコードで処理できます。連続する1が少ない場合、素直に順番に探索するのが最速です。連続する1が多い場合、枝切りをすることで速めることができます。ノイズをカウントせずに、1個以上連続で出現している1の総数がほしいのであれば、 http://labs.eecs.tottori-u.ac.jp/sd/Member/oyamada/OpenCV/html/py_tutorials/py_imgproc/py_morphological_ops/py_morphological_ops.html オープニングを行うことでノイズを消せます。後に、和を取ることによって1の数を得ることができます。

2018/04/16 07:11

ありがとうございます！！こっちもやってみます！！

行動規範の内容に同意します

ベストアンサー

0と1の境界線が欲しいなら、変化点を探して距離を測ればいいんじゃないかなぁ。変化点なんで、まず微分します。

python
1a = np.array([1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0])
2offset = 1
3diff_a = a[offset:] - a[:-offset]

1->0は-1の変化、0->1は+1の変化なので、それぞれの場所をnp.argwhereを使って抽出します。

pos1to0 = np.where(diff_a ==-1)[0] + offset # -> [8, 10]
pos0to1 = np.where(diff_a ==+1)[0] + offset # -> [9]

右端が必ず0で終わると仮定すれば、pos0to1を微調整してインデックスの差分を取れば1が連続する長さを測れます。

pos0to1 = np.append([0], pos0to1)
len11 = pos1to0 - pos0to1

pythonのfor-loopが遅いので行列のインデックスを駆使した方法です。理解してコーディングしないと収拾がつかなくなるので、適度な用法用量を守ってご参考にしてください。

追記

mkgreyさんのコードそのまま回したケース
[groupby] done in 0.39933 s
[forward-search] done in 0.12787 s
[backward-search] done in 0.05486 s

最大の長さを返すよう若干修正したバージョン


In [51]: slen = 1000         
    ...: tslen = 1200           
    ...: N = 1000      
    ...:                            
    ...: random.seed(0)                                  
    ...: s_org = [1]*(tslen-slen) + [random.randint(0, 1) for _ in range(1000)]
    ...:                                      
In [52]: with timer("numpy"):
    ...:     a = np.array(s_org)
    ...:     offset = 1
    ...:     diff_a = a[1:] - a[:-1]
    ...:     pos1to0 = np.where(diff_a ==-1)[0] + offset 
    ...:     pos0to1 = np.where(diff_a ==+1)[0] + offset
    ...:     pos0to1 = np.append([0], pos0to1)
    ...:     len11 = (pos1to0 - pos0to1).max()
    ...:         
[numpy] done in 0.00026 s
In [53]: assert len11 == 202

投稿2018/04/15 08:37

編集2018/04/15 12:17

総合スコア3601

2018/04/15 09:21

文句をつけるようで申し訳ないのですが、手元の環境でやってみたところfor文を回したほうが速かったです。私の実装に自信がないので、ご確認していただいてよろしいでしょうか。

2018/04/15 09:58

なるへそ。コメントありがとうございます。しばしお待ちを。

2018/04/15 12:18

mkgreiさん、比較してみました。同じものを比べてますかね？？np.arrayを前もって変換しておけば実際の計算時間はこれの半分程度で済みました。

2018/04/15 12:32

ご確認していただいてありがとうございます。 a = np.array([1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0]) で比較していました。 numpyの方が長さが増えても計算時間がほとんど変わらず、速いですね。

2018/04/15 12:57

numpyのテクニックですね(ちなみにMATLABも同様)。アルゴリズムを工夫しても実装であしを引っ張られてしまうのは辛いところ。

2018/04/15 13:15

whereを使えば、最長だけでなくすべての区間も取り出せるのがいいですね。 lookupなら素直に組んでもワンチャンあるかとも思いましたが、ダメですね。

2018/04/16 06:09

返信遅れてすいません分かりやすい解説ありがとうございます！！

2018/04/19 09:23

np.whereで-1を探しても出てこないのですが、なぜでしょうか

2018/04/19 09:35

値が変化してない・・・ということは無いですか？

2018/04/20 02:58

確認してみましたが値は変化していました…

2018/04/20 07:01

from PIL import Image import numpy as np import matplotlib.pyplot as plt np.set_printoptions(threshold=np.inf) #閾値 t=100 for i in range(0,1): filename='C://toisi//f_%05d.jpg'%(i) im = np.array(Image.open(filename).convert('L')) #二値化 im[im<t]=0 im[im>=t]=1 #print(im.shape, im.dtype) #print(im[1,1250]) #0の数をカウント for i in range(1535): a = im[i] offset = 1 diff_a = a[offset:] - a[:-offset] pos0to1 = np.where(diff_a ==+1)[0] + offset pos1to0 = np.where(diff_a ==-1)[0] + offset #print(pos0to1) #print(pos1to0) print(diff_a)

2018/04/20 07:06

diff_aをプリントしてみてみると [0,0,0,0,0,0,0,0,0,0,0,0,255,1,1,0,0,0,0,0,0,0] のような感じになってしまい、出るはずのない255という値がでてきてしまいました

LouiS0616

2018/04/20 07:22

dtypeがnp.uint8なのでしょう。アンダーフローが発生しています。この場合気にせず diff_a == 255 の箇所を利用すれば良いかと思います。

2018/04/20 07:31

フォローありがとうございます。画像を読み込んだ結果がuint8になっているためで、その型uint8同士の演算もuint8になって、結果アンダーフローして255に飛んだわけですか。im=im.astype(np.int)でint型にするか、LouiS0616さんのコメントにあるとおり255を利用するかですね。迂闊でした。

2018/04/20 09:03

ありがとうございます！！できました！少し聞きたいのですがpos0to1のときは[9]でpos1to0の中身が[]のように何もないときpos1to0-pos0to1でエラーが表示されてしまうのですがif文で無理やり入れてしまうなどの対処法でよかったでしょうか度々すいません。よろしくお願いいたします

2018/04/20 09:40

m = (pos1to0 - pos0to1[:len(pos1to0)]).max() どっちが短いか見極めて、スライスしてください。

2018/04/25 11:46

ありがとうございます！！

行動規範の内容に同意します

一列目はすべて１から始まっています

本来の目的から考えると、私なら左端を塗りつぶした画像を新たに作成します。
OpenCVを使えば内部はネイティブ処理してくれると思うので速度も期待できます。
以下例では左端の黒を白ですが、逆も同様の考えでできます。

Python
1import cv2
2img = cv2.imread( 'bin.bmp', cv2.IMREAD_GRAYSCALE) # 簡単のため2値ビットマップを読込
3print(img.shape)
4# y方向に左端が黒を白に塗りつぶし
5for r in range(img.shape[0]):
6    if img[r,0] == 0:
7        print('fill',r)
8        cv2.floodFill(img,None,(0,r),255)
9cv2.imwrite('out.jpg',img)

bin.bmp

out.jpg

投稿2018/04/13 03:48

8524ba23

総合スコア38352

初めて1が来たところから0が来るまでfor文やwhile文でカウントし、カウントした数を別のリストに格納する。if文で評価し、0だったらbreak文で切る

これでどうでしょうか?

投稿2018/04/12 13:13

退会済みユーザー

総合スコア0

2018/04/12 13:17

早速の回答ありがとうございます！一度それでやったのですが時間がかかってしまいました 1536×2048の画像で，1000個のファイル処理をする予定です

退会済みユーザー

2018/04/12 13:24

なるほど... もっと効率的な方法を考えておきます。

2018/04/12 13:29

すみませんありがとうございます！！

退会済みユーザー

2018/04/12 13:48 編集

今思いついたアルゴリズムを書きたいと思います。同じ数字が続くということはその二つの整数の差は"0"です。0～1または1～0が変わるということは差の絶対値は1です。この差の絶対値が1というところが変わる境目なのでこれを数えるのはどうでしょうか? 差が1のところと-1のところを別々にカウントするという意味です。 (データが多いのでこれも時間がかかるかもしれません)

2018/04/12 14:37

ありがとうございます！！やってみます！