注意：検証が不十分な部分があり、そのままでは動作しない記述です。
また、プログラムの例は、見本的なプログラムとは程遠いものであることをお許しください。
画像処理では旧来より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文でループを書かなくてはならないことが減るはずです。
とは言っても、ループを書かなくてはならないこともあるはずなので、次の記事を紹介しておく。
以下に示す2重ループのコードは、高速化を意識しない記述になっています。
この文章の中では、そのようなコードよりもOpenCVの標準の関数の方が格段に速い結果を与えていることを示すものです。その理由で、「画像処理で2重ループをなるべく書くな」ということを主張します。
悪いコードの見本
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){
グレースケール画像に対する処理
}
}
}
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])
#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);
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
差の絶対値を求めることは多いですね。
地道に計算時間を節約。
#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);
}
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()
[[ 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
を使おう。これらの関数を用いた場合には、マルチコアへの最適化がされている可能性が高い。
さらに、最大値・最小値で正規化する場合には、既に関数が標準で用意されている。
配列を重みつきで加算する。
これはBLASなどのライブラリには
SAXPY y = ax + y
などがあるのに対応している関数のようです。
dst(I)= scalesrc1(I) + src2(I)
cv::scaleAdd
cv2.scaleAdd
次の２つの記述は、ほぼ同等の結果になる。ただし端数の扱いが異なるため一致しない。
scaleAddを用いた方がpythonインタプリタ上の計測では格段に速い。
#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を足す例です。
書き方によって、処理時間が大幅に違うことがわかります。
実装が異なるので端数の扱いが違ってきますから、
処理結果の違いは確認してから、使いましょう。
#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;
}
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()
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](http://docs.opencv.org/3.0-beta/modules/core/doc/basic_structures.html#Mat& Mat::setTo(InputArray value, InputArray mask))
mat.setTo(0)
//次の記述は、その後に示す2重ループのfor文に相当します。
//記述が簡潔になるので、処理の見通しがよくなります。
numpy の場合
img3[img1 > 128] = 255
OpenCV C++の場合
img3.setTo(255, img1>128);
#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]
(http://docs.opencv.org/3.0-beta/modules/core/doc/basic_structures.html?highlight=copyto#void Mat::copyTo(OutputArray m) const)
//次の記述は、その後に示す2重ループのfor文に相当します。
//記述が簡潔になるので、処理の見通しがよくなります。
src.copyTo(dstimg, maskimg);
#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のコア数などによって結果が変わってくると思うので、この処理結果の時間については深入りしません。）
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()
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()
#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()の出番では？
画像処理で遅くなりがちの処理の１つに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 11or larger) and the direct algorithm for small kernels.
並び替えの順序を知りたい
sortIdx()はnumpy.argsort() に似たものです。
次のように置き換えることができます。
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]]
>>>
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の関数で対応する機能を見つけて書き換える必要があります。
追記：
２重ループを単一のループにする。
numpyのarrayには reshape()があります。
２次元配列が１次元の配列になるようにreshape()しておいてから、
単一のfor ループを実行させるという手があります。
そうすると２重ループよりは、ループの範囲の判定が単純になるため、速くする効果はあるらしいです。
しかし、明示的にループを書くよりは、ライブラリに用意されている最適化済みの関数を使うことだと思っています。