使用FFT变换实现图像卷积

Imageshop    2018-8-31

本文重点主要不在于FFT的SSE优化,而在于使用FFT实现快速卷积的相关技巧和过程。

关于FFT变换,有很多参考的代码,特别是对于长度为2的整数次幂的序列,实现起来也是非常简易的,而对于非2次幂的序列,就稍微有点麻烦了,matlab中是可以实现任意长度FFT的,FFTW也是可以的,而Opencv则有选择性的实现了某些长度序列的变换,查看Opencv的代码,可以发现其只有对是4的整数次幂的数据部分采用了SSE优化,比如4、16、64、256、1024这样的序列部分,因此基4的FFT是最快的,而剩余的部分则依旧是普通的C语言实现,我这里主要是把Opencv的代码抠出来了。

Opencv关于FFT实现的代码在Opencv 3.0\opencv\sources\modules\core\src\dxt.cpp中,代码写的特别复杂,扣取工作也做的相当艰苦,其基4的SSE优化核心代码如下所示:

//optimized radix-4 transformtemplate<>structDFT_VecR4<float>{intoperator()(Complex<float>* dst,intN,intn0,int& _dw0,constComplex<float>* wave)const{intn =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 i3w23 = _mm_loadh_pi(w23, (const__m64*)&wave[dw*3]);//w2, w3 = wr2 wi2 wr3 wi3t0=_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 i2w01 = _mm_loadl_pi(w01, (const__m64*)&wave[dw]);//w1 = wr1 wi1x02 = _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;returnn;
    }
};

其实也不复杂,总觉得里面的_mm_loadh_pi用的很好,这也得益于其数据源采用了Complex格式,这样数据的实部和虚部在内存中是连续的,用SSE加载数据是也就很方便了。

上面只是部分代码,为了配合该过程,还有很多的初始化,比如数据的shuffle等。详见Opencv的文件。

对于2维的FFT变换,我没有去扣CV的代码,而是直接先每行进行一维的FFT1D,然后对结果在进行列方向的FFT1D,由于FFT1D算法需要处理的序列必须是连续的内存,因此,需要对中间的结果进行转置,处理完后在转置回来,为了节省时间,这个转置也应该用SSE优化。

当2维的宽度和高度相同时,这个转置是不需要分配另外一份额外的内存的,这个叫In-Place转置,另外一个重要优点就是FFT1D算法也支持In-Place操作。

由于是对Complex数据进行转置,我们可以换个角度考虑,因为一个Complex正好和double数据占用同样的内存,这样直接利用和double相关的SSE指令就可以方便的实现Complex相关数据的转置。比如当宽度和高度一样时,刻直接使用下述方式。

inlinevoidInplace_TransposeSSE2X2D(double*Src,double*Dest,intLength)
{
    __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));
}

参考系统的Intrinsic的_MM_TRANSPOSE4_PS函数,可以发现double类型的转置方式不太一样,可以直接使用有关的unpack函数,而不使用shuffle,当然使用shuffle也是没有问题的。

那么最后我实现的FFT2D函数如下所示:

intFFT2D(Complex *Input, Complex *Output,intWidth,intHeight,boolInvert,intStartZeroLines =0,intEndZeroLines =0)
{if((Input == NULL) || (Output == NULL))returnIM_STATUS_NULLREFRENCE;if((IsPowerOfTwo(Width) ==false) || (IsPowerOfTwo(Height) ==false))returnIM_STATUS_INVALIDPARAMETER;if(Width ==Height)
    {void*Buffer =FFT_Init(Width);if(Buffer == NULL)returnIM_STATUS_OUTOFMEMORY;for(intY = StartZeroLines; Y < Height - EndZeroLines; Y++)
            FFT1D(Input+ Y * Width, Output + Y * Width, Buffer, Width, Invert);//先水平方向变换InPlace_TransposeD((double*)Output, Width);//再位转置for(intY =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)returnIM_STATUS_OUTOFMEMORY;for(intY = 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);returnIM_STATUS_OUTOFMEMORY;
        }
        TransposeD((double*)Output, (double*)Temp, Width, Height);
        FFT_Free(Buffer);

        Buffer= FFT_Init(Height);//必须要二次初始化if(Buffer == NULL)returnIM_STATUS_OUTOFMEMORY;for(intY =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);
    }returnIM_STATUS_OK;
}

注意在函数的最后我增加了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库,不过无论如何,这个速度还是可以接受和相当快的。

下面我们重点谈下基于FFT的图像卷积的实现,理论上如果图像a大小为N * M,卷积核b大小为 X * Y,则卷积实现的过程如下:

首先扩展数据,扩展后的大小为 (N + X - 1) * (M + Y - 1),将卷积核数据放置到扩展后的数据的左上角,其他元素填充0,得到bb, 对bb进行FFT2D正向变换得到B,然后也将图像数据放置到图像的左上角,其他元素填充为0,得到aa,对aa也进行FFT2D正向变换得到B,接着对A和B进行点乘得到C,最后对C进行逆向的FFT变换得到D,最后取D的中间部分有效数据就是卷积的结果。

举个例子,假设图像数据为:

卷积核为:

扩展后的图像数据为:

扩展后的卷积数据为:

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

中间部分的数据就是卷积的有效数据。

仔细观察,可以发现这样的卷积在边缘处是有问题,他把超出边缘部分的数据全部用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变换是没有必要的,可以优化掉。

优化后的部分代码如下:

intIM_FFTConv2(unsignedchar*Src, unsignedchar*Dest,intWidth,intHeight,intStride,float*Kernel,intKerWidth,intKerHeight)//阈值{intChannel = Stride /Width;if((Src == NULL) || (Dest == NULL))returnIM_STATUS_NULLREFRENCE;if((Width <=0) || (Height <=0))returnIM_STATUS_INVALIDPARAMETER;if((KerWidth >50) || (KerHeight >50))returnIM_STATUS_INVALIDPARAMETER;//卷积核越大,每次的有效计算量就越小了if((Channel !=1) && (Channel !=3) && (Channel !=4))returnIM_STATUS_INVALIDPARAMETER;intStatus =IM_STATUS_OK;constintTileWidth =256, TileHeight =256;intHalfW = KerWidth /2, HalfH = KerHeight /2;intValidW = TileWidth +1- KerWidth - (KerWidth -1);//第一个加1是因为卷积扩展的矩阵大小为N+M-1,第二个-1是因为为了取得有效数据,还要对边缘进行扩展,扩展的大小为KerHeight - 1,比如KerHeight=5,则每边需要扩展2个像素,一共扩展4个像素intValidH = TileHeight +1- KerHeight - (KerHeight -1);//默认的卷积的效果再卷积周边是用0来填充的,如果分块处理时,这明显是不能满足要求的,会带来明显的块于块之间的分界线intTileAmountX = (Width / ValidW) + (Width % ValidW ?1:0);//图像需要分成的块数intTileAmountY = (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;gotoFreeMemory;
    }
    Status= IM_GetOffsetPos(RowOffset, Width, HalfW, TileAmountX * ValidW + (KerWidth - HalfW) - Width, IM_EDGE_MIRROR);//左右对称if(Status != IM_STATUS_OK)gotoFreeMemory;
    Status= IM_GetOffsetPos(ColOffset, Height, HalfH, TileAmountY * ValidH + (KerHeight - HalfH) -Height, IM_EDGE_MIRROR);if(Status != IM_STATUS_OK)gotoFreeMemory;

    memset(Conv,0, TileWidth * TileHeight *sizeof(Complex));//卷积核的其他元素都为0,这里先整体赋值为0for(intY =0; Y < KerHeight; Y++)
    {intIndex = Y *KerWidth;for(intX =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)gotoFreeMemory;if(Channel ==1)//单通道时可以一次性处理2个块{

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

    }elseif(Channel ==4)
    {
        Complex*TileBG = Tile, *TileRA = Tile + TileWidth *TileHeight;for(intTileY =0; TileY < TileAmountY; TileY++)
        {for(intTileX =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(intY = ValidR.Top; Y < ValidR.Bottom + KerHeight -1; Y++)
                {byte*LinePS = Src + ColOffset[Y] *Stride;intIndex = (Y - ValidR.Top) *TileWidth;for(intX = 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)gotoFreeMemory;

                Status= FFT2D(TileRA, TileRA, TileWidth, TileHeight,false,0, TileHeight - (ValidR.Bottom + KerHeight -1-ValidR.Top));if(Status != IM_STATUS_OK)gotoFreeMemory;for(intY =0; Y < TileWidth * TileHeight; Y++)
                {floatTemp =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)gotoFreeMemory;

                FFT2D(TileRA, TileRA, TileWidth, TileHeight,true);if(Status != IM_STATUS_OK)gotoFreeMemory;for(intY = ValidR.Top; Y < ValidR.Bottom; Y++)
                {byte*LinePD = Dest + Y * Stride + ValidR.Left *4;intIndex = (Y - ValidR.Top + KerHeight -1) * TileWidth + KerWidth -1;for(intX = 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);returnStatus;
}

程序里我固定每个块的大小为256*256,这主要是考虑256是4的整数次幂,是能够完全用SSE优化掉的,同时也占用更少的内存。

如果考虑更好的处理,在最后一列以及最后一行那个块,因为有效元素的减少,可以考虑使用更小的块来计算,不过这会增加程序的复杂性。本例没有考虑了。

速度测试:

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

卷积核15*15 348ms

卷积核25*25400ms

卷积核49*49600ms

以前我写过基于SSE直接卷积,同一个机器上也汇报下测试速度:

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

卷积核15*15 550ms

卷积核25*25 1300ms

所以我以前那篇文章的有些结论是错误的。






暗色背景
手机扫码阅读本页