OpenCV 并行计算函数 parallel_for_ 的使用
Reference：
- OpenCV并行加速Parallel_for_与ParallelLoopBody教程
- opencv 并行计算函数 parallel_for_的使用
- 官网：Parallel Processing
- 高翔，张涛 《视觉SLAM十四讲》
相关文章：
在使用 OpenCV 的过程中，对图片的处理计算量还是很大的，所以在实施运行的程序中如何高效的计算会节省很多时间。现有的方法有很多，如 OpenMp, TBB, OpenCL 当然还有 Nvidia 的 CUDA。
CUDA 是个好东西，但是并不太适合毫秒级别的程序运行，单一张图片在 cpu 和 gpu 之间的传输时间就已经达到 300ms（使用 opencv 的cuda 库函数）；在TX2上直接对cuda进行编程，数据的传输也是在 50ms（不包含初始化）以上，根本不能拿来做实时的运算。所以如何在cpu上更加高效的计算变得尤为重要。偶然间发现了 OpenCV 的并行计算函数 parallel_for_，它整合了上述的多个组件。
对于一些基本的循环运算，如果我们直接使用循环，即便是使用指针，运算效率也不高，如果我们使用并行计算，会大大提升运算效率，OpenCV 里面的很多运算都是使用了并行加速的。在 OpenCV 3.2 中，并行框架按照以下顺序提供：
- 英特尔线程构建块（第三方库，应显式启用），如TBB(Thread Building Blocks)
- OpenMP（集成到编译器，应该被显式启用）
- APPLE GCD（系统范围广，自动使用（仅限APPLE））
- Windows RT并发（系统范围，自动使用（仅Windows RT））
- Windows并发（运行时的一部分，自动使用（仅限Windows） - MSVC ++> = 10））
- Pthreads（如果有的话）
可以看到，OpenCV 库中可以使用多个并行框架：
一些并行库是第三方库，必须在 CMake（例如TBB）中进行显式构建和启用，其他可以自动与平台（例如 APPLE GCD）一起使用，但是您应该可以使用这些库来访问并行框架直接或通过启用CMake中的选项并重建库；
第二个（弱）前提条件与要实现的任务更相关，因为并不是所有的计算都是合适的/可以被平行地运行。为了保持简单，可以分解成多个基本操作而没有内存依赖性（无可能的竞争条件）的任务很容易并行化。计算机视觉处理通常易于并行化，因为大多数时间一个像素的处理不依赖于其他像素的状态。
Parallel_for_
本文主要描述 Parallel_for_与 ParallelLoopBody 的使用方法，需要使用 <opencv2/core/utility.hpp> 这个头文件。
Parallel_for_ 的介绍为 parallel data processor，有两种使用方式：
void cv::parallel_for_ (const Range &range, const ParallelLoopBody &body, double nstripes=-1.)
static void cv::parallel_for_ (const Range &range, std::function< void(const Range &)> functor, double nstripes=-1.)
1. 先从构造函数的运算符重载说起
我们在调用函数的时候，实际上是使用了括号 ( ) 运算符，构造函数也是普通的函数，所以也用到了括号运算符，那如果想要重载这个括号 ( ) 运算符怎么做呢？先来看一个简单的小例子：
class Animal
{
public:
Animal(float weight_,int age_)
{
weight = weight_;
age = age_;
}
void operator()(float height) const //重载操作符（）
{
height = 100.0;
std::cout << "the age is : " << age << std::endl;
std::cout << "the weight is : " << weight << std::endl;
std::cout << "the height is : " << height << std::endl;
}
private:
float weight;
int age;
};
现在调用：
int main(int argc, char* argv[])
{
myanimal::Animal animal(50.0,25); //创建对象
animal(100.0); //通过对象调用重载的括号 () 运算符
getchar();
return 0;
}
/*
the age is : 25
the weight is : 50
the height is : 100
*/
总结：
（1）重载的括号运算符就像一个对象的方法一杨，依旧是通过对象去调用，调用的方式为 “对象名(参数列表)” 这样的形式；
（2）重载括号运算符的一般操作为 “返回类型 operator(参数列表)” ，后面的const可以不要，参数列表可以使任意的，
没有参数，则调用方式为：obj（）
一个参数，则调用方式为：obj（参数1）
多个参数，则调用方式为：obj（参数1，参数2，…）
2. Parallel_for_ 结合 ParallelLoopBody 使用的一般步骤
使用步骤一般遵循三步走的原则
（1）第一步：自定义一个类或者是一个结构体，使这个结构体或者是类继承自 ParallelLoopBody 类，如下：
class MyParallelClass : public ParallelLoopBody
{}
struct MyParallelStruct : public ParallelLoopBody
{}
（2）第二步：在自定义的类或者是结构体中，重写括号运算符（ ），注意：虽然前面讲括号运算符重载可以接受任意数量的参数，但是这里只能接受一个 Range 类型的参数（这是与一般的重载不一样的地方），因为后面的parallel_for_需要使用，如下：
void operator()(const Range& range)
{
//在这里面进行“循环操作”
}
（3）第三步：使用 parallel_for_ 进行并行处理
再看一下parallel_for_的函数原型
CV_EXPORTS void parallel_for_(const Range& range, //重载的括号运算符里面的参数，是一个cv::Range类型
const ParallelLoopBody& body, //自己实现的从ParallelLoopBody类继承的类或者是结构体对象
double nstripes=-1.);
使用方式如下：
parallel_for_(Range(start, end), MyParallelClass(构造函数列表));
//Range(start, end) 就是一个Range对象
//MyParallelClass(构造函数列表) 就是一个继承自ParallelLoopBody的类的对象
常规的使用方式应该像如下这样才对。这样写本身没什么问题，但是这样的执行方式，在括号重载运算符里面的内容是按照顺序执行的，并没有并发处理。
MyParallelClass obj = MyParallelClass(构造函数列表)); //构造对象
obj(Range(start, end)); //调用
而上面所述的使用方式:
parallel_for_(Range(start, end), MyParallelClass(构造函数列表));
没有显式的调用重载的括号运算符，但实际上是隐式调用的，并且以并发方式处理重载运算里面的内容。
3. Parallel_for_结合ParallelLoopBody的加速效果实验
3.1 自定义类实现
任务描述：我要定义两个Mat矩阵的逐元素乘积，如下所示
（1）自定义一个类继承自 ParallelLoopBody，并且重载括号运算
class ParallelAdd : public ParallelLoopBody//参考官方给出的answer，构造一个并行的循环体类
{
public:
ParallelAdd(Mat& _src1,Mat& _src2,Mat _result) //构造函数
{
src1 = _src1;
src2 = _src2;
result = _result;
CV_Assert((src1.rows == src2.rows) && (src1.cols == src2.cols));
rows = src1.rows;
cols = src1.cols;
}
void operator()(const Range& range) const //重载操作符（）
{
int step = (int)(result.step / result.elemSize1());//获取每一行的元素总个数（相当于cols*channels，等同于step1)
for (int col = range.start; col < range.end; ++col)
{
float * pData = (float*)result.col(col).data;
float * p1 = (float*)src1.col(col).data;
float * p2 = (float*)src2.col(col).data;
for (int row = 0; row < result.rows; ++row)
pData[row*step] = p1[row*step] * p2[row*step];
}
}
private:
Mat src1;
Mat src2;
Mat result;
int rows;
int cols;
};
可见重载的运算符里面是一个耗时操作，现在定义两种方式来实现两个 Mat 的逐元素乘积，一种是普通的逐元素处理，另一种是使用parallel进行并发处理，分别通过两个函数完成，如下：
//直接通过obj()形式调用，不采用并发处理
void testParallelClassWithFor(Mat _src1,Mat _src2,Mat result)
{
result = Mat(_src1.rows, _src1.cols, _src1.type());
int step = (int)(result.step / result.elemSize1());
int totalCols = _src1.cols;
ParallelAdd add = ParallelAdd(_src1, _src2, result);
add(Range(0, totalCols)); //直接调用，没有并发
}
void testParallelClassTestWithParallel_for_(Mat _src1,Mat _src2,Mat result)
{
result = Mat(_src1.rows, _src1.cols, _src1.type());
int step = (int)(result.step / result.elemSize1());
int totalCols = _src1.cols;
parallel_for_(Range(0, totalCols), ParallelAdd(_src1,_src2,result)); //隐式调用，并发
}
现在开始测试耗时对比:
int main(int argc, char* argv[])
{
Mat testInput1 = Mat::ones(6400, 5400, CV_32F);
Mat testInput2 = Mat::ones(6400, 5400, CV_32F);
Mat result1, result2, result3;
clock_t start, stop;
//****************测试耗时对比****************************
start = clock();
testParallelClassWithFor(testInput1, testInput2, result1);
stop = clock();
cout << "Running time using \'general for \':" << (double)(stop - start) / CLOCKS_PER_SEC * 1000 << "ms" << endl;
start = clock();
testParallelClassWithParallel_for_(testInput1, testInput2, result2);
stop = clock();
cout << "Running time using \'parallel for \':" << (double)(stop - start) / CLOCKS_PER_SEC * 1000 << "ms" << endl;
start = clock();
result3 = testInput1.mul(testInput2);
stop = clock();
cout << "Running time using \'mul function \':" << (double)(stop - start) / CLOCKS_PER_SEC * 1000 << "ms" << endl;
return 0;
}
/*
Running time using 'general for ':645ms
Running time using 'parallel for ':449ms
Running time using 'mul function ':70ms
*/
总结：我们可以看见，使用 parallel_for_ 并发的方式的确比直接调用快一些，快了将近 200ms，但是依旧没有使用 OpenCV 自带的标准函数 mul 函数速度快，因为，OpenCV 实现的函数库不仅仅经过了并行处理，还是用了更强大的底层优化，所以，只要是 OpenCV 自己带的方法，一般都是优先使用，除非自己写的比 OpenCV 的还牛逼一些。
3.2 自定义结构体实现
任务描述：现在定义一个并行运算的结构体，实现Mat逐元素的三次方运算
（1）自定义一个结构体继承自 ParallelLoopBody，并且重载括号运算，如下：
struct ParallelPow:ParallelLoopBody//构造一个供parallel_for使用的循环结构体
{
Mat* src; //结构体成员，一个Mat类型的指针
ParallelPow(Mat& _src)//struct 结构体构造函数
{
src = &_src;
}
void operator()(const Range& range) const
{
Mat& result = *src;
int step = (int)(result.step / result.elemSize1());
for (int col = range.start; col < range.end; ++col)
{
float* pData = (float*)result.col(col).data;
for (int row = 0; row < result.rows; ++row)
pData[row*step] = std::pow(pData[row*step], 3); //逐元素求立方
}
}
};
下面定义两个函数，一个是直接通过 for 循环逐元素进行立方运算，一个是通过 parallel_for_ 并发运算的,通过两个函数实现，如下所示：
void testParallelStructWithFor(Mat _src)
{
int totalCols = _src.cols;
ParallelPow obj = ParallelPow(_src);
obj(Range(0, totalCols));
}
void testParallelStructWithParallel_for(Mat _src)
{
int totalCols = _src.cols;
parallel_for_(Range(0, totalCols), ParallelPow(_src));
}
现在开始测试耗时对比:
int main(int argc, char* argv[])
{
Mat testInput1 = Mat::ones(6400, 5400, CV_32F);
Mat testInput2 = Mat::ones(6400, 5400, CV_32F);
Mat result1, result2, result3;
clock_t start, stop;
start = clock();
testParallelStructWithFor(testInput1);
stop = clock();
cout << "Running time using \'general for \':" << (double)(stop - start) / CLOCKS_PER_SEC * 1000 << "ms" << endl;
start = clock();
testParallelStructWithParallel_for(testInput1);
stop = clock();
cout << "Running time using \'parallel for \':" << (double)(stop - start) / CLOCKS_PER_SEC * 1000 << "ms" << endl;
start = clock();
testInput1.mul(testInput1).mul(testInput1);
stop = clock();
cout << "Running time using \'mul function \':" << (double)(stop - start) / CLOCKS_PER_SEC * 1000 << "ms" << endl;
return 0;
}
/*
Running time using 'general for ':881ms
Running time using 'parallel for ':195ms
Running time using 'mul function ':76ms
*/
可以发现，并行运算效率有着显著提升，但是相较于OpenCV的标准实现，依旧偏慢。因此在能够使用OpenCV标准函数实现的时候，尽量不要再自己定义运算，标准的OpenCV函数式经过优化了的，运算效率很高。
4. Parallel_for_不结合与结合ParallelLoopBody比较
事实上在c++11中，不一定要用一个类或结构体去继承并行计算循环体类（ParallelLoopBody），可以使用函数来替代。两种方式的比较在下面可见：
class Parallel_My : public ParallelLoopBody
{
public:
Parallel_My (Mat &img, const float x1, const float y1, const float scaleX, const float scaleY)
: m_img(img), m_x1(x1), m_y1(y1), m_scaleX(scaleX), m_scaleY(scaleY)
{
}
virtual void operator ()(const Range& range) const
{
for (int r = range.start; r < range.end; r++) //process of for loop
{
/***
***/
}
}
Parallel_My& operator=(const Parallel_My &) {
return *this;
};
private:
Mat &m_img;
float m_x1;
float m_y1;
float m_scaleX;
float m_scaleY;
};
int main()
{
//! [mandelbrot-transformation]
Mat Img(4800, 5400, CV_8U);
float x1 = -2.1f, x2 = 0.6f;
float y1 = -1.2f, y2 = 1.2f;
float scaleX = mandelbrotImg.cols / (x2 - x1);
float scaleY = mandelbrotImg.rows / (y2 - y1);
#ifdef CV_CXX11 //method one
parallel_for_(Range(0, Img.rows*tImg.cols), [&](const Range& range)
{
for (int r = range.start; r < range.end; r++) //这是需要并行计算的for循环
{
}
});
#else //method two
Parallel_My parallel_my0(Img, x1, y1, scaleX, scaleY);
parallel_for_(Range(0, Img.rows*Img.cols), parallel_my0);
#endif
}
5. 单层光流（补充内容）
上面讲的应该比较清楚了，可以看一个关于光流的示例：
class OpticalFlowTracker {
public:
OpticalFlowTracker(
const Mat &img1_,
const Mat &img2_,
const vector<KeyPoint> &kp1_,
vector<KeyPoint> &kp2_,
vector<bool> &success_,
bool inverse_ = true, bool has_initial_ = false) :
img1(img1_), img2(img2_), kp1(kp1_), kp2(kp2_), success(success_), inverse(inverse_),
has_initial(has_initial_) {}
void calculateOpticalFlow(const Range &range);
private:
const Mat &img1;
const Mat &img2;
const vector<KeyPoint> &kp1;
vector<KeyPoint> &kp2;
vector<bool> &success;
bool inverse = true;
bool has_initial = false;
};
void OpticalFlowSingleLevel(
const Mat &img1,
const Mat &img2,
const vector<KeyPoint> &kp1,
vector<KeyPoint> &kp2,
vector<bool> &success,
bool inverse, bool has_initial) {
kp2.resize(kp1.size());
success.resize(kp1.size());
OpticalFlowTracker tracker(img1, img2, kp1, kp2, success, inverse, has_initial);
parallel_for_(Range(0, kp1.size()),
std::bind(&OpticalFlowTracker::calculateOpticalFlow, &tracker, placeholders::_1));
}
void OpticalFlowTracker::calculateOpticalFlow(const Range &range) {
// parameters
int half_patch_size = 4;
int iterations = 10;
for (size_t i = range.start; i < range.end; i++) {
auto kp = kp1[i];
double dx = 0, dy = 0; // dx,dy need to be estimated
if (has_initial) {
dx = kp2[i].pt.x - kp.pt.x;
dy = kp2[i].pt.y - kp.pt.y;
}
double cost = 0, lastCost = 0;
bool succ = true; // indicate if this point succeeded
// Gauss-Newton iterations
Eigen::Matrix2d H = Eigen::Matrix2d::Zero(); // hessian
Eigen::Vector2d b = Eigen::Vector2d::Zero(); // bias
Eigen::Vector2d J; // jacobian
for (int iter = 0; iter < iterations; iter++) {
if (inverse == false) {
H = Eigen::Matrix2d::Zero();
b = Eigen::Vector2d::Zero();
} else {
// only reset b
b = Eigen::Vector2d::Zero();
}
cost = 0;
// compute cost and jacobian
for (int x = -half_patch_size; x < half_patch_size; x++)
for (int y = -half_patch_size; y < half_patch_size; y++) {
double error = GetPixelValue(img1, kp.pt.x + x, kp.pt.y + y) -
GetPixelValue(img2, kp.pt.x + x + dx, kp.pt.y + y + dy);; // Jacobian
if (inverse == false) {
J = -1.0 * Eigen::Vector2d(
0.5 * (GetPixelValue(img2, kp.pt.x + dx + x + 1, kp.pt.y + dy + y) -
GetPixelValue(img2, kp.pt.x + dx + x - 1, kp.pt.y + dy + y)),
0.5 * (GetPixelValue(img2, kp.pt.x + dx + x, kp.pt.y + dy + y + 1) -
GetPixelValue(img2, kp.pt.x + dx + x, kp.pt.y + dy + y - 1))
);
} else if (iter == 0) {
// in inverse mode, J keeps same for all iterations
// NOTE this J does not change when dx, dy is updated, so we can store it and only compute error
J = -1.0 * Eigen::Vector2d(
0.5 * (GetPixelValue(img1, kp.pt.x + x + 1, kp.pt.y + y) -
GetPixelValue(img1, kp.pt.x + x - 1, kp.pt.y + y)),
0.5 * (GetPixelValue(img1, kp.pt.x + x, kp.pt.y + y + 1) -
GetPixelValue(img1, kp.pt.x + x, kp.pt.y + y - 1))
);
}
// compute H, b and set cost;
b += -error * J;
cost += error * error;
if (inverse == false || iter == 0) {
// also update H
H += J * J.transpose();
}
}
// compute update
Eigen::Vector2d update = H.ldlt().solve(b);
if (std::isnan(update[0])) {
// sometimes occurred when we have a black or white patch and H is irreversible
cout << "update is nan" << endl;
succ = false;
break;
}
if (iter > 0 && cost > lastCost) {
break;
}
// update dx, dy
dx += update[0];
dy += update[1];
lastCost = cost;
succ = true;
if (update.norm() < 1e-2) {
// converge
break;
}
}
success[i] = succ;
// set kp2
kp2[i].pt = kp.pt + Point2f(dx, dy);
}
}
代码在 calculateOpticalFlow 函数中实现了单层光流函数，其中调用了 cv::parallel_for_ 并行调用 OpticalFlowTracker::calculateOpticalFlow 该函数计算指定范围内特征点的光流。这个并行 for 循环内部是 Intel tbb 库实现的，我们只需按照其接口，将函数本体定义出来，然后将函数作为 std::function 对象传递给它。
如果对 std::bind 不太熟悉，可以参考之前的文章：
2871
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
