图像处理中随心所欲核卷积(matlab中conv2函数)的长足实现。SSE图像算法优化系列十一:使用FFT变换实现图像卷积。

   
 卷积其实是图像处理着最为核心的操作,我们普遍的一部分算法比如:均值模糊、高斯模糊、锐化、Sobel、拉普拉斯、prewitt边缘检测等等一些及世界有关的算法,都可以透过卷积算法实现。只不过由这些算法的卷积矩阵的特殊性,一般不见面直接实现其,而是通过有些优化的招数为计算量变多少。但是来头情况下卷积矩阵的元素值无大规律或发特殊要求,无法通过正规手段优化,这个时节只得通过原始的措施实现。因此,如何快速的兑现图像的随机卷积矩阵操作为时有发生必不可少举行适度的研讨。

     
本文重点要害不在于FFT的SSE优化,而介于使FFT实现高效卷积的连带技能与进程。

     
目前,通过友好共享或自己寻找找到的如出一辙切开关于自由核算法优化的章有: Reshuffling:
A Fast Algorithm for Filtering with Arbitrary
Kernels,改文章称能提高原始程序速度之40%左右,但是旧的程序是怎样勾勒的啊尚未掌握。

      关于FFT变换,有过多参阅的代码,特别是对长度为2之整数次于幂的序列,实现起来也是非常简单的,而对非2次幂的行,就有些有硌麻烦了,matlab中凡是得兑现自由长度FFT的,FFTW也是好的,而Opencv则发选择性的贯彻了少数长度序列的更换,查看Opencv的代码,可以窥见那个只有对凡4底整数不良幂的数有行使了SSE优化,比如4、16、64、256、1024如此的班部分,因此基4的FFT是极度抢的,而余下的片则依然是平常的C语言实现,我这里根本是拿Opencv的代码抠出来了。

     
在matlab中生几个函数都跟图像卷积有关,比如imfilter就可兑现卷积,或者
conv2也实践,他们的速还是一定快之,比如3000*3000的灰度图,卷积矩阵大小为15*15,在I5的CPU上运行时刻使170ms左右,相当之给力。

     Opencv关于FFT实现的代码在Opencv
3.0\opencv\sources\modules\core\src\dxt.cpp中,代码写的专门复杂,扣取工作吧举行的一对一困难,其基4的SSE优化核心代码如下所示:

     
在Celery的博客中,也涉及了他的优化后底conv2和matlab相当甚至快为matlab,详见http://blog.csdn.net/celerychen2009/article/details/38852105。

// optimized radix-4 transform
template<> struct DFT_VecR4<float>
{
    int operator()(Complex<float>* dst, int N, int n0, int& _dw0, const Complex<float>* wave) const
    {
        int n = 1, i, j, nx, dw, dw0 = _dw0;
        __m128 z = _mm_setzero_ps(), x02=z, x13=z, w01=z, w23=z, y01, y23, t0, t1;
        Cv32suf t; t.i = 0x80000000;
        __m128 neg0_mask = _mm_load_ss(&t.f);
        __m128 neg3_mask = _mm_shuffle_ps(neg0_mask, neg0_mask, _MM_SHUFFLE(0,1,2,3));

        for( ; n*4 <= N; )
        {
            nx = n;
            n *= 4;
            dw0 /= 4;

            for( i = 0; i < n0; i += n )
            {
                Complexf *v0, *v1;

                v0 = dst + i;
                v1 = v0 + nx*2;

                x02 = _mm_loadl_pi(x02, (const __m64*)&v0[0]);
                x13 = _mm_loadl_pi(x13, (const __m64*)&v0[nx]);
                x02 = _mm_loadh_pi(x02, (const __m64*)&v1[0]);
                x13 = _mm_loadh_pi(x13, (const __m64*)&v1[nx]);

                y01 = _mm_add_ps(x02, x13);
                y23 = _mm_sub_ps(x02, x13);
                t1 = _mm_xor_ps(_mm_shuffle_ps(y01, y23, _MM_SHUFFLE(2,3,3,2)), neg3_mask);
                t0 = _mm_movelh_ps(y01, y23);
                y01 = _mm_add_ps(t0, t1);
                y23 = _mm_sub_ps(t0, t1);

                _mm_storel_pi((__m64*)&v0[0], y01);
                _mm_storeh_pi((__m64*)&v0[nx], y01);
                _mm_storel_pi((__m64*)&v1[0], y23);
                _mm_storeh_pi((__m64*)&v1[nx], y23);

                for( j = 1, dw = dw0; j < nx; j++, dw += dw0 )
                {
                    v0 = dst + i + j;
                    v1 = v0 + nx*2;

                    x13 = _mm_loadl_pi(x13, (const __m64*)&v0[nx]);
                    w23 = _mm_loadl_pi(w23, (const __m64*)&wave[dw*2]);
                    x13 = _mm_loadh_pi(x13, (const __m64*)&v1[nx]); // x1, x3 = r1 i1 r3 i3
                    w23 = _mm_loadh_pi(w23, (const __m64*)&wave[dw*3]); // w2, w3 = wr2 wi2 wr3 wi3

                    t0 = _mm_mul_ps(_mm_moveldup_ps(x13), w23);
                    t1 = _mm_mul_ps(_mm_movehdup_ps(x13), _mm_shuffle_ps(w23, w23, _MM_SHUFFLE(2,3,0,1)));
                    x13 = _mm_addsub_ps(t0, t1);
                    // re(x1*w2), im(x1*w2), re(x3*w3), im(x3*w3)
                    x02 = _mm_loadl_pi(x02, (const __m64*)&v1[0]); // x2 = r2 i2
                    w01 = _mm_loadl_pi(w01, (const __m64*)&wave[dw]); // w1 = wr1 wi1
                    x02 = _mm_shuffle_ps(x02, x02, _MM_SHUFFLE(0,0,1,1));
                    w01 = _mm_shuffle_ps(w01, w01, _MM_SHUFFLE(1,0,0,1));
                    x02 = _mm_mul_ps(x02, w01);
                    x02 = _mm_addsub_ps(x02, _mm_movelh_ps(x02, x02));
                    // re(x0) im(x0) re(x2*w1), im(x2*w1)
                    x02 = _mm_loadl_pi(x02, (const __m64*)&v0[0]);

                    y01 = _mm_add_ps(x02, x13);
                    y23 = _mm_sub_ps(x02, x13);
                    t1 = _mm_xor_ps(_mm_shuffle_ps(y01, y23, _MM_SHUFFLE(2,3,3,2)), neg3_mask);
                    t0 = _mm_movelh_ps(y01, y23);
                    y01 = _mm_add_ps(t0, t1);
                    y23 = _mm_sub_ps(t0, t1);

                    _mm_storel_pi((__m64*)&v0[0], y01);
                    _mm_storeh_pi((__m64*)&v0[nx], y01);
                    _mm_storel_pi((__m64*)&v1[0], y23);
                    _mm_storeh_pi((__m64*)&v1[nx], y23);
                }
            }
        }

        _dw0 = dw0;
        return n;
    }
};

     
由于matlab的代码中行使及了IPL库进行加速,目前本身写的Conv2函数还无法成功同其一定,对于另外审核速度大约为matlab的一半。

  其实为无复杂,总看里面的_mm_loadh_pi用之不可开交好,这也沾光于该数据源采用了Complex格式,这样数据的实部和虚部在内存中是连的,用SSE加载数据是啊即特别有利了。

     
简单的记录下我在举行卷积过程被因故到之优化吧。

  上面才是有代码,为了配合该过程,还有不少之初始化,比如数据的shuffle等。详见Opencv的文本。

     
原始的卷积的实现用四重循环,简单的表述如下:

  对于2维的FFT变换,我从未去看CV的代码,而是一直先每行进行同样维的FFT1D,然后针对结果在进展排方向的FFT1D,由于FFT1D算法需要处理的排必须是连续的内存,因此,需要针对中等的结果进行转置,处理终结晚当转置回来,为了节省时间,这个转置也相应用SSE优化。

for (Y = 0; Y < Height; Y++)
{
    for (X = 0; X < Width; X++)
    {
        Index = .....;
        Sum = 0;
        for (XX = 0; XX < ConvW; XX++)
        {
            for (YY = 0; YY < ConvH; YY++)
            {
                Index1 = ..... ;
                Index2 = ..... ;
                Sum += Conv[Index1] * Pixel[Index2];
            }
        }
        Dest[Index] = Sum / Weight;
    }
}

  当2维的涨幅与惊人一样时,这个转置是休欲分配另外一卖附加的内存的,这个被In-Place转置,另外一个最主要亮点就是是FFT1D算法为支持In-Place操作。

  当卷积矩阵较充分时,计算量将会见生要命,而且由于程序中的内存访问很频繁,cache
miss现象比较严重,因此效率远低下。

  由于是指向Complex数据进行转置,我们可以转换个角度考虑,因为一个Complex正好跟double数据占用同样的内存,这样直接运用以及double相关的SSE指令就可方便之兑现Complex相关数据的转置。比如当宽度和冲天一致不时,刻直接以下述方式。

   
 我之优化措施主要概括以下几独面: 

inline void Inplace_TransposeSSE2X2D(double *Src, double *Dest, int Length)
{
    __m128d S0 = _mm_loadu_pd(Src);                        
    __m128d S1 = _mm_loadu_pd((Src + Length));            
    __m128d D0 = _mm_loadu_pd(Dest);                
    __m128d D1 = _mm_loadu_pd((Dest + Length));            
    _mm_storeu_pd(Dest, _mm_unpacklo_pd(S0, S1));
    _mm_storeu_pd(Dest + Length, _mm_unpackhi_pd(S0, S1));
    _mm_storeu_pd(Src, _mm_unpacklo_pd(D0, D1));
    _mm_storeu_pd(Src + Length, _mm_unpackhi_pd(D0, D1));
}

   
 一:使用SSE进行乘法计算,由于SSE可以一次性进行4独单精度浮点数的精打细算,因此可发鲜明的速度提升。

  参考系统的Intrinsic的_MM_TRANSPOSE4_PS函数,可以窥见double类型的转置方式不极端相同,可以直接采用有关的unpack函数,而非应用shuffle,当然使用shuffle也是未曾问题之。

   
 二:通过当的处理方式,对每个取样点周边的卷积矩阵内的素进行汇总,使得各个动一个如素点不见面待打内存中展开大气底寻找工作。

  那么最终我实现之FFT2D函数如下所示:

     具体来说实现过程如下:

int FFT2D(Complex *Input, Complex *Output, int Width, int Height, bool Invert, int StartZeroLines = 0, int EndZeroLines = 0)
{
    if ((Input == NULL) || (Output == NULL))                                return IM_STATUS_NULLREFRENCE;
    if ((IsPowerOfTwo(Width) == false) || (IsPowerOfTwo(Height) == false))    return IM_STATUS_INVALIDPARAMETER;

    if (Width == Height)
    {
        void *Buffer = FFT_Init(Width);
        if (Buffer == NULL)    return IM_STATUS_OUTOFMEMORY;

        for (int Y = StartZeroLines; Y < Height - EndZeroLines; Y++)
            FFT1D(Input + Y * Width, Output + Y * Width, Buffer, Width, Invert);            //    先水平方向变换
        InPlace_TransposeD((double *)Output, Width);                                        //    再位转置

        for (int Y = 0; Y < Height; Y++)
            FFT1D(Output + Y * Width, Output + Y * Width, Buffer, Width, Invert);            //    再垂直方向变换
        InPlace_TransposeD((double *)Output, Width);
        FFT_Free(Buffer);
    }
    else
    {
        void *Buffer = FFT_Init(Width);
        if (Buffer == NULL)    return IM_STATUS_OUTOFMEMORY;

        for (int Y = StartZeroLines; Y < Height - EndZeroLines; Y++)
            FFT1D(Input + Y * Width, Output + Y * Width, Buffer, Width, Invert);
        Complex *Temp = (Complex *)malloc(Width * Height * sizeof(Complex));
        if (Temp == NULL)
        {
            FFT_Free(Buffer);
            return IM_STATUS_OUTOFMEMORY;
        }
        TransposeD((double *)Output, (double *)Temp, Width, Height);
        FFT_Free(Buffer);

        Buffer = FFT_Init(Height);                                                            //    必须要二次初始化
        if (Buffer == NULL)    return IM_STATUS_OUTOFMEMORY;

        for (int Y = 0; Y < Width; Y++)
            FFT1D(Temp + Y * Height, Temp + Y * Height, Buffer, Height, Invert);
        TransposeD((double *)Temp, (double *)Output, Height, Width);
        FFT_Free(Buffer);
        free(Temp);
    }
    return IM_STATUS_OK;
}

         
 1、为了用SSE的优势,首先以卷积矩阵进行调整,调整卷积矩阵一行的元素个数,使其为免小于原始值的4的平头加倍,并且于新的卷积矩阵的内存布局符合SSE相关函数的16字节本着同的要求。

  注意在函数的终极自己多了StartZeroLines
和EndZeroLines
两个参数,他们要是在展开实施方向的FFT1D是忽视最前的StartZeroLines和最后给之EndZeroLines
个全0的推行,因为全0的变换后面的结果还是0,没有必要展开测算,减少计算时间。在opencv中为是产生个nonzero_rows的,但是他单纯是指向前的行,而其实被,比如很多情况是,先使壮大数据,然后拿中数据放置在左上角,这样就出右下比赛有非零元素,所以opencv这样做就是无法从及加速作用了。

  实现代码如下:

  对扣取的代码进行了实在的测试,1024*1024之数,进行100次于正反变换,耗时3000ms,使用matlab进行相同的操作,耗时大体5500ms,并且观察任务管理器,在4核对PC上CPU使用率100%,说明外中采用了大半线程,不过起某些纵是matlab使用的是double类型的数量。听说matlab最新版本使用的就算是FFTW库,不过不管怎样,这个速度还是得承受以及一对一快之。

float *Conv16 = (float *)_mm_malloc(PadConvLine * ConvH * sizeof(float), 16);                        //    保存16字节对齐的卷积矩阵,以方便使用SSE

for(Y = 0; Y < ConvH; Y++) 
{
    memcpy (Conv16 + Y * PadConvLine, Conv->Data.F + Y * ConvW , ConvW * sizeof(float));            //    复制卷积矩阵的数据
    memset(Conv16 + Y * PadConvLine + ConvW, 0, (PadConvLine - ConvW) * sizeof(float));                //    把冗余部分的卷积数据设置为0
}

  下面我们要出口下基于FFT的图像卷积的贯彻,理论及使图像a大小也N
* M,卷积核b大小为 X * Y,则卷积实现的历程如下:

    其中PadConvLine = Pad4(ConvW)
以及Pad4的原型为: #define Pad4(bits) (((bits) + 3) / 4 * 4);

  首先扩展数据,扩展后底大小也 (N + X –
1) * (M + Y –
1),将卷积核数据放置到扩大后底数码的左上角,其他因素填充0,得到bb,
对bb进行FFT2D正朝着易得到B,然后呢拿图像数据放置到图像的左上角,其他因素填充为0,得到aa,对aa也拓展FFT2D正为易得到B,接着对A和B进行接触就得到C,最后对C进行逆向的FFT变换得到D,最后取D的中档部分中数据就是是卷积的结果。

   
注意_mm_malloc函数分配的内存中之价值是随机值,对于扩大的有的一定要填写充0,否则就见面损坏卷积的结果。

  举个例子,假要图像数据也:

   
那么只要我们也还要获了需要被卷积的组成部分数据的言语(卷积核肯定与卷积矩阵一样大小,且为理应是16字节针对同之),可以用如下的SSE的代码进行乘法计算:

       大红鹰葡京会娱乐 1

float MultiplySSE(float *Kernel, float *Conv, int Length)
{
    int Block;  
    const float *Data;                        // 将SSE变量上的多个数值合并时所用指针.
    float Sum = 0;
    if (Length > 16)                        //    可以进行四次SSE计算,测试表明,这个还是快些的    
    {
        const int BlockWidth = 4 * 4;        // 块宽. SSE寄存器能一次处理4个float,然后循环展开4次.
        Block = Length / BlockWidth;        // 块数.    
        float *KernelP = Kernel, *ConvP = Conv;                // SSE批量处理时所用的指针.

        __m128 Sum0 = _mm_setzero_ps();         // 求和变量。SSE赋初值0
        __m128 Sum1 = _mm_setzero_ps();
        __m128 Sum2 = _mm_setzero_ps();
        __m128 Sum3 = _mm_setzero_ps();

        for(int I = 0; I < Block; I++)
        {
            Sum0 = _mm_add_ps(Sum0, _mm_mul_ps(_mm_load_ps(KernelP), _mm_load_ps(ConvP)));                    // SSE单精浮点紧缩加法
            Sum1 = _mm_add_ps(Sum1, _mm_mul_ps(_mm_load_ps(KernelP + 4), _mm_load_ps(ConvP + 4)));
            Sum2 = _mm_add_ps(Sum2, _mm_mul_ps(_mm_load_ps(KernelP + 8), _mm_load_ps(ConvP + 8)));
            Sum3 = _mm_add_ps(Sum3, _mm_mul_ps(_mm_load_ps(KernelP + 12), _mm_load_ps(ConvP + 12)));
            KernelP += BlockWidth;
            ConvP += BlockWidth;
        }

        Sum0 = _mm_add_ps(Sum0, Sum1);    // 两两合并(0~1).
        Sum2 = _mm_add_ps(Sum2, Sum3);    // 两两合并(2~3).
        Sum0 = _mm_add_ps(Sum0, Sum2);    // 两两合并(0~2).

        Data = (const float *)&Sum0;
        Sum = Data[0] + Data[1] + Data[2] + Data[3];

        Length = Length - Block * BlockWidth;            // 剩余数量.
    }
    if (Length != 0)
    {
        const int BlockWidth = 4;                        //    程序已经保证了数量必然是4的倍数 
        Block = Length / BlockWidth;        
        float *KernelP = Kernel, *ConvP = Conv;                
        __m128 Sum0 = _mm_setzero_ps();        

        for(int I = 0; I < Block; I++)
        {
            Sum0 = _mm_add_ps(Sum0, _mm_mul_ps(_mm_load_ps(KernelP), _mm_load_ps(ConvP)));        
            KernelP += BlockWidth;
            ConvP += BlockWidth;
        }

        Data = (const float *)&Sum0;
        Sum += Data[0] + Data[1] + Data[2] + Data[3];
    }
    return Sum;
}

  卷积核为:

  当卷积矩阵(扩充后)的因素数量超过16时,我们下了4行程并行的SSE乘法实现,我当I3的CPU上测试时,2路SSE和4路SSE已经没有什么异常的别了,而于I5的CPU上则4路还是有比较明显的增进,因此采取4路SSE同时运行。当然1路SSE肯定还是比2路慢。另外,如果元素的多少有限16还是过16而无可知叫16规整除,那么余下的有些由于以前的扩展,剩余元素数量也势必是4底倍数,因此好为此单路的SSE实现。
这吗是编码上的技艺。

         大红鹰葡京会娱乐 2

   
 2、前面提到了用给卷积的有些数据,这一部分哪高效的取得呢。观察最原始之4重循环,其里面的2重就算为获得需要吃卷积的片,但是这里实在产生无数题目。第一:由于卷积取样时得有一些取样点的坐标在本来图像的有效限制外,因此须开展判定,耗时。第二:同样为采取SSE,也要把取样的多少在和壮大的卷积矩阵一样大小的内存中。这里我事先贴出自己的代码在拓展分解具体的兑现:

  扩展后底图像数据吧:

IS_RET __stdcall Conv2(TImage *Src, TMatrix *Conv, TImage *Dest, EdgeMode Edge)
{
    if (Src == NULL || Dest == NULL || Conv == NULL) return IS_RET_ERR_PARA;
    if (Src->Width != Dest->Width || Src->Height != Dest->Height || Src->BitCount != Dest->BitCount || Src->Stride != Dest->Stride) return IS_RET_ERR_PARA;
    if (Src->Scan0 == NULL || Dest->Scan0 == NULL || Conv->Data.F == NULL) return IS_RET_ERR_MEM;
    if (Conv->Width < 1 || Conv->Height < 1) return IS_RET_ERR_PARA;

    int Width = Src->Width, Height = Src->Height, Stride = Src->Stride;
    int ConvW = Conv->Width, ConvH = Conv->Height;
    unsigned char *PtSrc = Src->Scan0, *PtDest = Dest->Scan0;


    if (Src->BitCount == 24)
    {

    }
    else
    {
        int Left = ConvW / 2, Top = ConvH / 2, Right = ConvW - Left - 1, Bottom = ConvH - Top - 1, ExpHeight = Height + ConvH - 1;        //    注意核中心那个元素不用扩展,比如核的宽度为3,则只要左右各扩展一个像素就可以了
        int PadConvLine = Pad4(ConvW), Length = PadConvLine * ConvH;
        int X, Y, IndexD, IndexE, IndexK, ExpStride;
        float *CurKer, Inv, Sum = 0;
        unsigned char *PtExp, *PtDest;

        TImage *Expand;
        IS_RET Ret = GetPadImage(Src, &Expand, Left, Right, Top, Bottom, Edge);                                //    得到扩展后的数据,可以提速和方便编程,但是多占用一份内存
        if (Ret != IS_RET_OK) return Ret;

        PtExp = Expand->Scan0; PtDest = Dest->Scan0; ExpStride = Expand->Stride;

        for (X = 0; X < ConvH * ConvW; X ++) Sum += Conv->Data.F[X];
        Inv = (Sum == 0 ? 1: 1 / Sum);                                                                        //    如果卷积举证的和为0,则设置为1

        float *Conv16 = (float *)_mm_malloc(PadConvLine * ConvH * sizeof(float), 16);                        //    保存16字节对齐的卷积矩阵,以方便使用SSE
        float *Kernel = (float *)_mm_malloc(PadConvLine * ExpHeight * sizeof(float), 16);                    //    保存16字节对齐的卷积核矩阵,以方便使用SSE

        for(Y = 0; Y < ConvH; Y++) 
        {
            memcpy (Conv16 + Y * PadConvLine, Conv->Data.F + Y * ConvW , ConvW * sizeof(float));            //    复制卷积矩阵的数据
            memset(Conv16 + Y * PadConvLine + ConvW, 0, (PadConvLine - ConvW) * sizeof(float));                //    把冗余部分的卷积数据设置为0
        }

        for (Y = 0; Y < ExpHeight; Y++)
        {
            IndexE = Y * ExpStride;
            CurKer = Kernel + Y * PadConvLine;                        //    计算第一列所有像素将要取样的卷积核数据
            for (X = 0; X < ConvW; X++)
            {
                CurKer[X] = PtExp[IndexE++];
            }
        }
        for (X = 0 ; X < Width ; X ++)
        {
            if (X != 0)                                                //    如果不是第一列,需要更新卷积核的数据
            {
                memcpy(Kernel, Kernel + 1, (PadConvLine * ExpHeight - 1) * sizeof(float));    //    往前移动一个数据
                IndexK = ConvW - 1 ;
                IndexE = IndexK + X;
                for (Y = 0; Y < ExpHeight; Y++)
                {
                    Kernel[IndexK] = PtExp[IndexE];        //    只要刷新下一个元素
                    IndexK += PadConvLine;
                    IndexE += ExpStride;
                }
            }

            CurKer = Kernel;    IndexD = X;
            for (Y = 0; Y < Height; Y ++)                            //    沿列的方向进行更新
            {
                PtDest[IndexD] = Clamp((int)( MultiplySSE(Conv16, CurKer, Length) * Inv + 0.5));        //    直接把函数放在这里也没有啥提速的,注意改函数不会被内联的
                CurKer += PadConvLine;
                IndexD += Stride;
            }
        }
        _mm_free(Conv16);
        _mm_free(Kernel);
        FreeImage(Expand);
        return IS_RET_OK;
    }
}

        大红鹰葡京会娱乐 3

   
 对于第一单问题,解决之艺术要命简答,即用空间更换时间,新建一适合(Width +
ConvW – 1, Height + ConvH
-1)大小的图像,然后四周的ConvW及ConvH的例如素用边缘之价或者边缘镜像的值填充,正中间的尽管就此原来的图复制过来,这样操作后开展抽样时不再原图取样,而于马上福扩展的希冀被抽样,就避免了坐标判断等if语句的跳转耗时了,上GetPadImage即实现了改观功能。

  扩展后的卷积数据也:

     第二单问题虽然需出自然之兑现技能,我们分配一片PadConvLine
* (Height + ConvH – 1)
大小的内存,然后计算原图第一排像素串联起来的用卷积的片的数码,这同局部代码如上述44-52实施所展示。有矣这么的多寡,如果要计算第一排列的卷积结果,则生简单了,每超过了相同排则拿给卷积的数起点增加PadConvLine个因素,在调用上述MultiplySSE函数获得卷积结果。接着则计算第二列像从的卷积值,此时用总体创新就等同排列像素串联起来的急需给卷积的数据,更新也酷粗略,就是把本来的数量完全往左移动一个像素,这个好为此memcpy快速实现,然后在填充充入新入的老元素,就ok了,接着就是再次调用MultiplySSE函数,如此重复下去。

        大红鹰葡京会娱乐 4

   
 经过编码测试,对于3000*3000之灰度图,15*15底审批在I5的CPU上之测试平均结果为360ms,比matlab的徐了大体上。

  进行上述操作:D =
ifft2(fft2(aa).*fft2(bb)),得到:

   
 最后验明正身某些,很多丁还说用FFT可以快速的兑现卷积,并且是O(1)的,我较同意后半句,但是前半句是纯属的发出题目之,至少在甄别小于50*50常,FFT实现之卷积不见面比直接实现块。要理解FFT的计算量其实是非常非常的。

        大红鹰葡京会娱乐 5

   
 大红鹰葡京会娱乐 6

  中间部分的数码就是是卷积的实用数据。

****************************笔者:
laviewpbt   时间: 2014.11.27    联系QQ:  33184777
转载请保留本行信息**********************

  仔细观察,可以发现这样的卷积在边缘处是产生问题,他管过边缘有的多寡总体用0代替了,因此这种卷积还是无全面的,一栽解决办法就是第一将原图沿着边缘扩展,扩展的点子得以是双重或镜像,扩展的高低当然就是和卷积核的大小有关了,一般遵循卷积核中心对称分布,这时,边缘扩展后底图像大小也  (N

 

  • X – 1) * (M + Y – 1)。

  也就是说合理之过程是进行简单次扩大,在进展FFT变换前的尾声数额维度为 (N

  • X – 1 + X – 1) * (M + Y – 1 + Y – 1),在拓展逆变换得到的结果遭到从第(X
  • 1, Y – 1)坐标处取有效值。

  我们在精心审视下上述过程,第一、内存占用方面,需要2独 (N

  • X – 1 + X – 1) * (M + Y – 1 + Y – 1) *
    sizeof(Complex)大小的内存,这老惊人之,特备是对此移动装备,第二,由于我们目前之FFT都是不得不处理2的平头不成幂的排,如果因为全图也目标,则有或会见多很怪的失效计算,比如N
  • X – 1 + X – 1
    万一当1025,则需要调动有关2048,尺寸约好,这种无效计算的可能性越来越怪,而当时为我们面前希望使2之整数不好幂的FFT最抢之初衷就是矛盾了。

  一栽缓解办法就是是分块计算,比如我们将图像分成很多独满足条件 (NN+
X – 1 + X – 1)  = 256 和  (MM + Y – 1 + Y – 1) = 256底块,其中NN *
MM就是图像分块大之轻重缓急,这样对卷积核那一些的FFT就是同等差256*256之第二维卷积,并且是公共的,这较原来的算计同一片大死之(N

  • X – 1 + X – 1) * (M + Y – 1 + Y –
    1)要节约成千上万时,同时,只占了颇粗之同块内存。当卷积核的分寸不超越50时常,每次有效之盘算的块NN
    * MM相对于完整的2D
    FFT计算的话占比要相当高之。这样可使得的滑坡1025尺码直接成为了2048这么的FFT计算。

  另外注意一点,FFT卷积是虚部和实部的图是同的,也就是说我们得以又展开个别个未思关元素的测算,比如对于32号图像,可以把一个块的Blue分量填充到实部,把Green分量填充到虚部,这样一次性就完了了2个重的乘除。而对于灰度图,则足以同时计算两只片,即首先独片的灰度值填充到实部,第二个片的灰度填充到虚部。这样一次性就到位了2只片的算计。

  还有,在每个片的底层,有X – 1 + X
-1只实施是只要填为0的,因此这些的执行之一律维FFT变换是无必要的,可以优化掉。

  优化后底局部代码如下:

int IM_FFTConv2(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, float *Kernel, int KerWidth, int KerHeight)                //    阈值
{
    int Channel = Stride / Width;
    if ((Src == NULL) || (Dest == NULL))                        return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))                            return IM_STATUS_INVALIDPARAMETER;
    if ((KerWidth >50) || (KerHeight > 50))                        return IM_STATUS_INVALIDPARAMETER;            // 卷积核越大,每次的有效计算量就越小了
    if ((Channel != 1) && (Channel != 3) && (Channel != 4))        return IM_STATUS_INVALIDPARAMETER;

    int Status = IM_STATUS_OK;

    const int TileWidth = 256, TileHeight = 256;

    int HalfW = KerWidth / 2, HalfH = KerHeight / 2;
    int ValidW = TileWidth + 1 - KerWidth - (KerWidth - 1);                                //    第一个加1是因为卷积扩展的矩阵大小为N+M-1,第二个-1是因为为了取得有效数据,还要对边缘进行扩展,扩展的大小为KerHeight - 1,比如KerHeight=5,则每边需要扩展2个像素,一共扩展4个像素
    int ValidH = TileHeight + 1 - KerHeight - (KerHeight - 1);                            //    默认的卷积的效果再卷积周边是用0来填充的,如果分块处理时,这明显是不能满足要求的,会带来明显的块于块之间的分界线

    int TileAmountX = (Width / ValidW) + (Width % ValidW ? 1 : 0);                        //    图像需要分成的块数
    int TileAmountY = (Height / ValidH) + (Height % ValidH ? 1 : 0);

    Complex *Conv = (Complex *)malloc(TileWidth * TileHeight * sizeof(Complex));        //    需要将卷积核扩展到TileWidth和TileHeight大小                    
    Complex *Tile = (Complex *)malloc(TileWidth * TileHeight * sizeof(Complex) * 3);    //    每个小块对应的数据,当为24位模式时需要3份内存,灰度只需一份

    int *RowOffset = (int *)malloc((TileAmountX * ValidW + KerWidth) * sizeof(int));    //    每个小块取样时的坐标偏移,这样在中间的块也可以取到周边合理的只,在边缘处则位镜像值
    int *ColOffset = (int *)malloc((TileAmountY * ValidH + KerHeight) * sizeof(int));

    if ((Conv == NULL) || (Tile == NULL) || (RowOffset == NULL) || (ColOffset == NULL))
    {
        Status = IM_STATUS_OUTOFMEMORY;
        goto FreeMemory;
    }
    Status = IM_GetOffsetPos(RowOffset, Width, HalfW, TileAmountX * ValidW + (KerWidth - HalfW) - Width, IM_EDGE_MIRROR);            //    左右对称
    if (Status != IM_STATUS_OK) goto FreeMemory;
    Status = IM_GetOffsetPos(ColOffset, Height, HalfH, TileAmountY * ValidH + (KerHeight - HalfH) - Height, IM_EDGE_MIRROR);
    if (Status != IM_STATUS_OK) goto FreeMemory;

    memset(Conv, 0, TileWidth * TileHeight * sizeof(Complex));                //    卷积核的其他元素都为0,这里先整体赋值为0
    for (int Y = 0; Y < KerHeight; Y++)
    {
        int Index = Y * KerWidth;
        for (int X = 0; X < KerWidth; X++)
        {
            Conv[Y * TileWidth + X].Real = Kernel[Index + X];                //    卷积核需要放置在左上角
        }
    }

    Status = FFT2D(Conv, Conv, TileWidth, TileHeight, false, 0, TileHeight - KerHeight);                //    对卷积核进行FFT变换,注意行方向上下部都为0,这样可以节省部分计算时间    
    if (Status != IM_STATUS_OK)    goto FreeMemory;

    if (Channel == 1)                                                                                    //    单通道时可以一次性处理2个块
    {

    }
    else if (Channel == 3)                                                                // 3通道时可以利用两个块的信息分别填充到BG  RB  GR 序列里,这样就更快了
    {

    }
    else if (Channel == 4)
    {
        Complex *TileBG = Tile, *TileRA = Tile + TileWidth * TileHeight;
        for (int TileY = 0; TileY < TileAmountY; TileY++)
        {
            for (int TileX = 0; TileX < TileAmountX; TileX++)
            {
                IM_Rectangle SrcR, ValidR;
                IM_SetRect(&SrcR, TileX * ValidW, TileY * ValidH, TileX * ValidW + ValidW, TileY * ValidH + ValidH);
                IM_SetRect(&ValidR, 0, 0, Width, Height);
                IM_IntersectRect(&ValidR, ValidR, SrcR);
                for (int Y = ValidR.Top; Y < ValidR.Bottom + KerHeight - 1; Y++)
                {
                    byte *LinePS = Src + ColOffset[Y] * Stride;
                    int Index = (Y - ValidR.Top) * TileWidth;
                    for (int X = ValidR.Left; X < ValidR.Right + KerWidth - 1; X++)
                    {
                        byte *Sample = LinePS + RowOffset[X] * 4;
                        TileBG[Index].Real = Sample[0];
                        TileBG[Index].Imag = Sample[1];
                        TileRA[Index].Real = Sample[2];
                        TileRA[Index].Imag = Sample[3];
                        Index++;
                    }
                }
                Status = FFT2D(TileBG, TileBG, TileWidth, TileHeight, false, 0, TileHeight - (ValidR.Bottom + KerHeight - 1 - ValidR.Top));
                if (Status != IM_STATUS_OK)    goto FreeMemory;

                Status = FFT2D(TileRA, TileRA, TileWidth, TileHeight, false, 0, TileHeight - (ValidR.Bottom + KerHeight - 1 - ValidR.Top));
                if (Status != IM_STATUS_OK)    goto FreeMemory;

                for (int Y = 0; Y < TileWidth * TileHeight; Y++)
                {
                    float Temp = TileBG[Y].Real;
                    TileBG[Y].Real = TileBG[Y].Real * Conv[Y].Real - TileBG[Y].Imag * Conv[Y].Imag;        
                    TileBG[Y].Imag = Temp * Conv[Y].Imag + TileBG[Y].Imag * Conv[Y].Real;
                    Temp = TileRA[Y].Real;
                    TileRA[Y].Real = TileRA[Y].Real * Conv[Y].Real - TileRA[Y].Imag * Conv[Y].Imag;
                    TileRA[Y].Imag = Temp * Conv[Y].Imag + TileRA[Y].Imag * Conv[Y].Real;
                }
                FFT2D(TileBG, TileBG, TileWidth, TileHeight, true);
                if (Status != IM_STATUS_OK)    goto FreeMemory;

                FFT2D(TileRA, TileRA, TileWidth, TileHeight, true);
                if (Status != IM_STATUS_OK)    goto FreeMemory;

                for (int Y = ValidR.Top; Y < ValidR.Bottom; Y++)
                {
                    byte *LinePD = Dest + Y * Stride + ValidR.Left * 4;
                    int Index = (Y - ValidR.Top + KerHeight - 1) * TileWidth + KerWidth - 1;
                    for (int X = ValidR.Left; X < ValidR.Right; X++)
                    {
                        LinePD[0] = IM_ClampToByte(TileBG[Index].Real);                        //    取有效的数据
                        LinePD[1] = IM_ClampToByte(TileBG[Index].Imag);
                        LinePD[2] = IM_ClampToByte(TileRA[Index].Real);
                        LinePD[3] = IM_ClampToByte(TileRA[Index].Imag);
                        LinePD += 4;
                        Index++;
                    }
                }
            }
        }
    }
FreeMemory:
    if (RowOffset != NULL) free(RowOffset);
    if (ColOffset != NULL) free(ColOffset);
    if (Conv != NULL) free(Conv);
    if (Tile != NULL) free(Tile);
    return Status;
}

 

  程序里我定位每个片的高低为256*256,这至关重要是考虑256是4底平头次幂,是力所能及统统用SSE优化掉的,同时为占有更不见之内存。

  如果考虑重新好的拍卖,在结尾一排列以及尾声一执不行块,因为中元素的滑坡,可以考虑使用重复小之丘来计算,不过当下会增加程序的繁杂。本例没有设想了。

  速度测试:

  3000*3000的灰度图    卷积核5*5       
 328ms

       
             卷积核15*15        348ms

              
 卷积核25*25       400ms

              
 卷积核49*49       600ms

  以前我勾勒了因SSE直接卷积,同一个机器上呢申报下测试速度:

  3000*3000底灰度图    卷积核5*5       
 264ms

       
             卷积核15*15       550ms

              
 卷积核25*25       1300ms

  所以我以前那么篇文章的有点结论是漏洞百出的。

  测试工程的地点:http://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar

 

 大红鹰葡京会娱乐 7

相关文章

admin

网站地图xml地图