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

Question

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の平均値としています。

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

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

Comments

[WathMorks]

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

[退会済みユーザー]

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

[Q71]

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

Answer 1

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

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

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

Comments

[退会済みユーザー]

ご返信ありがとうございます。
画像単位で処理するとはどういうことなのでしょうか？よろしくお願いします

[退会済みユーザー]

画像単位で処理するというのは、ある画像に対して行うべき処理を全て行ってから次の画像の処理に移ることを指しているつもりです

質問者さんのソースでは画像座標を決めてから全フレームを処理しているので同じ画像を少なくとも画像の画素数回数分参照することになってしまいます。

この手の画像処理の速度が遅くなる原因は大抵の場合メモリアクセス効率であり、画像座標固定してフレーム間走査をし、次の座標に移るという方式は非常にアクセス効率が悪いため低速になります。
それよりは、フレーム1の全画素に対して処理を行った後次のフレームに移る方が断然高速になります。

RGB平均は単一画像で完結する処理なので画像単位で処理できますし、相関係数も工夫をすれば計算途中のデータを保持することで画像単位で処理することが可能です。
（分散、共分散を計算するために和、二乗和、積和を計算するというのがそれに該当します）

[退会済みユーザー]

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

[退会済みユーザー]

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

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

[退会済みユーザー]

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

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

Answer 2

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()

Comments

[退会済みユーザー]

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

[tiitoi]

追記しました。

Answer 3

画像で試してみました。（全て各ピクセルの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()

Comments

[meg_]

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