人家吧发过不少飞快高斯模糊算法.俺呢发过不少很快高斯模糊算法.原和链接。

void CalGaussianCoeff(float sigma, float * a0, float * a1, float * a2, float * a3, float * b1, float * b2, float * cprev, float * cnext) {
    float alpha, lamma, k;

    if (sigma < 0.5f)
        sigma = 0.5f;
    alpha = (float)exp((0.726) * (0.726)) / sigma;
    lamma = (float)exp(-alpha);
    *b2 = (float)exp(-2 * alpha);
    k = (1 - lamma) * (1 - lamma) / (1 + 2 * alpha * lamma - (*b2));
    *a0 = k; *a1 = k * (alpha - 1) * lamma;
    *a2 = k * (alpha + 1) * lamma;
    *a3 = -k * (*b2);
    *b1 = -2 * lamma;
    *cprev = (*a0 + *a1) / (1 + *b1 + *b2);
    *cnext = (*a2 + *a3) / (1 + *b1 + *b2);
}

void gaussianHorizontal(unsigned char * bufferPerLine, unsigned char * lpRowInitial, unsigned char  * lpColumn, int width, int height, int Channels, int Nwidth, float a0a1, float a2a3, float b1b2, float  cprev, float cnext)
{
    int HeightStep = Channels*height;
    int WidthSubOne = width - 1;
    if (Channels == 3)
    {
        float prevOut[3];
        prevOut[0] = (lpRowInitial[0] * cprev);
        prevOut[1] = (lpRowInitial[1] * cprev);
        prevOut[2] = (lpRowInitial[2] * cprev);
        for (int x = 0; x < width; ++x) {
            prevOut[0] = ((lpRowInitial[0] * (a0a1)) - (prevOut[0] * (b1b2)));
            prevOut[1] = ((lpRowInitial[1] * (a0a1)) - (prevOut[1] * (b1b2)));
            prevOut[2] = ((lpRowInitial[2] * (a0a1)) - (prevOut[2] * (b1b2)));
            bufferPerLine[0] = prevOut[0];
            bufferPerLine[1] = prevOut[1];
            bufferPerLine[2] = prevOut[2];
            bufferPerLine += Channels;
            lpRowInitial += Channels;
        }
        lpRowInitial -= Channels;
        lpColumn += HeightStep * WidthSubOne;
        bufferPerLine -= Channels;
        prevOut[0] = (lpRowInitial[0] * cnext);
        prevOut[1] = (lpRowInitial[1] * cnext);
        prevOut[2] = (lpRowInitial[2] * cnext);

        for (int x = WidthSubOne; x >= 0; --x) {
            prevOut[0] = ((lpRowInitial[0] * (a2a3)) - (prevOut[0] * (b1b2)));
            prevOut[1] = ((lpRowInitial[1] * (a2a3)) - (prevOut[1] * (b1b2)));
            prevOut[2] = ((lpRowInitial[2] * (a2a3)) - (prevOut[2] * (b1b2)));
            bufferPerLine[0] += prevOut[0];
            bufferPerLine[1] += prevOut[1];
            bufferPerLine[2] += prevOut[2];
            lpColumn[0] = bufferPerLine[0];
            lpColumn[1] = bufferPerLine[1];
            lpColumn[2] = bufferPerLine[2];
            lpRowInitial -= Channels;
            lpColumn -= HeightStep;
            bufferPerLine -= Channels;
        }
    }
    else if (Channels == 4)
    {
        float prevOut[4];

        prevOut[0] = (lpRowInitial[0] * cprev);
        prevOut[1] = (lpRowInitial[1] * cprev);
        prevOut[2] = (lpRowInitial[2] * cprev);
        prevOut[3] = (lpRowInitial[3] * cprev);
        for (int x = 0; x < width; ++x) {
            prevOut[0] = ((lpRowInitial[0] * (a0a1)) - (prevOut[0] * (b1b2)));
            prevOut[1] = ((lpRowInitial[1] * (a0a1)) - (prevOut[1] * (b1b2)));
            prevOut[2] = ((lpRowInitial[2] * (a0a1)) - (prevOut[2] * (b1b2)));
            prevOut[3] = ((lpRowInitial[3] * (a0a1)) - (prevOut[3] * (b1b2)));

            bufferPerLine[0] = prevOut[0];
            bufferPerLine[1] = prevOut[1];
            bufferPerLine[2] = prevOut[2];
            bufferPerLine[3] = prevOut[3];
            bufferPerLine += Channels;
            lpRowInitial += Channels;
        }
        lpRowInitial -= Channels;
        lpColumn += HeightStep * WidthSubOne;
        bufferPerLine -= Channels;

        prevOut[0] = (lpRowInitial[0] * cnext);
        prevOut[1] = (lpRowInitial[1] * cnext);
        prevOut[2] = (lpRowInitial[2] * cnext);
        prevOut[3] = (lpRowInitial[3] * cnext);

        for (int x = WidthSubOne; x >= 0; --x) {
            prevOut[0] = ((lpRowInitial[0] * a2a3) - (prevOut[0] * b1b2));
            prevOut[1] = ((lpRowInitial[1] * a2a3) - (prevOut[1] * b1b2));
            prevOut[2] = ((lpRowInitial[2] * a2a3) - (prevOut[2] * b1b2));
            prevOut[3] = ((lpRowInitial[3] * a2a3) - (prevOut[3] * b1b2));
            bufferPerLine[0] += prevOut[0];
            bufferPerLine[1] += prevOut[1];
            bufferPerLine[2] += prevOut[2];
            bufferPerLine[3] += prevOut[3];
            lpColumn[0] = bufferPerLine[0];
            lpColumn[1] = bufferPerLine[1];
            lpColumn[2] = bufferPerLine[2];
            lpColumn[3] = bufferPerLine[3];
            lpRowInitial -= Channels;
            lpColumn -= HeightStep;
            bufferPerLine -= Channels;
        }
    }
    else if (Channels == 1)
    {
        float prevOut = (lpRowInitial[0] * cprev);

        for (int x = 0; x < width; ++x) {
            prevOut = ((lpRowInitial[0] * (a0a1)) - (prevOut  * (b1b2)));
            bufferPerLine[0] = prevOut;
            bufferPerLine += Channels;
            lpRowInitial += Channels;
        }
        lpRowInitial -= Channels;
        lpColumn += HeightStep*WidthSubOne;
        bufferPerLine -= Channels;

        prevOut = (lpRowInitial[0] * cnext);

        for (int x = WidthSubOne; x >= 0; --x) {
            prevOut = ((lpRowInitial[0] * a2a3) - (prevOut  * b1b2));
            bufferPerLine[0] += prevOut;
            lpColumn[0] = bufferPerLine[0];
            lpRowInitial -= Channels;
            lpColumn -= HeightStep;
            bufferPerLine -= Channels;
        }
    }
}

void gaussianVertical(unsigned char * bufferPerLine, unsigned char * lpRowInitial, unsigned char * lpColInitial, int height, int width, int Channels, float a0a1, float a2a3, float b1b2, float  cprev, float  cnext) {

    int WidthStep = Channels*width;
    int HeightSubOne = height - 1;
    if (Channels == 3)
    {
        float prevOut[3];
        prevOut[0] = (lpRowInitial[0] * cprev);
        prevOut[1] = (lpRowInitial[1] * cprev);
        prevOut[2] = (lpRowInitial[2] * cprev);

        for (int y = 0; y < height; y++) {
            prevOut[0] = ((lpRowInitial[0] * a0a1) - (prevOut[0] * b1b2));
            prevOut[1] = ((lpRowInitial[1] * a0a1) - (prevOut[1] * b1b2));
            prevOut[2] = ((lpRowInitial[2] * a0a1) - (prevOut[2] * b1b2));
            bufferPerLine[0] = prevOut[0];
            bufferPerLine[1] = prevOut[1];
            bufferPerLine[2] = prevOut[2];
            bufferPerLine += Channels;
            lpRowInitial += Channels;
        }
        lpRowInitial -= Channels;
        bufferPerLine -= Channels;
        lpColInitial += WidthStep * HeightSubOne;
        prevOut[0] = (lpRowInitial[0] * cnext);
        prevOut[1] = (lpRowInitial[1] * cnext);
        prevOut[2] = (lpRowInitial[2] * cnext);
        for (int y = HeightSubOne; y >= 0; y--) {
            prevOut[0] = ((lpRowInitial[0] * a2a3) - (prevOut[0] * b1b2));
            prevOut[1] = ((lpRowInitial[1] * a2a3) - (prevOut[1] * b1b2));
            prevOut[2] = ((lpRowInitial[2] * a2a3) - (prevOut[2] * b1b2));
            bufferPerLine[0] += prevOut[0];
            bufferPerLine[1] += prevOut[1];
            bufferPerLine[2] += prevOut[2];
            lpColInitial[0] = bufferPerLine[0];
            lpColInitial[1] = bufferPerLine[1];
            lpColInitial[2] = bufferPerLine[2];
            lpRowInitial -= Channels;
            lpColInitial -= WidthStep;
            bufferPerLine -= Channels;
        }
    }
    else if (Channels == 4)
    {
        float prevOut[4];

        prevOut[0] = (lpRowInitial[0] * cprev);
        prevOut[1] = (lpRowInitial[1] * cprev);
        prevOut[2] = (lpRowInitial[2] * cprev);
        prevOut[3] = (lpRowInitial[3] * cprev);

        for (int y = 0; y < height; y++) {
            prevOut[0] = ((lpRowInitial[0] * a0a1) - (prevOut[0] * b1b2));
            prevOut[1] = ((lpRowInitial[1] * a0a1) - (prevOut[1] * b1b2));
            prevOut[2] = ((lpRowInitial[2] * a0a1) - (prevOut[2] * b1b2));
            prevOut[3] = ((lpRowInitial[3] * a0a1) - (prevOut[3] * b1b2));
            bufferPerLine[0] = prevOut[0];
            bufferPerLine[1] = prevOut[1];
            bufferPerLine[2] = prevOut[2];
            bufferPerLine[3] = prevOut[3];
            bufferPerLine += Channels;
            lpRowInitial += Channels;
        }
        lpRowInitial -= Channels;
        bufferPerLine -= Channels;
        lpColInitial += WidthStep*HeightSubOne;
        prevOut[0] = (lpRowInitial[0] * cnext);
        prevOut[1] = (lpRowInitial[1] * cnext);
        prevOut[2] = (lpRowInitial[2] * cnext);
        prevOut[3] = (lpRowInitial[3] * cnext);
        for (int y = HeightSubOne; y >= 0; y--) {
            prevOut[0] = ((lpRowInitial[0] * a2a3) - (prevOut[0] * b1b2));
            prevOut[1] = ((lpRowInitial[1] * a2a3) - (prevOut[1] * b1b2));
            prevOut[2] = ((lpRowInitial[2] * a2a3) - (prevOut[2] * b1b2));
            prevOut[3] = ((lpRowInitial[3] * a2a3) - (prevOut[3] * b1b2));
            bufferPerLine[0] += prevOut[0];
            bufferPerLine[1] += prevOut[1];
            bufferPerLine[2] += prevOut[2];
            bufferPerLine[3] += prevOut[3];
            lpColInitial[0] = bufferPerLine[0];
            lpColInitial[1] = bufferPerLine[1];
            lpColInitial[2] = bufferPerLine[2];
            lpColInitial[3] = bufferPerLine[3];
            lpRowInitial -= Channels;
            lpColInitial -= WidthStep;
            bufferPerLine -= Channels;
        }
    }
    else if (Channels == 1)
    {
        float prevOut = 0;
        prevOut = (lpRowInitial[0] * cprev);
        for (int y = 0; y < height; y++) {
            prevOut = ((lpRowInitial[0] * a0a1) - (prevOut * b1b2));
            bufferPerLine[0] = prevOut;
            bufferPerLine += Channels;
            lpRowInitial += Channels;
        }
        lpRowInitial -= Channels;
        bufferPerLine -= Channels;
        lpColInitial += WidthStep*HeightSubOne;
        prevOut = (lpRowInitial[0] * cnext);
        for (int y = HeightSubOne; y >= 0; y--) {
            prevOut = ((lpRowInitial[0] * a2a3) - (prevOut * b1b2));
            bufferPerLine[0] += prevOut;
            lpColInitial[0] = bufferPerLine[0];
            lpRowInitial -= Channels;
            lpColInitial -= WidthStep;
            bufferPerLine -= Channels;
        }
    }
}
//本人博客:http://tntmonks.cnblogs.com/ 转载请注明出处.
void  GaussianBlurFilter(unsigned char * input, unsigned char * output, int Width, int Height, int Stride, float GaussianSigma) {

    int Channels = Stride / Width;
    float a0, a1, a2, a3, b1, b2, cprev, cnext;

    CalGaussianCoeff(GaussianSigma, &a0, &a1, &a2, &a3, &b1, &b2, &cprev, &cnext);

    float a0a1 = (a0 + a1);
    float a2a3 = (a2 + a3);
    float b1b2 = (b1 + b2); 

    int bufferSizePerThread = (Width > Height ? Width : Height) * Channels;
    unsigned char * bufferPerLine = (unsigned char*)malloc(bufferSizePerThread);
    unsigned char * tempData = (unsigned char*)malloc(Height * Stride);
    if (bufferPerLine == NULL || tempData == NULL)
    {
        if (tempData)
        {
            free(tempData);
        }
        if (bufferPerLine)
        {
            free(bufferPerLine);
        }
        return;
    }
    for (int y = 0; y < Height; ++y) {
        unsigned char * lpRowInitial = input + Stride * y;
        unsigned char * lpColInitial = tempData + y * Channels;
        gaussianHorizontal(bufferPerLine, lpRowInitial, lpColInitial, Width, Height, Channels, Width, a0a1, a2a3, b1b2, cprev, cnext);
    }
    int HeightStep = Height*Channels;
    for (int x = 0; x < Width; ++x) {
        unsigned char * lpColInitial = output + x*Channels;
        unsigned char * lpRowInitial = tempData + HeightStep * x;
        gaussianVertical(bufferPerLine, lpRowInitial, lpColInitial, Height, Width, Channels, a0a1, a2a3, b1b2, cprev, cnext);
    }

    free(bufferPerLine);
    free(tempData);
}

4.DCT (Discrete Cosine Transform)

原稿链接

邮箱地址是:
gaozhihan@vip.qq.com

7.Deriche (Recursive filter)

背景知识和题材陈述

突发性,对循环进行矢量化处理并不足以提升算法的特性。 英特尔® C++
编译器可能会见提醒“可以矢量化但效率可能坏小”。
但仅仅因为循环可以进行矢量化,并无表示生成的代码比无开展矢量化的巡回再迅速。
如果矢量化无法晋升性,您可选择查明其偷的原故。
通常,为了博便捷的 SIMD 代码,要求再设计数据布局与算法。
许多气象下,无论是否开展矢量化,有利于 SIMD
的优化措施是促成性能提升的机要缘由。 不过,通过升级算法的效率,SIMD
性能以见面明显加强。

本文将介绍示例源代码及其四只其他版本的大循环,说明我们也增长 SIMD
效率所作出的转。 图 1 用故作本文参考和已下载的源代码。 有关版本 0 –
3 的部分是本文的主干部分。 额外的版 4 部分用介绍能够消除 SIMD
转化开销的高等级 SDLT 特性。

图 1: 本子编号图例以及可用源代码中关于代码集修改的相应描述。
版本编号相同含有修改顺序。

渴求数搜集和分散的算法会对标量和 SIMD 的习性造成影响。
而且,如果你拥有募集(分散)链,会更为回落性能。
如果循环中寓间接访问(或非单位步长内存访问4),如图
2
所示,编译器可能会见变收集指令(显式收集指令或多独拟数据收集之通令)。
而且由于用间接访问大型指令,收集指令的多少将趁着数据基元数量之加而上指数级增长。
例如,如果指令“A”包含 4 次加倍,间接访问这同指令以会晤生成 4 单采访指令。
有几情况下,算法中见面不可避免地冒出间接访问。
不过,如发或,您应该找来解决方式避免间接访问。
避免出现收集(或分散)等无效行为,可明白提升 SIMD 性能。

此外,数据对齐也可升级 SIMD 性能。 如果循环在匪对齐 SIMD
数据通道的数目上运行,性能会所有下滑。

图 2
间接内存寻址可能是采访或者分散行为,而且支持循环索引用于找另一样寻觅引。
收集是索引负载。 分散是索引保存。

我们透过一个简单易行的 3D
网格变形算法示例说明几码可用来升级生成代码效率的艺,从而也标量和 SIMD
提供优势。 在图 3 中,3D
网格的每个终端都生一个附件,其中富含可影响顶点变形的数量。
每个附件间接引用 4 独接合点。 附件和接合点保存于 1D 阵列中。

图 3: 3D 网格变形示例算法。

10.FA (Fast Anisotropic)

3.Vliet-Young-Verbeek (Recursive filter)

本子 2(第 2 组成部分): 数据填

此刻,Opt 报告展望性能将大幅升级。
然而,我们以前期的巨型附件循环重新整理成多独小型的子循环,并且我们注意到,在骨子里循环执行进程遭到,处理短行程计数时,性能没达到最佳状态。
就短行程计数而言,未完全矢量化的脱离循环或剩余循环中见面吃大量执时。
图 11
所提供的示范说明非对旅多少会导致在离循环、主循环或剩余循环中推行。
如果迭代空中的起始索引或终止索引(或双方)不是 SIMD
矢量通道数量之翻番,这种情景就算会并发。
理想状态下,我们想保有执行时都起于主 SIMD 循环中。

图 11: SIMD 循环解析。 当编译器执行矢量化时,会也 3 类循环(主 SIMD
循环、剥离循环和多余循环)生成代码。 本图中出一个 4
矢量通道示例,其中循环迭代空间限制也 3-18。 主循环将一次性处理 4
单因素(从 SIMD 通道边界 4 开始,到 15 为止),剥离循环将处理元素
3,剩余循环将拍卖 16-18.

图 12: 英特尔® VTune™ Amplifier XE (2016)
可用于查看相应汇编代码中日耗情况。 查看英特尔 VTune Amplifier XE
中之(已施行)汇编,及其滚动条时,蓝色横条表示执行时间。
通过辨认汇编中之退循环、主循环和剩余循环,可以确定矢量化主循环外部消耗了多长时间(如发)。

之所以,除了整理附件之外,填充附件数据因要该变成 SIMD
矢量通道数量的倍数,也不过升级 SIMD 性能。 图 13
解释了填数据阵列如何支撑所有执行都发生在主 SIMD
循环中,这是一样栽名特新优精状态。
结果也许每发生距离,但填数据一般是同等项大实用的技艺。

图 13: 填充数据阵列。 在 4
矢量通道示例中,显示了附件经过整理,可分为两组子循环。 (左侧)子循环 1
遭,附件 0-3 在主循环中处理,剩下零星单因素(4 同 5)由剩余循环处理。
子循环 2 中,仅来一个吧 3 的路计数,全部 3 只均是因为剥离循环处理。
(右侧)我们填充了所有子循环,以要其成 4 长达 SIMD
通道的翻番,从而支持具有附件均是因为矢量化循环来处理。

8.ebox (Extended Box)

http://dx.doi.org/10.1109/ICPR.1998.711192

本子 2(第 1 有些): 通过整理对数码进行预处理为管教数量统一

于咱们的言传身教中,Opt-报告报告了 48 次收集及呼应的“间接访问”备注。
我们要专门关爱间接访问备注,因为拖欠报告被几乎统统是这种备注。
经过进一步考察,我们发现其分别针对应 4
个(矢量化循环中间接访问的)接合点的 4×3 只矩阵值,总共 48 次收集。
我们明白收集(或分散)会影响属性。 那咱们该如何化解就同一问题吧。
必须开展募集,还是应想方法避免这些收集行为?

例如,我们问自己“是否来能够打循环体内部访问,并会晋级及循环外部的联合数?”
最初的答案是“否”,然后我们问“是否能够又规划算法,以便生成为循环不变式的统一数?”

图 8: 整理算法数据。 左侧,循环迭代接合点索引完全两样之附件。
右侧,附件经过整理,每个(内层)子循环都发出同样的合并衔接合点数据集。

一经图 8 所示,许多独门附件都只是共享相同之接合点索引值。
通过整治附件,所有共享相同索引的附件可分割以同样组,这样用产生机会循环附件的子集,其中接合点均匀(循环不变式)分布于子循环的迭代空间。
这样将力促将接合点数据提升至矢量化内层循环的外部。
这样内层矢量化循环将无见面冒出其他采访行为。

图 9: 版本 2: 重新设计算法为创办统一(循环不变式)数据。

祈求 9
显示了应用整理数据阵列后转变的代码,它以共享统一数的因素分以平等组,其中原始循环转化成可避免收集行为之外层和内层(矢量化)循环。
IndiceSet
的 mIndiceSetArray 阵列能够跟整理后阵列中之起始索引和停止索引。
因此我们来一个外层循环和一个内层循环。
而且,由于数量经过又排序,所以需要续加 workIdIndex 以钉原始位置,以便写起结果。

今日 Opt 报告(见图 10)不再报告接合点造成的 48
次索引掩码负载(或收集行为)。 而且“预计”英特尔® AVX 性能拿提升 2.35
倍增
。 在我们的案例中,实际性能提升了 2.30 倍

图 10: 版本 2: 关于通过集合衔接合点数据再度设计后底大循环的英特尔® C++
编译器 Opt-报告。

于图 10 中,我们应留神,Opt-报告仍报告了 8 次“收集”或“掩码步长负载”。
造成这种结果的原由是访问 mSortedAttatchments 阵列的结构阵列内存布局。
理想图景下,我们期望实现“非掩码对同单元步长”负载。
稍后我们以证明如何借助 SDLT 实现就无异于对象。
同样值得注意的凡,Opt-报告(见图 10)中起了 3 次分散。
这是以我们重新排列了输入数据,因此待盖科学的逐一以结果写有至输出(如图
9 第 29 行所示)。 但分散 3 个价如果好于收集 48
只价,因为光待少量开发就只是免大量资产。

http://dx.doi.org/10.1007/978-3-642-24785-9\_38

http://dx.doi.org/10.1109/78.295213

简介

为发挥
SIMD1 的无限老作用,除了对那开展矢量化处理2外,我们尚需要作出任何努力。可以品味吧循环添加 #pragma
omp
simd
3,查看编译器是否成功进行矢量化,如果性能兼备升级,则达到满意状态。
然而,可能性能根本无见面提升,甚至还会降。
无论处在何种情况,为了最特别限度发挥 SIMD
执行的优势并促成性能提升,通常需要还规划算法和多少布局,以便生成的
SIMD 代码尽可能快。
另外还可收到额外的职能,即标量(非矢量化)版代码会表现得重好。

正文将由此一个 3D
动画算法示例,逐步介绍除添加“#pragma”之外我们尚只是下什么方法。
在当时同样历程遭到,其他技术及道也可是也汝的生一样步矢量化工作提供帮扶。
我们还合并了算法和 SIMD 数据布局模板 (SDLT) — 英特尔® C++
编译器的一律宗特征,以增强多少布局和 SIMD 的效率。
本文的具备源代码均只是下载,并带有此处未提及的其他详细信息。

  注:支持通道数分别吗 1 ,3 ,4.

版本 1: SIMD

至于本 1,我们针对循环进行矢量化处理。 在演示中,通过抬高“#pragma omp
simd
”即可成功而循环实现矢量化(见图
5),因为其满足矢量化标准(例如,没有函数调用,单上前只有,直线式代码5)。
另外,它本 SDLT
的矢量化策略,即限制对象为拉编译器成功就私有化。6唯独,我们承诺注意,在广大景下,简单加加编译指示会招编译错误或错非常成代码。7寻常需进行代码重构,以要循环达到可矢量化状态。

图 5: 版本 1: 修改版 0的第 8 行(见图
4)以对循环进行矢量化处理。

希冀 6 显示了英特尔® C++ 编译器 (ICC) XE 关于本 1 循环的 Opt-报告。
对英特尔® 高级矢量扩展(英特尔®
AVX)9构建而言,大家可以看,Opt-报告显示就循环实现了矢量化,性能预计只能提升
5%。 但在我们的案例被,版本 1 的莫过于性能比版本 0 低 15%。 无论
Opt-报告展望大多异常程度之属性提升,都当进行测试为了解实际性能。

此外,图 6 显示了由全 4 只 Joint 构成的变矩阵中,每次加倍都起 48
项“间接访问”掩码索引负载。 同时相应地挺成了 48 独“间接访问”备注,图 7
列举了里的一个。 我们无应该忽视
Opt-报告备注;而应查明由并尝试解决这些题目。

图 6: 版本 1: 关于循环的英特尔® C++ 编译器 Opt-报告。

图 7: 版本 1: 英特尔® C++ 编译器 Opt-报告,间接访问备注。

即便循环实现了矢量化,间接访问造成的大量征集行为还会拦 SIMD
提升性。

  注:支持通道数分别吗 1 ,3 ,4.

5.box (Box filter)

版本 0: 算法

在图 4 中,算法循环迭代“附件”阵列。 每个附件包含 4
个在“接合点”阵列中间接拜的接合点索引值。 而且每个接合点包含一个是因为 12
次加倍组成的易矩阵 (3×4)。 因此每次循环迭代要求采访 48 次倍(12
次加倍乘以 4 个接合点)。 收集如此多加倍会降低 SIMD 性能。
因此,如果缩减或者避免这种集,SIMD 性能将会晤明确提升。

图 4: 版本 0: 每次循环迭代包含 48 次收集的演示算法。

终极以本人2.20GHz底CPU上,单对单线程,不动流动(SIMD)指令,达到了,处理一千六百万诸如从的彩色照片仅用700毫秒左右.

http://www.jstor.org/stable/2158018

版本 4:sdlt::uniform_soa_over1d

每当本子 4 的算法中,我们尚发现了多其它的性提升会。
请注意,从一个子循环到下一个子循环,4 单接合点数据遭到的 3
只都只是提出供内层循环中集合看。 而且大家还答应询问,为每次上 SIMD
循环准备联合数会来比充分之开发,而且我们以针对统一数的历次外层循环迭代过程遭到还设肩负这种资本。

即便 SIMD 循环而言,根据 SIMD
指令集的不比,10迭代初始前准备合并数会发出。
对每个统一值来说,编译器可以 1) 将标量加载到寄存器,2)
将寄存器中之标量值传播到 SIMD 寄存器的有通道,然后 3) 将 SIMD
寄存器保存在堆栈上的新职务,以供 SIMD 循环体使用。
就长行程计数而言,这同样开支可轻松摊销。
但就短行程计数而言,它或许会见潜移默化属性。 在本子 3
中,每次外层循环迭代都见面发出关于 4 个接合点的开,总共 48
次倍(每个接合点 12 次加倍)。

图 16: 探寻行程计数: 英特尔® Advisor XE (2016)
具有相同件能够为循环执行提供行程计数的实用功能。
它可助大家轻松识别短行程计数和长行程计数。

于这种气象下,SDLT 可支持通过规定何时来出,明确管理这种 SIMD
数据转发,而非会见自行出出成本。 SDLT
的这种高级功能为 sdlt::uniform_soa_over1d。 它促进排除循环中之
SIMD 转化开销,并扶持用户控制来出的光阴。 其规律是为 SIMD
就绪型格式保存循环不变式数据,以便 SIMD 循环不经过转化就可一直看数。
它支持统一数实现有更新与重复用,从而帮助提升性。

为证明 SIMD 数据转发开销何时起,以及 SDLT
如何帮助排除这等同开发,我们当图 17 和 18 中提供了一个伪代码示例。 图 17
显示每次外层循环迭代(第 18 行)和每次加倍(通过集合数看)(12
次加倍)都见面生出。 图 18
显示了 sdlt::uniform_soa_over1d 的用,它能够如这种支付就来同样蹩脚(第
6 行),从而降低一体化资产。 这种高档功能会吗特定情景带来巨大优势。
用户应进行考查。 实际结果也许截然不同。

图 17: 进入第 12 行的 SIMD 循环之前,对每个统一值来说,编译器可以
1) 将标量加载到寄存器,2) 将寄存器中的标量值传播及 SIMD
寄存器的备通道,然后 3) 将 SIMD 寄存器保存在堆栈上的新岗位,以供 SIMD
循环体使用。 就长行程计数而言,这同支出可轻松摊销。
但就短行程计数而言,它恐怕会见潜移默化性。

图 18: 通过应用
sdlt::uniform_soa_over1d,确定何时有出,可明明管理这种 SIMD
数据转发,而非会见自行出成本。 SDLT 的这项高级功能推进清除循环中之
SIMD 转化开销,并拉扯而决定来出的年华。 其规律是坐 SIMD
就绪型格式保存循环不变式数据,以便 SIMD 循环不经过转化就不过径直看数。

因此首先步是在短行程计数情况下更是升级性,我们再次规划算法,以便更利用自外围循环索引 i至 i+1 的
4 个接合点中的老三个,如图 19 所显示。 使用 SDLT 功能而帮忙解除为子循环准备
SIMD 数据时所共的开。

图 19: 版本 3: 面向里 3 个接合点(共 4
个)的联合数只是重新用于下一子循环
可有的更新统一数,以最好老限度地减少负荷(或外层循环中的搜集)。
而且使用 sdlt::uniform_soa_over1d 以 SIMD
就绪型格式保存统一数,还只是尽老限度地下降涉及所有子循环的 SIMD
转化开销。

由此重构循环以便将合并数更利用于历子循环,而且只有需要进行局部更新,平均来说,我们特需要创新
1/4 的统一数。 因此,设置用于 SIMD 循环的集合数常常所生的开支将跌
75%。

故基于英特尔即时卖代码,俺对该展开了改写以及优化.

余写算法追求干净整洁,高效简明,换言之即是休使其他硬件加速方案,实现简单便捷,以适应不同硬件环境.

 

 

一体化代码:

脚注

1 单指令多数据 (SIMD) 指利用数据层并行化,单条指令以处理多单数据。
它同习俗“标量操作”(使用单条指令处理单个数据)正好相反。

2 矢量化指计算机程序于标量实施转化为矢量(或 SIMD)实施。

pragma
simd
: https://software.intel.com/zh-cn/node/583427。 pragma
omp
simd
: https://software.intel.com/zh-cn/node/583456

4 非单位步长内存访问指以循环连续增量的进程被,从内存的非相邻位置访问数。
这样见面严重影响性。
相反,以单位步长(或相继)形式拜访内存会显著提高效率。

5 以供应参考: https://software.intel.com/sites/default/files/8c/a9/CompilerAutovectorizationGuide.pdf

6 SDLT 基元可限对象,有助于编译器成功做到 SIMD
循环中针对有变量的私有化,即每个 SIMD 通道还有一个私有变量实例。
为满足当下无异业内,对象要是 Plain Old Data
(POD),拥有行内对象成员,没有嵌套阵列,而且尚未虚构函数。

7 在矢量化过程中,开发人员应该试验各种不同之编译指示(比如 simd、omp
simd、ivdep 和矢量 {always [assert]}),并使用 Opt-报告。

8 我们啊基于 Linux* 的英特尔® C++ 编译器 16.0 (2016)
添加了命行选项“-qopt-report=5 –qopt-report-phase=vec”,以稀成 Opt-报告
(*.optrpt)。

9 以英特尔® C++ 编译器 16.0 生成英特尔® 高级矢量扩展(英特尔® AVX)
指令时,可拿选项 “-xAVX” 添加至编译命令行。

10AVX512 指令集带有传播负载指令,可落迭代上马之前准备联合数时所发出的
SIMD 开销。

i 以性检测过程被关系的软件及其性质只有当英特尔微处理器的架构下方能博得优化。

诸如 SYSmark* 和 MobileMark*
等测试都系因特定计算机体系、硬件、软件、操作系统和作用。
上述任何因素的更改都发出或导致测试结果的变。
请参考其他消息及性测试(包括结合其他产品使用时之运行性能)以对目标产品进行宏观评估。
更多信息要访问 http://www.intel.com/content/www/cn/zh/benchmarks/intel-product-performance.html。

布置: 英特尔® 至强® 处理器 E5-2699 v3(45M 高速缓存,2.30 GHz). CPU:
两个 18 核 C1,2.3 GHz。 非根本: 2.3 GHz. 英特尔® 快速通道互联技术: 9.6
GT/秒。 RAM: 128 GB,DDR4-2133 MHz (16 x 8 GB)。 磁盘: 7200 RPM SATA
磁盘。 800 GB 固态盘。 关闭超线程技术,关闭睿频加速。 Red Hat Enterprise
Linux* Server 7.0 版。 3.10.0-123.el7.x86_64。

 

思了相思,还是以代码共享出来,供大家参考学习.

事先,俺呢发过不少快速高斯模糊算法.

解决方案

矢量化成功后,性能可能会见有着升级,也恐怕不见面。
无论何种状况,对循环进行矢量化处理还应只是是优化过程的起点,而不是归根到底点。
相反,我们得使用工具(例如,Opt-报告、汇编代码,英特尔® VTune™
Amplifier XE、英特尔® Advisor
XE)帮助调查效率低下的原委,并经过实施迎刃而解方案来改善 SIMD 代码。

  GaussianBlurFilter(输入图像数据,输出图像数据,宽度,高度,通道数,强度)

IIR Gaussian Blur Filter Implementation using Intel® Advanced Vector
Extensions [PDF
513KB]
source: gaussian_blur.cpp [36KB]

结论

图 20: 基于英特尔® 交强® CPU 处理器 E5-2699 v3 (代号
Haswell)构建的英特尔® 高级矢量扩展所实现的特性提升i

对代码进行矢量化处理只是实现 SIMD 性能加速的初步,而休是最终。
矢量化之后该采取可用资源同工具(例如优化报告、英特尔 VTune Amplifier
XE、英特尔 Advisor XE)了解就生成代码的效率。
经过分析后,可以发现对标量和 SIMD 代码有用之优化会。
然后动并测试各项技能,无论是大技巧还是本文提到的技术。
你或许用重新考虑算法和数码布局,以极其可怜限度地增进代码(尤其是已经成形的汇编代码)的效率。

通过演示我们得以看到,对算法实施数据预处理为解所有间接访问(收集)的版
2 收到了无限充分作用。 版本 3 使用
SDLT,通过非掩码对一头单位步长负载改进内存访问,并填写充数据为对齐 SIMD
通道边界,显著升级了性能。 在短行程计数场景中,我们下 SDLT
高级功能,最要命限度地回落了合数开销的完全资产。

俺目前只是把opencv当资料库来拘禁,并不认为opencv可以用于绝大多数底生意项目.

题外话:

 

5.box (Box filter)

使了英特尔计算机的流(SIMD)指令,算法处理速度极其惊人.

版本 3: SDLT 容器

出于我们重新设计了算法为避免收集行为并判升级 SIMD
性能,现在我们得采用 SDLT 进一步提高 SIMD 代码的频率。
到现毕,所有负载都也“掩码”,且无对伙同。
理想状态下,我们希望负载为未掩码、对同、单位步长负载。 我们下 SDLT
基元和容器来兑现即时同对象。 SDLT 有助于使 SIMD
循环中的片变量成功完成私有化,即每个 SIMD 通道还起一个民用变量实例。
SDLT 容器和访问器将机关处理数量易与指向旅。

图 14 中,源代码显示了集成 SDLT 所待的改动。
主要的变更是为指令 AttachmentSorted 声明SDLT_PRIMITIVE,然后用用来附件阵列的输入数据容器从 std::vector 容器(结构阵列
(AOS) 数据布局)转化成 SDLT 容器。 编程人员可下 SDLT
访问器上之演算符 [],就比如它是 C 阵列或 std::vector。 最初我们用
SDLT 阵列结构 (SOA) 容器 (sdlt::soa1d_container),但阵列结构阵列
(ASA) 容器 (sdlt::asa1d_container) 也只是眼看升级性。 转换(使用
typedef)SDLT
容器类型因测试是否获得最佳性能的法非常简单,因此建议大家利用。 图 14
中,我们引入了 SDLT_SIMD_LOOP 宏指令,即 ICC 16.2 (SDLT v2)
中之“预览”特性,且兼容 ASA 和 SOA 两种植容器类。

图 14: 版本 3。 集成 SDLT 容器(第 1–3 行,以及第 7
行)和访问器(第 8 和 19 行);同样应用 SDLT_SIMD_LOOP
宏指令的“预览”特性(第 17 和 23 行)。 仅显示版本 2 的不同之处。

图 15: 版本 3: 使用 SDLT 基元和容器时的英特尔® C++ 编译器
Opt-报告。

图 15 中,Opt-报告预测本 3 将促成 1.88 倍的属性提升。
但请记住,这才是预估值,而不是实在性能提升。
事实上,我们的案例实现了 3.17 倍的其实性能提升。
另外,回忆一下,版本 2 的 Opt 报告(图 10)报告了“掩码步长”负载。
现在(图 15)的载重为“非掩码”、“对共同”、“单位步长”。
这是极其出彩的特性状态,可经过应用 SDLT
容器改进数据布局并加强内存访问效率来落实。

前面为闹网友咨询了这算法的实现问题.

2.SII (Stacked integral images)

参考资料

  1. 下载包含源代码的言传身教:
    https://software.intel.com/sites/default/files/managed/de/3f/animation-simd-sdlt-whitepaper.tar.gz
  2. SDLT 文档(包含有代码示例):
    https://software.intel.com/zh-cn/code-samples/intel-compiler/intel-compiler-features/intel-sdlt
  3. SIGGRAPH 2015: DreamWorks 动画 (DWA): 如何依靠 SIMD
    将肌肤变形性能提升 4 加倍:
    http://www.slideshare.net/IntelSoftware/dreamwork-animation-dwa
  4. “先试后买”英特尔® Parallel Studio XE 评估版:
    http://software.intel.com/intel-parallel-studio-xe/try-buy
  5. 面向合格学生、教育工作者、学术研究人员与开源工作者的英特尔® Parallel
    Studio XE 免费版本:
    https://software.intel.com/zh-cn/qualify-for-free-software
  6. 英特尔® VTune™ Amplifier 2016:
    https://software.intel.com/zh-cn/intel-vtune-amplifier-xe
  7. 英特尔® Advisor:
    https://software.intel.com/zh-cn/intel-advisor-xe

要要一律步一个脚印去实现部分最好核心的算法,扎实的根底才是构建上层建筑的中心条件.

世武功,唯快不清除。
本文只是抛砖引玉一下,若有其它有关题材或者需要吗得邮件联系个人探讨。

2.SII (Stacked integral images)

前为发生网友咨询了是算法的实现问题.

4.DCT (Discrete Cosine Transform)

有关IIR相关知识,参阅 百过词条 “IIR数字滤波器”

http://dx.doi.org/10.1109/ICPR.1998.711192

 

http://arxiv.org/abs/1107.4958

要么如一如既往步一个脚印错过落实部分最中心的算法,扎实的基本功才是构建上层建筑的中心条件.

https://software.intel.com/zh-cn/articles/iir-gaussian-blur-filter-implementation-using-intel-advanced-vector-extensions

http://dx.doi.org/10.1109/TPAMI.1986.4767776

过多网友直强调使用opencv,opencv的确十分强劲,但是若想要起再度老的迈入空间以及开创力.

末了以余2.20GHz的CPU上,单对单线程,不采用流动(SIMD)指令,达到了,处理一千六百万诸如从的彩色照片仅用700毫秒左右.

运了英特尔计算机的流(SIMD)指令,算法处理速度极其惊人.

http://arxiv.org/abs/1107.4958

 

  GaussianBlurFilter(输入图像数据,输出图像数据,宽度,高度,通道数,强度)

http://mathinfo.univ-reims.fr/IMG/pdf/Fast\_Anisotropic\_Gquss\_Filtering\_-\_GeusebroekECCV02.pdf

随常规,还是贴个效果图于直观.

事先犯的几乎单算法,在本人2.2GHz的CPU上耗时还见面过1秒.

我写算法追求干净卫生,高效简明,换言之便是无采用任何硬件加速方案,实现简单快速,以适应不同硬件环境.

即人家经过实测,IIR是兼职效果以及性能的对的计,也是半径无关(即模糊不同强度耗时基本不变换)的实现.

至于IIR相关知识,参阅 百渡过词条 “IIR数字滤波器”

6.AM(Alvarez, Mazorra)

脚下本人经过实测,IIR是兼职效果以及性能的是的法门,也是半径无关(即模糊不同强度耗时基本未转移)的实现.

前面犯之几只算法,在身2.2GHz的CPU上耗时还见面跨1秒.

倘鲜明,快速高斯模糊有很多贯彻方式:

1.FIR (Finite impulse response)

调用方法:

http://dx.doi.org/10.1016/0165-1684(95)00020-E

余一般认为,只要处理一千六百万如从彩色图片,在2.2GHz的CPU上才对单线程超过1秒的算法,都是难过的.

8.ebox (Extended Box)

http://dx.doi.org/10.1016/0165-1684(95)00020-E

邮箱地址是:
gaozhihan@vip.qq.com

……

http://dx.doi.org/10.1007/978-3-642-24785-9\_38

若本文帮到您,厚颜无耻求微信扫码打个赏.

http://baike.baidu.com/view/3088994.htm

前,俺呢发过不少飞速高斯模糊算法.

若是显,快速高斯模糊有诸多贯彻方式:

贯彻高斯模糊的点子虽然多,但是当算法而言,核心要是简单高效.

俺一般认为,只要处理一千六百万诸如从彩色图片,在2.2GHz的CPU上仅核单线程超过1秒的算法,都是难过的.

http://dx.doi.org/10.1109/78.295213

https://zh.wikipedia.org/wiki/%E9%AB%98%E6%96%AF%E6%A8%A1%E7%B3%8A

http://mathinfo.univ-reims.fr/IMG/pdf/Fast\_Anisotropic\_Gquss\_Filtering\_-\_GeusebroekECCV02.pdf

据老,还是贴个效果图于直观.

10.FA (Fast Anisotropic)

https://zh.wikipedia.org/wiki/%E9%AB%98%E6%96%AF%E6%A8%A1%E7%B3%8A

 

花特尔官方实现之即刻卖:

无数网友直接强调使用opencv,opencv的确挺强劲,但是倘若想要生还老的腾飞空间以及开创力.

http://dx.doi.org/10.1109/TPAMI.1986.4767776

个人目前只是将opencv当资料库来拘禁,并不认为opencv可以用来绝大多数底商业项目.

1.FIR (Finite impulse response)

若本文帮到您,厚颜无耻求微信扫码打个赏.

http://www.jstor.org/stable/2158018

9.IIR (Infinite Impulse Response)

6.AM(Alvarez, Mazorra)

纪念了纪念,还是以代码共享出来,供大家参考学习.

IIR Gaussian Blur Filter Implementation using Intel® Advanced Vector
Extensions [PDF
513KB]
source: gaussian_blur.cpp [36KB]

http://hal.inria.fr/docs/00/07/47/78/PDF/RR-1893.pdf

7.Deriche (Recursive filter)

完全代码:

……

http://dx.doi.org/10.1109/ROBOT.2010.5509400

https://software.intel.com/zh-cn/articles/iir-gaussian-blur-filter-implementation-using-intel-advanced-vector-extensions

大地武功,唯快不脱。
正文只是抛砖引玉一下,若发生另有关题材要么要求也可邮件联系我探讨。

实现高斯模糊的措施虽然多,但是当算法而言,核心重点是概括高效.

http://hal.inria.fr/docs/00/07/47/78/PDF/RR-1893.pdf

 

 

故基于英特尔即时卖代码,俺对该展开了改写以及优化.

9.IIR (Infinite Impulse Response)

http://dx.doi.org/10.1109/ROBOT.2010.5509400

http://baike.baidu.com/view/3088994.htm

花特尔官方实现之即卖:

3.Vliet-Young-Verbeek (Recursive filter)

void CalGaussianCoeff(float sigma, float * a0, float * a1, float * a2, float * a3, float * b1, float * b2, float * cprev, float * cnext) {
    float alpha, lamma, k;

    if (sigma < 0.5f)
        sigma = 0.5f;
    alpha = (float)exp((0.726) * (0.726)) / sigma;
    lamma = (float)exp(-alpha);
    *b2 = (float)exp(-2 * alpha);
    k = (1 - lamma) * (1 - lamma) / (1 + 2 * alpha * lamma - (*b2));
    *a0 = k; *a1 = k * (alpha - 1) * lamma;
    *a2 = k * (alpha + 1) * lamma;
    *a3 = -k * (*b2);
    *b1 = -2 * lamma;
    *cprev = (*a0 + *a1) / (1 + *b1 + *b2);
    *cnext = (*a2 + *a3) / (1 + *b1 + *b2);
}

void gaussianHorizontal(unsigned char * bufferPerLine, unsigned char * lpRowInitial, unsigned char  * lpColumn, int width, int height, int Channels, int Nwidth, float a0a1, float a2a3, float b1b2, float  cprev, float cnext)
{
    int HeightStep = Channels*height;
    int WidthSubOne = width - 1;
    if (Channels == 3)
    {
        float prevOut[3];
        prevOut[0] = (lpRowInitial[0] * cprev);
        prevOut[1] = (lpRowInitial[1] * cprev);
        prevOut[2] = (lpRowInitial[2] * cprev);
        for (int x = 0; x < width; ++x) {
            prevOut[0] = ((lpRowInitial[0] * (a0a1)) - (prevOut[0] * (b1b2)));
            prevOut[1] = ((lpRowInitial[1] * (a0a1)) - (prevOut[1] * (b1b2)));
            prevOut[2] = ((lpRowInitial[2] * (a0a1)) - (prevOut[2] * (b1b2)));
            bufferPerLine[0] = prevOut[0];
            bufferPerLine[1] = prevOut[1];
            bufferPerLine[2] = prevOut[2];
            bufferPerLine += Channels;
            lpRowInitial += Channels;
        }
        lpRowInitial -= Channels;
        lpColumn += HeightStep * WidthSubOne;
        bufferPerLine -= Channels;
        prevOut[0] = (lpRowInitial[0] * cnext);
        prevOut[1] = (lpRowInitial[1] * cnext);
        prevOut[2] = (lpRowInitial[2] * cnext);

        for (int x = WidthSubOne; x >= 0; --x) {
            prevOut[0] = ((lpRowInitial[0] * (a2a3)) - (prevOut[0] * (b1b2)));
            prevOut[1] = ((lpRowInitial[1] * (a2a3)) - (prevOut[1] * (b1b2)));
            prevOut[2] = ((lpRowInitial[2] * (a2a3)) - (prevOut[2] * (b1b2)));
            bufferPerLine[0] += prevOut[0];
            bufferPerLine[1] += prevOut[1];
            bufferPerLine[2] += prevOut[2];
            lpColumn[0] = bufferPerLine[0];
            lpColumn[1] = bufferPerLine[1];
            lpColumn[2] = bufferPerLine[2];
            lpRowInitial -= Channels;
            lpColumn -= HeightStep;
            bufferPerLine -= Channels;
        }
    }
    else if (Channels == 4)
    {
        float prevOut[4];

        prevOut[0] = (lpRowInitial[0] * cprev);
        prevOut[1] = (lpRowInitial[1] * cprev);
        prevOut[2] = (lpRowInitial[2] * cprev);
        prevOut[3] = (lpRowInitial[3] * cprev);
        for (int x = 0; x < width; ++x) {
            prevOut[0] = ((lpRowInitial[0] * (a0a1)) - (prevOut[0] * (b1b2)));
            prevOut[1] = ((lpRowInitial[1] * (a0a1)) - (prevOut[1] * (b1b2)));
            prevOut[2] = ((lpRowInitial[2] * (a0a1)) - (prevOut[2] * (b1b2)));
            prevOut[3] = ((lpRowInitial[3] * (a0a1)) - (prevOut[3] * (b1b2)));

            bufferPerLine[0] = prevOut[0];
            bufferPerLine[1] = prevOut[1];
            bufferPerLine[2] = prevOut[2];
            bufferPerLine[3] = prevOut[3];
            bufferPerLine += Channels;
            lpRowInitial += Channels;
        }
        lpRowInitial -= Channels;
        lpColumn += HeightStep * WidthSubOne;
        bufferPerLine -= Channels;

        prevOut[0] = (lpRowInitial[0] * cnext);
        prevOut[1] = (lpRowInitial[1] * cnext);
        prevOut[2] = (lpRowInitial[2] * cnext);
        prevOut[3] = (lpRowInitial[3] * cnext);

        for (int x = WidthSubOne; x >= 0; --x) {
            prevOut[0] = ((lpRowInitial[0] * a2a3) - (prevOut[0] * b1b2));
            prevOut[1] = ((lpRowInitial[1] * a2a3) - (prevOut[1] * b1b2));
            prevOut[2] = ((lpRowInitial[2] * a2a3) - (prevOut[2] * b1b2));
            prevOut[3] = ((lpRowInitial[3] * a2a3) - (prevOut[3] * b1b2));
            bufferPerLine[0] += prevOut[0];
            bufferPerLine[1] += prevOut[1];
            bufferPerLine[2] += prevOut[2];
            bufferPerLine[3] += prevOut[3];
            lpColumn[0] = bufferPerLine[0];
            lpColumn[1] = bufferPerLine[1];
            lpColumn[2] = bufferPerLine[2];
            lpColumn[3] = bufferPerLine[3];
            lpRowInitial -= Channels;
            lpColumn -= HeightStep;
            bufferPerLine -= Channels;
        }
    }
    else if (Channels == 1)
    {
        float prevOut = (lpRowInitial[0] * cprev);

        for (int x = 0; x < width; ++x) {
            prevOut = ((lpRowInitial[0] * (a0a1)) - (prevOut  * (b1b2)));
            bufferPerLine[0] = prevOut;
            bufferPerLine += Channels;
            lpRowInitial += Channels;
        }
        lpRowInitial -= Channels;
        lpColumn += HeightStep*WidthSubOne;
        bufferPerLine -= Channels;

        prevOut = (lpRowInitial[0] * cnext);

        for (int x = WidthSubOne; x >= 0; --x) {
            prevOut = ((lpRowInitial[0] * a2a3) - (prevOut  * b1b2));
            bufferPerLine[0] += prevOut;
            lpColumn[0] = bufferPerLine[0];
            lpRowInitial -= Channels;
            lpColumn -= HeightStep;
            bufferPerLine -= Channels;
        }
    }
}

void gaussianVertical(unsigned char * bufferPerLine, unsigned char * lpRowInitial, unsigned char * lpColInitial, int height, int width, int Channels, float a0a1, float a2a3, float b1b2, float  cprev, float  cnext) {

    int WidthStep = Channels*width;
    int HeightSubOne = height - 1;
    if (Channels == 3)
    {
        float prevOut[3];
        prevOut[0] = (lpRowInitial[0] * cprev);
        prevOut[1] = (lpRowInitial[1] * cprev);
        prevOut[2] = (lpRowInitial[2] * cprev);

        for (int y = 0; y < height; y++) {
            prevOut[0] = ((lpRowInitial[0] * a0a1) - (prevOut[0] * b1b2));
            prevOut[1] = ((lpRowInitial[1] * a0a1) - (prevOut[1] * b1b2));
            prevOut[2] = ((lpRowInitial[2] * a0a1) - (prevOut[2] * b1b2));
            bufferPerLine[0] = prevOut[0];
            bufferPerLine[1] = prevOut[1];
            bufferPerLine[2] = prevOut[2];
            bufferPerLine += Channels;
            lpRowInitial += Channels;
        }
        lpRowInitial -= Channels;
        bufferPerLine -= Channels;
        lpColInitial += WidthStep * HeightSubOne;
        prevOut[0] = (lpRowInitial[0] * cnext);
        prevOut[1] = (lpRowInitial[1] * cnext);
        prevOut[2] = (lpRowInitial[2] * cnext);
        for (int y = HeightSubOne; y >= 0; y--) {
            prevOut[0] = ((lpRowInitial[0] * a2a3) - (prevOut[0] * b1b2));
            prevOut[1] = ((lpRowInitial[1] * a2a3) - (prevOut[1] * b1b2));
            prevOut[2] = ((lpRowInitial[2] * a2a3) - (prevOut[2] * b1b2));
            bufferPerLine[0] += prevOut[0];
            bufferPerLine[1] += prevOut[1];
            bufferPerLine[2] += prevOut[2];
            lpColInitial[0] = bufferPerLine[0];
            lpColInitial[1] = bufferPerLine[1];
            lpColInitial[2] = bufferPerLine[2];
            lpRowInitial -= Channels;
            lpColInitial -= WidthStep;
            bufferPerLine -= Channels;
        }
    }
    else if (Channels == 4)
    {
        float prevOut[4];

        prevOut[0] = (lpRowInitial[0] * cprev);
        prevOut[1] = (lpRowInitial[1] * cprev);
        prevOut[2] = (lpRowInitial[2] * cprev);
        prevOut[3] = (lpRowInitial[3] * cprev);

        for (int y = 0; y < height; y++) {
            prevOut[0] = ((lpRowInitial[0] * a0a1) - (prevOut[0] * b1b2));
            prevOut[1] = ((lpRowInitial[1] * a0a1) - (prevOut[1] * b1b2));
            prevOut[2] = ((lpRowInitial[2] * a0a1) - (prevOut[2] * b1b2));
            prevOut[3] = ((lpRowInitial[3] * a0a1) - (prevOut[3] * b1b2));
            bufferPerLine[0] = prevOut[0];
            bufferPerLine[1] = prevOut[1];
            bufferPerLine[2] = prevOut[2];
            bufferPerLine[3] = prevOut[3];
            bufferPerLine += Channels;
            lpRowInitial += Channels;
        }
        lpRowInitial -= Channels;
        bufferPerLine -= Channels;
        lpColInitial += WidthStep*HeightSubOne;
        prevOut[0] = (lpRowInitial[0] * cnext);
        prevOut[1] = (lpRowInitial[1] * cnext);
        prevOut[2] = (lpRowInitial[2] * cnext);
        prevOut[3] = (lpRowInitial[3] * cnext);
        for (int y = HeightSubOne; y >= 0; y--) {
            prevOut[0] = ((lpRowInitial[0] * a2a3) - (prevOut[0] * b1b2));
            prevOut[1] = ((lpRowInitial[1] * a2a3) - (prevOut[1] * b1b2));
            prevOut[2] = ((lpRowInitial[2] * a2a3) - (prevOut[2] * b1b2));
            prevOut[3] = ((lpRowInitial[3] * a2a3) - (prevOut[3] * b1b2));
            bufferPerLine[0] += prevOut[0];
            bufferPerLine[1] += prevOut[1];
            bufferPerLine[2] += prevOut[2];
            bufferPerLine[3] += prevOut[3];
            lpColInitial[0] = bufferPerLine[0];
            lpColInitial[1] = bufferPerLine[1];
            lpColInitial[2] = bufferPerLine[2];
            lpColInitial[3] = bufferPerLine[3];
            lpRowInitial -= Channels;
            lpColInitial -= WidthStep;
            bufferPerLine -= Channels;
        }
    }
    else if (Channels == 1)
    {
        float prevOut = 0;
        prevOut = (lpRowInitial[0] * cprev);
        for (int y = 0; y < height; y++) {
            prevOut = ((lpRowInitial[0] * a0a1) - (prevOut * b1b2));
            bufferPerLine[0] = prevOut;
            bufferPerLine += Channels;
            lpRowInitial += Channels;
        }
        lpRowInitial -= Channels;
        bufferPerLine -= Channels;
        lpColInitial += WidthStep*HeightSubOne;
        prevOut = (lpRowInitial[0] * cnext);
        for (int y = HeightSubOne; y >= 0; y--) {
            prevOut = ((lpRowInitial[0] * a2a3) - (prevOut * b1b2));
            bufferPerLine[0] += prevOut;
            lpColInitial[0] = bufferPerLine[0];
            lpRowInitial -= Channels;
            lpColInitial -= WidthStep;
            bufferPerLine -= Channels;
        }
    }
}
//本人博客:http://tntmonks.cnblogs.com/ 转载请注明出处.
void  GaussianBlurFilter(unsigned char * input, unsigned char * output, int Width, int Height, int Stride, float GaussianSigma) {

    int Channels = Stride / Width;
    float a0, a1, a2, a3, b1, b2, cprev, cnext;

    CalGaussianCoeff(GaussianSigma, &a0, &a1, &a2, &a3, &b1, &b2, &cprev, &cnext);

    float a0a1 = (a0 + a1);
    float a2a3 = (a2 + a3);
    float b1b2 = (b1 + b2); 

    int bufferSizePerThread = (Width > Height ? Width : Height) * Channels;
    unsigned char * bufferPerLine = (unsigned char*)malloc(bufferSizePerThread);
    unsigned char * tempData = (unsigned char*)malloc(Height * Stride);
    if (bufferPerLine == NULL || tempData == NULL)
    {
        if (tempData)
        {
            free(tempData);
        }
        if (bufferPerLine)
        {
            free(bufferPerLine);
        }
        return;
    }
    for (int y = 0; y < Height; ++y) {
        unsigned char * lpRowInitial = input + Stride * y;
        unsigned char * lpColInitial = tempData + y * Channels;
        gaussianHorizontal(bufferPerLine, lpRowInitial, lpColInitial, Width, Height, Channels, Width, a0a1, a2a3, b1b2, cprev, cnext);
    }
    int HeightStep = Height*Channels;
    for (int x = 0; x < Width; ++x) {
        unsigned char * lpColInitial = output + x*Channels;
        unsigned char * lpRowInitial = tempData + HeightStep * x;
        gaussianVertical(bufferPerLine, lpRowInitial, lpColInitial, Height, Width, Channels, a0a1, a2a3, b1b2, cprev, cnext);
    }

    free(bufferPerLine);
    free(tempData);
}

题外话:

调用方法:

相关文章