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 不太熟悉,可以参考之前的文章: