Edited at

画像処理で2重ループをなるべく書くな

More than 1 year has passed since last update.

注意:検証が不十分な部分があり、そのままでは動作しない記述です。

 また、プログラムの例は、見本的なプログラムとは程遠いものであることをお許しください。

画像処理では旧来より2次元配列へのアクセスという発想で2重ループで書かれる処理の内容が多い。

OpenCVのC++で使う場合でも2重ループでの画素に対するアクセスの書き方しだいで、処理速度が変わることが知られています。

これがnumpyを使う場合には、2重ループでの画素に対するアクセスすると、そのアクセスのエラーチェックのために耐え難いほど遅くなることが分かっています(ほぼpython言語を使いつつ、ループを使いながら高速化したい場合にはCythonを調べてください。)。

ですから、この文章では「画像処理で2重ループをなるべく書くな」と主張します。


 2重ループ処理のほとんどは、ライブラリの関数として実装されています。

 2重ループで書きたくなる処理のうちほとんどは、既にライブラリの関数として実装されています。これはnumpy.arrayだけの問題ではなく、OpenCVのcv::Mat型についても当てはまります。cv::Mat型に対して用意されている関数やメソッドを調べてみましょう。行いたい作業をキーワードとして、ライブラリ名と組み合わせて検索してみてください。かなりの処理が既に実装されているはずです。2重ループの中でif()を使う必要があるのかどうか疑ってかかってみてください。cv::max(A, B)を使うだけでも、if文や2重ループの出番がなくなります。コードは簡潔になりますし、たいがいの場合、自分でループを書くよりも速い結果になるはずです。


2重ループを書くよりも、OpenCVにある関数を使うべき理由


  • OpenCVで用意されている関数の方が、自前のコードよりも多くの場合速い。

  • OpenCVで用意されている関数の方が、マルチコアやGPUなどへの対応がよい。

  • OpenCVで用意されている関数の方が、OpenBLASなどの 最適化された線形代数のライブラリを効率的に利用できる。

  • OpenCVで用意されている関数の方が、SIMD命令などを対応していることがある。

  • OpenCVで用意されている関数を使うと、2重ループの中でif文を書く必要がなくなる場合が多く、明快さと高速性とを両立させることことがある。

numpyは、arrayの添え字に

a[::2]のように、スライスを上手に使うことや

a[a <= 5] = 0

のように論理式で配列の範囲を指定するやりかたなどを使うことをすると

for文でループを書かなくてはならないことが減るはずです。

numpyで明示的にループを書くと極端に遅くなる

 とは言っても、ループを書かなくてはならないこともあるはずなので、次の記事を紹介しておく。

 cv::Mat::forEachを使った高速なピクセル操作

以下に示す2重ループのコードは、高速化を意識しない記述になっています。

この文章の中では、そのようなコードよりもOpenCVの標準の関数の方が格段に速い結果を与えていることを示すものです。その理由で、「画像処理で2重ループをなるべく書くな」ということを主張します。


悪いコードの見本


.cpp


for(int i = 0; i < img.rows; i++){
for(int j = 0; j < img.cols; j++){
if(img.channels() == 3){
カラー画像に対する処理
}else if(img.channels() == 1){
グレースケール画像に対する処理
}
}
}


.cpp

if(img.channels() == 3){

for(int i = 0; i < img.rows; i++){
for(int j = 0; j < img.cols; j++){
カラー画像に対する処理
}
}
}else if(img.channels() == 1){
for(int i = 0; i < img.rows; i++){
for(int j = 0; j < img.cols; j++){
グレースケール画像に対する処理
}
}
}


cv::maxをつかえばいい事例

2つの配列の中から大きな画素を選ぶ。

cv2.max

max(src1, src2[, dst]) -> dst

numpy.maximum(x1, x2[, out])


ex_max.cpp

#include <opencv2/opencv.hpp>


int main(int agrc, char* argv[]){
std::string name1 = "lena.png";
std::string name2 = "home.jpg";

cv::Mat img1=cv::imread(name1, 0);
cv::Mat img2=cv::imread(name2, 0);

int64 time0 = cv::getTickCount();
cv::Mat img3=cv::max(img1, img2);
int64 time1 = cv::getTickCount();

std::cout << (time1-time0)/cv::getTickFrequency() << std::endl;

cv::Mat img4=cv::Mat::zeros(img1.size(), img1.type());
int64 time2 = cv::getTickCount();
for(int i=0; i<img1.rows; i++){
for(int j=0; j<img1.cols;j++){
img4.at<uchar>(i,j) = std::max(img1.at<uchar>(i,j), img2.at<uchar>(i,j));
}
}
int64 time3 = cv::getTickCount();
std::cout << (time3-time2)/cv::getTickFrequency() << std::endl;

cv::imwrite("img3.png", img3);
cv::imwrite("img4.png", img4);
}


$ ./ex_max

0.000244705

0.00290059

$


cv::magnitudeを使えばいい事例

sqrt()を使う必要がなくなる。

cv2.magnitude

cv::Sobel(img, gx, img.depth(), 1, 0, 3);     

cv::Sobel(img, gy, img.depth(), 0, 1, 3);
cv::Mat mag(gx.size(), gx.type());
cv::magnitude(gx, gy, mag);


ex_magnitude.py

import cv2

import numpy as np
x=np.linspace(0, 180, 11)
si=np.sin(np.pi*x/180.0)
co=np.cos(np.pi*x/180.0)
r=cv2.magnitude(si, co)
print r


差の絶対値を求めることは多いですね。

地道に計算時間を節約。


ex_absdiff.cpp

#include <opencv2/opencv.hpp>


int main(int agrc, char* argv[]){
std::string name1 = "../../data/lena.png";
std::string name2 = "../../data/home.jpg";

cv::Mat img1=cv::imread(name1, 0);
cv::Mat img2=cv::imread(name2, 0);

int64 time0 = cv::getTickCount();
cv::Mat img3;
cv::absdiff(img1, img2, img3);
int64 time1 = cv::getTickCount();

std::cout << (time1-time0)/cv::getTickFrequency() << std::endl;

cv::Mat img4=cv::Mat::zeros(img1.size(), img1.type());
int64 time2 = cv::getTickCount();
for(int i=0; i<img1.rows; i++){
for(int j=0; j<img1.cols;j++){
img4.at<uchar>(i,j) = std::abs(img1.at<uchar>(i,j) - img2.at<uchar>(i,j));
}
}
int64 time3 = cv::getTickCount();
std::cout << (time3-time2)/cv::getTickFrequency() << std::endl;

cv::imwrite("ex_absdiff.png", img3);
cv::imwrite("ex_absdiff_loop.png", img4);
}



ex_absdiff.py

import cv2

import numpy as np
x = np.linspace(0, 180, 1000000)
si = np.sin(np.pi*x/180.0)
co = np.cos(np.pi*x/180.0)
t1 = cv2.getTickCount()
for i in range(100):
r = cv2.absdiff(si, co)
t2 = cv2.getTickCount()
print r
print (t2-t1)/cv2.getTickFrequency()

t3 = cv2.getTickCount()
for i in range(100):
r2 = np.abs(si-co)
t4 = cv2.getTickCount()
print r2
print (t4-t3)/cv2.getTickFrequency()



実行結果.txt

[[ 1.        ]

[ 0.99999686]
[ 0.99999372]
...,
[ 1.00000628]
[ 1.00000314]
[ 1. ]]
0.524811459022
[ 1. 0.99999686 0.99999372 ..., 1.00000628 1.00000314 1. ]
0.915332887212


画像を最大値で正規化したいことは多いですね。

画像を最大値で正規化するために、ループで最大最小を求めてから、画像を正規化してしまう。

まず最大最小を求めるのにループを書く必要がない。

cv::minMaxLoc

cv2.minMaxLoc

を使おう。これらの関数を用いた場合には、マルチコアへの最適化がされている可能性が高い。

さらに、最大値・最小値で正規化する場合には、既に関数が標準で用意されている。

C++: void normalize(InputArray src, OutputArray dst, double alpha=1, double beta=0, int norm_type=NORM_L2, int dtype=-1, InputArray mask=noArray() )


配列を重みつきで加算する。

これはBLASなどのライブラリには

SAXPY y = a*x + y

などがあるのに対応している関数のようです。

dst(I)= scale*src1(I) + src2(I)

cv::scaleAdd

cv2.scaleAdd

次の2つの記述は、ほぼ同等の結果になる。ただし端数の扱いが異なるため一致しない。

scaleAddを用いた方がpythonインタプリタ上の計測では格段に速い。


ex_scaleAdd.cpp

#include <opencv2/opencv.hpp>


int main(int agrc, char* argv[]){
std::string name1 = "../../data/lena.png";
std::string name2 = "../../data/home.jpg";

cv::Mat img1=cv::imread(name1, 0);
cv::Mat img2=cv::imread(name2, 0);

int64 time0 = cv::getTickCount();
cv::Mat img3;
cv::scaleAdd(img1, 0.5, img2, img3);
int64 time1 = cv::getTickCount();

std::cout << (time1-time0)/cv::getTickFrequency() << std::endl;

cv::Mat img4=cv::Mat::zeros(img1.size(), img1.type());
int64 time2 = cv::getTickCount();
for(int i=0; i<img1.rows; i++){
for(int j=0; j<img1.cols;j++){
img4.at<uchar>(i,j) = cv::saturate_cast<uchar>(0.5*img1.at<uchar>(i,j) + img2.at<uchar>(i,j));
}
}
int64 time3 = cv::getTickCount();
std::cout << (time3-time2)/cv::getTickFrequency() << std::endl;

cv::imwrite("ex_scaleAdd.png", img3);
cv::imwrite("ex_scaleAdd_loop.png", img4);
}


img3 = cv2.scaleAdd(img2, 0.5, img1)

img4 = 0.5*img2+img1


配列を重みつきで加算する(2)。

次の例は、画像の引き算を行って、差がないときのバイアスのカウントとして128を足す例です。

書き方によって、処理時間が大幅に違うことがわかります。

実装が異なるので端数の扱いが違ってきますから、

処理結果の違いは確認してから、使いましょう。


ex_addWeighted.cpp

#include <opencv2/opencv.hpp>


int main(int agrc, char* argv[]){
std::string name1 = "../../data/lena.png";

cv::Mat img1=cv::imread(name1, 0);
cv::Mat img2;
cv::flip(img1, img2, 1);
cv::Mat img3;

int64 time0 = cv::getTickCount();
cv::addWeighted(img2, 0.5, img1, -0.5, 128, img3);
int64 time1 = cv::getTickCount();

cv::imwrite("addWeighted.png", img3);

cv::Size size=cv::Size(img1.cols, img1.rows);
cv::Mat img4 = cv::Mat::zeros(size, CV_8UC1);

int64 time2 = cv::getTickCount();
for(int y=0; y< img1.rows; y++){
for(int x=0; x< img1.cols; x++){
img4.at<uchar>(y,x) = 0.5*img2.at<uchar>(y,x) - 0.5* img1.at<uchar>(y,x) +128;
}
}
int64 time3 = cv::getTickCount();
cv::imwrite("addWeighted_loop.png", img4);

std::cout << (time1-time0)/cv::getTickFrequency() << std::endl;
std::cout << (time3-time2)/cv::getTickFrequency() << std::endl;

}



ex_addWeighted.py

import cv2

import numpy as np
import pylab as plt

img1 = cv2.imread("lena.png", 0)
img2 = img1[:, ::-1]+0

t1 = cv2.getTickCount()
img3 = cv2.addWeighted(img2, 0.5, img1, -0.5, 128)
t2 = cv2.getTickCount()

cv2.imwrite("img3.png", img3)
t3 = cv2.getTickCount()
img4 = 0.5*img2-0.5*img1 +128
t4 = cv2.getTickCount()
cv2.imwrite("img4.png", img4)

t5 = cv2.getTickCount()
img5 = (np.array(img2, dtype=np.int16) -np.array(img1,dtype=np.int16))/2 +128
t6 = cv2.getTickCount()
cv2.imwrite("img5.png", img4)

print (t2-t1)/cv2.getTickFrequency()
print (t4-t3)/cv2.getTickFrequency()
print (t6-t5)/cv2.getTickFrequency()

plt.figure(1)
plt.clf()
plt.subplot(2, 3, 1)
plt.imshow(img3)
plt.colorbar()
plt.subplot(2, 3, 2)
plt.imshow(img4)
plt.colorbar()
plt.subplot(2, 3, 3)
plt.imshow(img5)
plt.colorbar()
plt.subplot(2, 3, 5)
plt.imshow(img4-img3)
plt.colorbar()
plt.subplot(2, 3, 6)
plt.imshow(img5-img3)
plt.colorbar()
plt.show()



実行結果例.txt

0.000149026401895

0.00955493244547
0.00421297227616



条件付きで加算できる。

cv::add()

C++: void add(InputArray src1, InputArray src2, OutputArray dst, InputArray mask=noArray(), int dtype=-1)

Python: cv2.add(src1, src2[, dst[, mask[, dtype]]]) → dst

これも、2重ループ内のif文の出番をなくすことにつながる。


単位行列を設定する。

cv::setIdentity

cv2.setIdentity

>>> a=np.ones((3,3))

>>> cv2.setIdentity(a)
>>> a
array([[ 1., 0., 0.],
[ 0., 1., 0.],
[ 0., 0., 1.]])
>>>


条件式によって値を設定する。

setTo

mat.setTo(0)

//次の記述は、その後に示す2重ループのfor文に相当します。

//記述が簡潔になるので、処理の見通しがよくなります。

numpy の場合

img3[img1 > 128] = 255

OpenCV C++の場合

img3.setTo(255, img1>128);


ex_setTo.cpp

#include <opencv2/opencv.hpp>


int main(int agrc, char* argv[]){
std::string name1 = "lena.png";

cv::Mat img1=cv::imread(name1, 0);

cv::Mat img3 = cv::Mat::zeros(img1.size(), img1.type());
cv::Mat img4=cv::Mat::zeros(img1.size(), img1.type());

int64 time0 = cv::getTickCount();
img3.setTo(255, img1>128);

int64 time1 = cv::getTickCount();

std::cout << (time1-time0)/cv::getTickFrequency() << std::endl;

int64 time2 = cv::getTickCount();

for(int y =0; y<img1.rows; y++){
for(int x=0; x< img1.cols; x++){
if(img1.at<uchar>(y, x)>128){
img4.at<uchar>(y, x)=255;
}
}
}
int64 time3 = cv::getTickCount();
std::cout << (time3-time2)/cv::getTickFrequency() << std::endl;

cv::imwrite("img3.png", img3);
cv::imwrite("img4.png", img4);
}


$ ./ex_max

0.000244705

0.00290059

$

このように2重ループを使わない記述の方が、格段に高速です。


条件式によって値をコピーする。

copyTo

//次の記述は、その後に示す2重ループのfor文に相当します。

//記述が簡潔になるので、処理の見通しがよくなります。

src.copyTo(dstimg, maskimg);


ex_copyTo.cpp

#include <opencv2/opencv.hpp>


int main(int agrc, char* argv[]){
std::string name1 = "../../data/lena.png";
std::string name2 = "../../data/home.jpg";

cv::Mat img1=cv::imread(name1, 1);
cv::Mat img2=cv::imread(name2, 1);

cv::Mat gray1=cv::imread(name1, 0);

cv::Mat img3 = img2.clone();
cv::Mat img4 = img2.clone();

int64 time0 = cv::getTickCount();
img1.copyTo(img3, gray1>64);

int64 time1 = cv::getTickCount();

std::cout << (time1-time0)/cv::getTickFrequency() << std::endl;

int64 time2 = cv::getTickCount();

for(int y =0; y<img1.rows; y++){
for(int x=0; x< img1.cols; x++){
if(gray1.at<uchar>(y, x)>64){
img4.at<cv::Vec3b>(y, x)=img1.at<cv::Vec3b>(y, x);
}
}
}
int64 time3 = cv::getTickCount();
std::cout << (time3-time2)/cv::getTickFrequency() << std::endl;

cv::imwrite("ex_copyTo.png", img3);
cv::imwrite("ex_copyTo_loop.png", img4);
}


copyTo()の利用例

OpenCVでクロマキー


複数の画像の平均画像を作るには

Mat::convertTo()

- uint8の画像(配列)で型変換してから2枚の画像の重み付け平均を出す。uint8の画像のまま足し算するuint8からオーバーフローしてしまう。それを防止するためにfloatに型変換する。その後で重み付けの和をとり、元々のuint8に換算する。そのような無駄が標準で用意されている関数を使うと取り除ける。

image.convertTo(tmp, CV_32F, 1.0/255);

convertScaleAbs()


 reduce()の利用

reduce() を使うと2次元の配列をベクトルにする操作が、マルチコアやGPUで効果的に処理できそうに見えます。

(マルチコアの設定がどれだけ有効に働いているのかの検証はまだできていません。

numpyのビルドとcv2のビルドの設定による部分で動作状況が変わってくるとみられるので、

np.mean() とcv2.mean() について公平な評価になっていないことに注意してください。

cv2.mean()がcv2.reduce()を用いて自作したfloatMean()よりも格段に速いことに驚いた。

標準化されている関数は、かなり最適化されているのだろう。CPUやOSの種類、OSのビット数、CPUのコア数などによって結果が変わってくると思うので、この処理結果の時間については深入りしません。)


test_reduce.py

import cv2

import numpy as np
def intMean(testMat):
gg = cv2.reduce(testMat, 0, cv2.cv.CV_REDUCE_AVG, dtype=cv2.CV_32S)
jj = 0
for i in range(gg.shape[1]):
jj += gg[0, i]

jj /= gg.shape[1]
return jj

def floatMean(testMat):
gg = cv2.reduce(testMat, 0, cv2.cv.CV_REDUCE_AVG, dtype=cv2.CV_32F)
jj = 0
for i in range(gg.shape[1]):
jj += gg[0, i]

jj /= float(gg.shape[1])
return jj

testMat = cv2.imread("lena.png", 0);
[h, w] = testMat.shape[:2]

t1 = cv2.getTickCount()
for x in range(100):
jj = intMean(testMat)
t2 = cv2.getTickCount()
print jj

t3 = cv2.getTickCount()
for x in range(100):
jj = floatMean(testMat)
t4 = cv2.getTickCount()
print jj

t5 = cv2.getTickCount()
for x in range(100):
r = np.mean(testMat)
t6 = cv2.getTickCount()
print r

t7 = cv2.getTickCount()
for x in range(100):
r = cv2.mean(testMat)
t8 = cv2.getTickCount()
print r

print "time intMean = ", (t2-t1)/cv2.getTickFrequency()
print "time floatMean= ", (t4-t3)/cv2.getTickFrequency()
print "time np.mean= ", (t6-t5)/cv2.getTickFrequency()
print "time cv2.mean= ", (t8-t7)/cv2.getTickFrequency()



実行結果例.txt

100

100.312770844
100.312770844
(100.31277084350586, 0.0, 0.0, 0.0)
time intMean = 0.0234855946179
time floatMean= 0.0263815737886
time np.mean= 0.0440058280862
time cv2.mean= 0.00841697830397


標準偏差が必要になるときには平均も必要ですね。

cv::meanStdDev()

cv2.meanStdDev()


ex_mean.cpp

#include <opencv2/opencv.hpp>


int main(int agrc, char* argv[]){
std::string name1 = "lena.png";

cv::Mat img1=cv::imread(name1, 0);
int64 time0 = cv::getTickCount();
cv::Mat mean;
cv::Mat stddev;
cv::meanStdDev(img1, mean, stddev );
int64 time1 = cv::getTickCount();

std::cout << "mean= " << mean << std::endl;
std::cout << "stddev= " << stddev << std::endl;
std::cout << "mean= " << mean.at<double>(0,0) << std::endl;
std::cout << "stddev= " <<stddev.at<double>(0,0) << std::endl;

int64 time2 = cv::getTickCount();
double sum=0;
for(int y=0; y<img1.rows; y++){
for(int x=0; x< img1.cols; x++){
sum += img1.at<uchar>(y, x);
}
}

double meanv = sum /(img1.rows*img1.cols);
int64 time3 = cv::getTickCount();
std::cout <<"meanv=" << meanv << std::endl;

std::cout << (time1-time0)/cv::getTickFrequency() << std::endl;
std::cout << (time3-time2)/cv::getTickFrequency() << std::endl;

}


$ ./ex_mean

mean= [100.3127708435059]

stddev= [44.94230414402137]

mean= 100.313

stddev= 44.9423

meanv=100.313

0.000210359

0.00106834

$

平均の他に標準偏差も求めているのに、自前の2重ループよりも格段に速い結果になっています。


  • 0ではない値を数えるときに、2重ループを書いてしまうことはありがちですね。
    cv::countNonZero()
    cv2.countNonZero()


4重ループ? 実はfilter2D()の出番では?

画像処理で遅くなりがちの処理の1つに4重ループがあります。

もし、あなたのコードの中に4重ループがあったら、filter2D()で置き換えられるものかどうか検討してみてください。以下の記述が

http://docs.opencv.org/3.0-beta/modules/imgproc/doc/filtering.html

の中にありますように、kernelが大きいときにはDFT(離散フーリエ変換)を使って効率的な計算をしてくれます。


The function uses the DFT-based algorithm in case of sufficiently large kernels (~11 x 11 or larger) and the direct algorithm for small kernels.




並び替えの順序を知りたい

sortIdx()はnumpy.argsort() に似たものです。

次のように置き換えることができます。


ex_argsort.py

import numpy as np

a=np.array([ [ 6, 22, 14, 4, 22],
[10, 11, 18, 1, 4],
[1, 13, 7, 13, 0],
[23, 18, 1, 17, 3],
[ 3, 18, 0, 14, 23]])

print a

print "ROW|ASCENDING"
print np.argsort(a, axis=1)
print "ROW|DESCENDING"
print np.argsort(a, axis=1)[:, ::-1]
print "COLUMN|ASCENDING"
print np.argsort(a, axis=0)
print "COLUMN|DESCENDING"
print np.argsort(a, axis=0)[::-1, :]


[[ 6 22 14  4 22]

[10 11 18 1 4]
[ 1 13 7 13 0]
[23 18 1 17 3]
[ 3 18 0 14 23]]
ROW|ASCENDING
[[3 0 2 1 4]
[3 4 0 1 2]
[4 0 2 1 3]
[2 4 3 1 0]
[2 0 3 1 4]]
ROW|DESCENDING
[[4 1 2 0 3]
[2 1 0 4 3]
[3 1 2 0 4]
[0 1 3 4 2]
[4 1 3 0 2]]
COLUMN|ASCENDING
[[2 1 4 1 2]
[4 2 3 0 3]
[0 3 2 2 1]
[1 4 0 4 0]
[3 0 1 3 4]]
COLUMN|DESCENDING
[[3 0 1 3 4]
[1 4 0 4 0]
[0 3 2 2 1]
[4 2 3 0 3]
[2 1 4 1 2]]
>>>


ex_sortIdx.py

import numpy as np

import cv2

a=np.array([ [ 6, 22, 14, 4, 22],
[10, 11, 18, 1, 4],
[1, 13, 7, 13, 0],
[23, 18, 1, 17, 3],
[ 3, 18, 0, 14, 23]])

print a
print
print "COLUMN|ASCENDING"
print cv2.sortIdx(a, cv2.SORT_EVERY_ROW | cv2.SORT_ASCENDING)
print
print "COLUMN|DESCENDING"
print cv2.sortIdx(a, cv2.SORT_EVERY_ROW | cv2.SORT_DESCENDING)
print

print "COLUMN|ASCENDING"
print cv2.sortIdx(a, cv2.SORT_EVERY_COLUMN | cv2.SORT_ASCENDING)
print
print "COLUMN|DESCENDING"
print cv2.sortIdx(a, cv2.SORT_EVERY_COLUMN | cv2.SORT_DESCENDING)





cv::Mat のメソッドは、cv2では利用できない。

numpy.arrayのメソッド、numpyの関数で対応する機能を見つけて書き換える必要があります。

追記:


2重ループを単一のループにする。

numpyのarrayには reshape()があります。

2次元配列が1次元の配列になるようにreshape()しておいてから、

単一のfor ループを実行させるという手があります。

そうすると2重ループよりは、ループの範囲の判定が単純になるため、速くする効果はあるらしいです。

しかし、明示的にループを書くよりは、ライブラリに用意されている最適化済みの関数を使うことだと思っています。