質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

89.62%

計算時間を圧倒的に短くする方法が知りたいです。

受付中

回答 3

投稿

  • 評価
  • クリップ 2
  • VIEW 295

santamarianeed

score 15

30秒の動画において、、ある1ピクセルを指定し、他の全てのピクセルの強度の相関係数を算出し、ヒートマップで可視化する作業を行ってます。ピクセルの個数が300×300=90000なので計算時間がものすごいかかってしまいます。

# -*- coding: utf-8 -*-
import cv2
import matplotlib.pyplot as plt
import numpy as np

cap = cv2.VideoCapture("test.mkv")
count = cap.get(cv2.CAP_PROP_FRAME_COUNT)
fps = cap.get(cv2.CAP_PROP_FPS)

# 先に全フレームを読み込む。
frames = []
while True:
    ret, frame = cap.read()
    if not ret:
        break


    frames.append(frame)

frames = np.array(frames)
numE=count
print(numE)
########################################################################
# 場所を決める
#点滅の場所
xpx=int(raw_input("点滅1 X:"))
ypx=int(raw_input("点滅1 Y:"))
########################################################################
# フレーム [numS,numE] の範囲で各フレームの [ymin, ymax]x[xmin, xmax] の画素の平均を計算する。
frame_no = np.arange(count)
print("frame number",count)
########## read the target pixel intensity
intGt=np.zeros((0)) # the target intensity
print("read target position x=%d y=%d intensity"%(xpx,ypx))
for i in frame_no:
    # フレーム frame_no を取得する。
    cap.set(cv2.CAP_PROP_POS_FRAMES, i)
    ret, frame = cap.read()
    if not ret:
        print('Failed to grab frame.')
        break
    # 平均画素値を計算する。
    G=np.average(frame[xpx,ypx])
    intGt=np.append(intGt,G)

print("read the target intensity!")

# read first frame (frame zero)
# フレーム frame_no を取得する。
cap.set(cv2.CAP_PROP_POS_FRAMES, 0)
ret, frame = cap.read()
XS=np.shape(frame)[0]
YS=np.shape(frame)[1]
# show the target pixel intesity
fig, ax = plt.subplots()
ax.invert_yaxis()
#im=ax.imshow(std,cmap="gray")
ax.set_title("(water)")
plt.plot(intGt,".")
plt.ylim([0,255])
#plt.imshow(std, cmap="jet",vmin=0,vmax=32)
#plt.imshow(std, cmap="jet")
#plt.colorbar()
plt.show()



raw_input("stop")
import time as t
corrMap=np.zeros((0))

x=y=0
for y in range(YS):
    for x in range(XS):
        t0=t.time() # start clock
        intG = np.zeros((0)) # reset intensity
        for i in frame_no:
            # フレーム frame_no を取得する。
            cap.set(cv2.CAP_PROP_POS_FRAMES, i)
            ret, frame = cap.read()
            if not ret:
                print('Failed to grab frame.')
                break
            # 平均画素値を計算する。
            G=np.average(frame[x,y])
            intG=np.append(intG,G)
            # correlation between target and current pixel
        CC=np.corrcoef(intGt,intG) # correlation coefficient
        print(np.shape(CC),CC[0,1])
        print("now at position x=%d y=%d, cc=%4.2g"%(x,y,CC[0,1]))
        corrMap=np.append(corrMap,CC[0,1])
        print("it took %4.2g secs"%(t.time()-t0))
        print("you need %4.2g mins"%((t.time()-t0)*XS*YS/60))

###########
corrMap=corrMap.reshape(shape(frame))
raw_input("stop")
# ヒートマップで描画する。
fig, ax = plt.subplots()
ax.invert_yaxis()
im=ax.imshow(corrMap,cmap="gray")
ax.set_title("(water)")
#plt.plot(intG,".")
#plt.imshow(std, cmap="jet",vmin=0,vmax=32)
#plt.imshow(std, cmap="jet")
#plt.colorbar()
plt.show()


上記のようなプログラムを書いて、シミュレーション時間を計算してみましたが多分20日くらい?かかる感じです。強度はRGBの平均値としています。

何か解決策はあるのでしょうか?

ちなみに求めたいヒートマップのイメージは添付の通りです。よろしくお願いします

イメージ説明

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

質問への追記・修正、ベストアンサー選択の依頼

  • GradMaisonTokyo

    2020/01/16 20:15

    「ピクセルの強度の相関係数」の定義は何ですか?

    キャンセル

  • santamarianeed

    2020/01/17 14:25

    強度は1フレームおきにplotされます。1フレームは30fpsなので0.03秒おきにplotされます。私はこの0.03秒での相関関係を見ています

    キャンセル

  • Q71

    2020/01/23 18:52

    メモリには、フレーム,高さ,幅,色 の順で並びます。xを全部走査してからy、yの走査が済んでからフレームを変更するようにすると、メモリへのアクセスが良くなります。
    ピクセル間の強度とは、1フレームの中で計算するのでしょうか。それともフレーム間で計算するのでしょうか。

    キャンセル

回答 3

+1

何やってるんだろうと思ったら画素単位のループの内側で毎フレームRGB平均とってフレーム間データ取り出して正規化相互相関とってるんですね。そら重いはずですわ。

pythonはあまり詳しくないので画像処理のやり方としての回答になりますが
frame間走査を最も内側にするとメモリアクセスが滅茶苦茶とびとびになるので、各画素でframe間走査するのはやめて画像単位で処理しましょう。
まずframe間操作は一番外側にして各フレームで画像単位でRGB平均を計算し、フレーム間の平均・分散・共分散を計算して正規化相互相関計算を行います。
フレーム間の和・二乗和・相関計算対象との積の和を用意してやれば1回のframe間走査で処理できると思います。

何か勘違いがあったらすいませんが、こんな感じでどうでしょうか。

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2020/01/23 17:50

    すなわち、画像座標の指定(写真の指定した1ピクセル)は一番最後の段階でやった方がベターってことで大丈夫ですか?

    キャンセル

  • 2020/01/23 18:18

    最後=内側のループ という認識でしたら大丈夫です。

    厳密ではありませんが簡単にその方が良い理由を説明します。
    前述の通り物凄く計算が複雑でない限りは大抵は計算時間よりもメモリアクセス時間が大きくなります。
    メモリアクセスの回数が減らせば良い訳ですが、CPUはメモリからデータを引っ張ってくる際に、対象データのあるメモリ位置に連続するデータを一緒にまとめて引っ張って来ます。
    画像の場合連続するデータは画像X方向に隣接する画素になります(画像右端の場合は次の行に連続)
    つまり画像X方向のループを一番内側にすれば一緒に引っ張って来たデータが使えるのでメモリアクセスの回数が減るということになります。
    同じ理屈でXのループの外側はYのループ、その外側はフレームのループにするのがループの順序としては最適です。

    キャンセル

  • 2020/01/23 18:27 編集

    多分ですけど他の方の回答に示されてるような画像同士での演算処理をうまく組み合わせればX,Yのループはわざわざ書かなくてもいけるかも知れませんね(もちろん内部的には行われてる)

    tiitoiさんが別の切り口の速度問題に言及されてるのでので用が足りたらそれだけで十分かもです(pythonの細かい仕様まで把握してなくてすいません)

    キャンセル

0

画像で試してみました。(全て各ピクセルのRGBの平均としています)

import cv2
import matplotlib.pyplot as plt

img = cv2.imread("Flower1.jpg")

shape = img.shape

data = (img[:,:,0] + img[:,:,1] + img[:,:,2]) / 3

fig, ax = plt.subplots()
plt.imshow(data, cmap=plt.cm.Blues)
plt.colorbar()

plt.show()


イメージ説明
イメージ説明

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2020/01/16 20:39

    上記画像は求めているものとは違うかもしれませんが、numpyの計算はforループ使わない方法でやった方が良いです。

    キャンセル

0

300x300 の30秒ぐらいの動画であれば、多分、全部メモリに乗り切るので、まず全フレームを読み込んで、(フレーム数, 画像の高さ, 画像の幅, 3) の numpy 配列を作りましょう。
numpy 配列を作ったら、mean_imgs = np.mean(imgs, axis=-1) でチャンネルごとの平均をとります。
相関係数計算のところは、numpy.corrcoef() の引数に1次元配列を渡さないといけないので、仕方なく for 文で各画素ごとに計算しています。
下記の環境、コードでだいたい1分ぐらいで計算できました。

質問のコードが遅い原因

  • np.append() は、追加するたびに配列を作り直すので、リストと同じように使うと、とても遅くなります。基本的にリストに append() して、最後に numpy.array() で numpy 配列に変換したほうが早いです。

  • cap.set(cv2.CAP_PROP_POS_FRAMES, i) は遅いので、各フレーム順番にとってく場合は、cap.set(cv2.CAP_PROP_POS_FRAMES, 0) で1回だけ位置を最初のフレームに戻せばよいでしょう。

試した環境

  • 576x768 の動画 795 フレーム分
  • Ubuntu 16.04
  • Intel(R) Core(TM) i7-6800K CPU @ 3.40GHz
  • 66.96 秒

コード

import time
from urllib.request import urlretrieve

import cv2
import numpy as np

start = time.time()

# サンプル用動画をダウンロードする。
video_path = "sample.avi"
video_url = "https://github.com/opencv/opencv/raw/master/samples/data/vtest.avi"
urlretrieve(video_url, video_path)

# 動画の全フレームを読み込む。
cap = cv2.VideoCapture(filepath)
# 1フレームずつ取得する。
imgs = []
while True:
    ret, frame = cap.read()
    if not ret:
        break
    imgs.append(frame)

# チャンネルごとに平均をとる。
mean_imgs = np.mean(imgs, axis=-1)
print(mean_imgs.shape)  # (795, 576, 768)

# ピクセルごとに相関係数を計算する。
x0, y0 = 10, 10  # 対称とする点

h, w = mean_imgs.shape[1:3]  # 画像の幅、高さ
vec1 = mean_imgs[:, y0, x0]

corr_img = np.empty((h, w))
for y in range(h):
    for x in range(w):
        vec2 = mean_imgs[:, y, x]
        # 相関行列を計算する。
        C = np.corrcoef(vec1, vec2)
        corr_img[y, x] = C[0, 1]  # x, y の相関係数

print(f"Done! {time.time() - start}s")  # 66.96258449554443s

fig, ax = plt.subplots()
ax.imshow(corr_img)

plt.show()

イメージ説明

追記

1フレームの強度の振れ幅がピクセルごとに違うので、正規化を行いたいと思います。

チャンネルごとに平均をとって、(フレーム数, 高さ, 幅) の配列を作ったあと、axis=0 方向、つまり、画素ごとに [-1, 1] の範囲に正規化を行えばよいです。

import time
from urllib.request import urlretrieve

import cv2
import matplotlib.pyplot as plt
import numpy as np

start = time.time()

# サンプル用動画をダウンロードする。
video_path = "sample.avi"
video_url = "https://github.com/opencv/opencv/raw/master/samples/data/vtest.avi"
urlretrieve(video_url, video_path)

# 動画の全フレームを読み込む。
cap = cv2.VideoCapture(video_path)
# 1フレームずつ取得する。
imgs = []
while True:
    ret, frame = cap.read()
    if not ret:
        break
    imgs.append(frame)

# チャンネルごとに平均をとる。
mean = np.mean(imgs, axis=-1)
del imgs  # メモリ節約のため、不要になったので、削除する。
print(mean.shape)  # (795, 576, 768)

# 範囲を [-1, 1] にする。
mins, maxs = mean.min(axis=0), mean.max(axis=0)

# メモリ節約のため、演算はインラインで行う。
mean -= mins
mean *= 2 / (maxs - mins)
mean -= 1

# ピクセルごとに相関係数を計算する。
x0, y0 = 10, 10  # 対称とする点

h, w = mean.shape[1:3]  # 画像の幅、高さ
vec1 = mean[:, y0, x0]

corr_img = np.empty((h, w))
for y in range(h):
    for x in range(w):
        vec2 = mean[:, y, x]
        # 相関行列を計算する。
        C = np.corrcoef(vec1, vec2)
        corr_img[y, x] = C[0, 1]  # x, y の相関係数

print(f"Done! {time.time() - start}s")  # 66.96258449554443s

fig, ax = plt.subplots()
ax.imshow(corr_img)

plt.show()

イメージ説明

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • この投稿は削除されました

  • 2020/02/03 11:58

    1フレームの強度の振れ幅がピクセルごとに違うので、正規化を行いたいと思います。
    全ピクセルに2×((y(t)-ymin/ymax-ymin))-1を行って、この計算によって得られた時系列強度を使って、同様corrlation mapに表示させるにはどうしたら良いでしょうか?

    キャンセル

  • 2020/02/04 00:24

    追記しました。

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 89.62%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる