目标

本教程的目标是向你展示如何使用OpenCV parallel_for_ 框架轻松地并行化你的代码。为了说明这个概念，我们将编写一个程序来绘制曼德勃罗集合，利用几乎所有可用的CPU负载。完整的教程代码在此。如果你想了解更多关于多线程的信息，你将不得不参考参考书或课程，因为本教程旨在保持简单。

前提条件

第一个前提条件是使用并行框架构建OpenCV。在OpenCV 3.2中，按以下顺序提供了以下并行框架：

英特尔线程构建块 (Intel Threading Building Blocks) (第三方库，应显式启用)
C= 并行C/C++编程语言扩展 (第三方库，应显式启用)
OpenMP (集成到编译器中，应显式启用)
APPLE GCD (系统范围，自动使用 (仅限APPLE))
Windows RT 并发 (系统范围，自动使用 (仅限Windows RT))
Windows 并发 (运行时的一部分，自动使用 (仅限Windows - MSVC++ >= 10))
Pthreads (如果可用)

正如你所看到的，OpenCV库可以使用多个并行框架。一些并行库是第三方库，必须在CMake中显式构建和启用(例如TBB，C=)；其他的库可通过平台自动获得(例如APPLE GCD)，但你应该能够直接或通过在CMake中启用选项并重新构建库来访问并行框架。

第二个(较弱的)前提条件与你想要完成的任务更相关，因为并非所有计算都适合/可以适应以并行方式运行。为保持简单，可以将任务拆分为多个没有内存依赖关系(没有可能的竞争条件)的基本操作的任务很容易并行化。计算机视觉处理通常很容易并行化，因为大多数情况下，一个像素的处理不依赖于其他像素的状态。

简单示例：绘制曼德勃罗集合

我们将使用绘制曼德勃罗集合的示例来展示如何从常规的顺序代码轻松地调整代码以并行化计算。

理论

曼德勃罗集合的定义是为了向数学家本华·曼德勃罗致敬而由数学家阿德里安·杜瓦迪命名的。它在数学领域之外也很有名，因为图像表示是分形类的一个例子，分形是一个在每个尺度上都显示重复模式的数学集合(更重要的是，曼德勃罗集合是自相似的，因为整个形状可以在不同的尺度上重复出现)。要了解更多信息，可以查看相应的维基百科文章。在这里，我们只介绍绘制曼德勃罗集合的公式(来自上面提到的维基百科文章)。

曼德勃罗集合是复平面中\( c \)的值的集合，对于二次映射的0的轨道：
\[\begin{cases} z_0 = 0 \\ z_{n+1} = z_n^2 + c \end{cases}\]

保持有界。也就是说，如果从\( z_0 = 0 \)开始并重复应用迭代，那么复数\( c \)是曼德勃罗集合的一部分，无论\( n \)有多大，\( z_n \)的绝对值都保持有界。这也可以表示为：
\[\limsup_{n\to\infty}|z_{n+1}|\leqslant2\]

伪代码

生成曼德勃罗集合表示的一个简单算法称为“逃逸时间算法”。对于渲染图像中的每个像素，我们使用递推关系测试复数在最大迭代次数下是否是有界的。不属于曼德勃罗集合的像素将很快逃逸，而我们假设在固定的最大迭代次数后像素位于集合中。较高的迭代次数将产生更详细的图像，但计算时间将相应增加。我们使用“逃逸”所需的迭代次数来描绘图像中的像素值。

对于屏幕上的每个像素(Px, Py)，执行：
{
x0 = 像素的缩放x坐标(缩放为位于曼德勃罗X尺度(-2, 1)内)
y0 = 像素的缩放y坐标(缩放为位于曼德勃罗Y尺度(-1, 1)内)
x = 0.0
y = 0.0
iteration = 0
max_iteration = 1000
while (x*x + y*y < 2*2 AND iteration < max_iteration) {
xtemp = x*x - y*y + x0
y = 2*x*y + y0
x = xtemp
iteration = iteration + 1
  }
color = palette[iteration]
plot(Px, Py, color)
}

为了将伪代码和理论联系起来，我们有：

\( z = x + iy \)
\( z^2 = x^2 + i2xy - y^2 \)
\( c = x_0 + iy_0 \)

在这个图中，我们回忆一下，复数的实部在x轴上，虚部在y轴上。你可以看到，如果我们在特定位置放大，整个形状可以重复可见。

实现

逃逸时间算法实现

int mandelbrot(const complex<float> &z0, const int max)
{
complex<float> z = z0;
    for (int t = 0; t < max; t++)
    {
        如果 (z.real() * z.real() + z.imag() * z.imag() > 4.0f) 返回 t;
        z = z * z + z0;
    }
 
    返回 max;
}

这里，我们使用 std::complex 模板类来表示复数。此函数执行测试以检查像素是否在集合中，并返回“逃逸”迭代次数。

顺序 Mandelbrot 实现

void sequentialMandelbrot(Mat &img, const float x1, const float y1, const float scaleX, const float scaleY)
{
    for (int i = 0; i < img.rows; i++)
    {
        for (int j = 0; j < img.cols; j++)
        {
            float x0 = j / scaleX + x1;
            float y0 = i / scaleY + y1;
 
complex<float> z0(x0, y0);
            uchar value = (uchar) mandelbrotFormula(z0);
img.ptr<uchar>(i)[j] = value;
        }
    }
}

在这个实现中，我们顺序迭代渲染图像中的像素，执行测试以检查像素是否可能属于 Mandelbrot 集。

另一件事是将像素坐标转换为 Mandelbrot 集空间，方法是

    Mat mandelbrotImg(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);

最后，为了将灰度值赋给像素，我们使用以下规则：

如果像素达到最大迭代次数（假设像素在 Mandelbrot 集中），则像素为黑色；
否则，我们根据逃逸迭代次数分配灰度值，并将其缩放以适应灰度范围。

int mandelbrotFormula(const complex<float> &z0, const int maxIter=500) {
    int value = mandelbrot(z0, maxIter);
    如果(maxIter - value == 0)
    {
        返回 0;
    }
 
    返回 cvRound(sqrt(value / (float) maxIter) * 255);
}

使用线性比例变换不足以感知灰度变化。为了克服这个问题，我们将使用平方根比例变换来增强感知（借鉴自 Jeremy D. Frens 在其博客文章中的方法）： \( f \left( x \right) = \sqrt{\frac{x}{\text{maxIter}}} \times 255 \)

绿色曲线对应于简单的线性比例变换，蓝色曲线对应于平方根比例变换，您可以观察到在这些位置查看斜率时，最低值将如何得到提升。

并行 Mandelbrot 实现

查看顺序实现，我们可以注意到每个像素都是独立计算的。为了优化计算，我们可以通过利用现代处理器的多核架构来并行执行多个像素计算。为了轻松实现这一点，我们将使用 OpenCV cv::parallel_for_ 框架。

class ParallelMandelbrot : public ParallelLoopBody
{
public:
ParallelMandelbrot (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 CV_OVERRIDE
    {
        for (int r = range.start; r < range.end; r++)
        {
            int i = r / m_img.cols;
            int j = r % m_img.cols;
 
            float x0 = j / m_scaleX + m_x1;
            float y0 = i / m_scaleY + m_y1;
 
complex<float> z0(x0, y0);
            uchar value = (uchar) mandelbrotFormula(z0);
m_img.ptr<uchar>(i)[j] = value;
        }
    }
 
ParallelMandelbrot& operator=(const ParallelMandelbrot &) {
        返回 *this;
    };
 
private:
    Mat &m_img;
    float m_x1;
    float m_y1;
    float m_scaleX;
    float m_scaleY;
};

首先，声明一个从 cv::ParallelLoopBody 继承的自定义类，并覆盖 virtual void operator ()(const cv::Range& range) const。

operator ()中的范围表示单个线程将处理的像素子集。此分割是自动完成的，以平均分配计算负载。我们必须将像素索引坐标转换为二维[row, col]坐标。还要注意，我们必须保留对mat图像的引用，以便能够就地修改图像。

并行执行使用以下方式调用：

ParallelMandelbrot parallelMandelbrot(mandelbrotImg, x1, y1, scaleX, scaleY);

parallel_for_(Range(0, mandelbrotImg.rows*mandelbrotImg.cols), parallelMandelbrot);

这里，范围表示要执行的总操作数，也就是图像中的总像素数。要设置线程数，可以使用：cv::setNumThreads。您也可以在cv::parallel_for_中使用nstripes参数指定分割数量。例如，如果您的处理器有4个线程，则设置cv::setNumThreads(2)或设置nstripes=2应该与默认情况下使用所有可用的处理器线程但只将工作负载拆分为两个线程相同。

注意: C++ 11标准允许通过摆脱ParallelMandelbrot类并将其替换为lambda表达式来简化并行实现。

    parallel_for_(Range(0, mandelbrotImg.rows*mandelbrotImg.cols), [&](const Range& range){
for (int r = range.start; r < range.end; r++)
        {
int i = r / mandelbrotImg.cols;
int j = r % mandelbrotImg.cols;
 
float x0 = j / scaleX + x1;
float y0 = i / scaleY + y1;
 
complex<float> z0(x0, y0);
uchar value = (uchar) mandelbrotFormula(z0);
mandelbrotImg.ptr<uchar>(i)[j] = value;
        }
    });

结果

您可以在这里找到完整的教程代码here。并行实现的性能取决于您使用的CPU类型。例如，在4核/8线程CPU上，您可以期望大约6.9X的加速。有很多因素可以解释为什么我们没有达到近8X的加速。主要原因可能是由于：

创建和管理线程的开销；
并行运行的后台进程；
每个核心4个硬件核心和2个逻辑线程与8个硬件核心之间的差异。

教程代码生成的图像（您可以修改代码以使用更多迭代，并根据逃逸迭代分配像素颜色，并使用调色板来获得更美观的图像）：

目录

目标

前提条件

简单示例：绘制曼德勃罗集合

理论

伪代码

实现

逃逸时间算法实现

顺序 Mandelbrot 实现

并行 Mandelbrot 实现

结果