图像卷积的fft实现验证(python)

1 篇文章 0 订阅
1 篇文章 0 订阅

1、Caffe的卷积操作时间主要在矩阵乘法,假设一个m*n卷积核,且输入通道数为1,输出特征图大小为h*w,则乘法个数m*n*h*w,这里的优化仅限于对矩阵的乘法优化,因此,只要选择适合的矩阵计算库就可以了。

2、若使用FFT来计算图像卷积。其主要步骤如下。

假设输入图像的大小为len=h*w,卷积核大小k_len=m*n;通常len>>k_len;

  1. 对输入图像A做FFT,其算法的时间复杂度为o(lenlglen);

  2. 对卷积核B做FFT,但是由于卷积核与输入图像尺寸不一样,需要将卷积核扩展,即将卷积核倒置后,补len-k_len个0。

  3. 将A、B傅里叶变换的结果相乘,即对应位相乘获得结果C。乘法个数为len,时间复杂度为o(len)

  4. 对C做IFFT,得到结果D,在D中每隔k_len的值实部取出来,就是图像卷积的结果。因为图像卷积其实就是对应位相乘,所以需要每隔k_len取值,时间复杂度为o(len)

假设卷积核个数>1,需要对卷积核重新扩展后重复2)3)4)步骤,并与上一个卷积核图像卷积的值对应位相加就能获得。

验证正确性:

输入图像的卷积顺序-1,1,3,8 2,43,1,3 1,3,1,3 54,-2,-3,-1

卷积核1,2,-1,3

图像卷积结果:22,96,15,50.

此处注意,计算结果是用卷积核与输入图像进行乘法累加。

使用FFT结果:22,96,15,50。

此处注意,将卷积核逆置。

结果正确

 

 

暂时未能看出此方法的优越性,从空间上看,需要对卷积核进行扩展,其空间大小与输入图像的尺寸大小一样。时间上分析,仅二者FFT对应相乘的时间的乘法个数就和矩阵乘法个数的数量级是一样的了(当len>>k_len时)。适用于卷积核尺寸较大的情景。

图2出处:https://core.ac.uk/download/pdf/24989291.pdf

虽然没理解为什么性能提升这么多,但是该论文所做的实验证明了FFT性能很好,虽然复杂度推导跟我推导的差不多~此论文用了cuFFT库

 

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

至于多通道,多卷积核的,经python验证,卷积结果正确,其大概步骤如下:

水平有限,出错之处,希望得到各位的指正。

--------------------------------找到当时写的fft图像卷积,忘了对不对了---------------------------------------------------------------------------

template <typename Dtype>
void BaseConvolutionLayer<Dtype>::forward_cpu_gemm(const Dtype* input,
    const Dtype* weights, Dtype* output, bool skip_im2col) {
    const Dtype* col_buff = input;
    //get an array from a channel imag
    //get an array from a kernel,kernel_h*kernel_w=N;
    //make fft and ifft,res=ifft(),res[N,2N,3N....]is the channel conv result 
    //get another channel imag,and repeat
    //output=sum(conv result)
  const Dtype * data_im = input;
  const int channels = conv_in_channels_;
  const int height = conv_input_shape_.cpu_data()[1];
  const int width = conv_input_shape_.cpu_data()[2];
  const int kernel_h = kernel_shape_.cpu_data()[0];
  const int kernel_w = kernel_shape_.cpu_data()[1];
  const int pad_h = pad_.cpu_data()[0];
  const int pad_w = pad_.cpu_data()[1];
  const int stride_h = stride_.cpu_data()[0];
  const int stride_w = stride_.cpu_data()[1];
  const int dilation_h = dilation_.cpu_data()[0];
  const int dilation_w = dilation_.cpu_data()[1];
  const int output_h = (height + 2 * pad_h -
    (dilation_h * (kernel_h - 1) + 1)) / stride_h + 1;
  const int output_w = (width + 2 * pad_w -
    (dilation_w * (kernel_w - 1) + 1)) / stride_w + 1;
  const int channel_size = height * width;
  const int kernel_size = kernel_h * kernel_w;
  const int out_size = output_h * output_w;
  

    //only for test
   /*int channels=1;
   int channel_size=9;
   int out_size=4;
   int kernel_size=4;
   Dtype my_test[]={2,3,1,7,6,3,2,1,9};
   Dtype *my_test_p=my_test;
   Dtype my_kernel[]={4,2,1,0};
   int kernel_h=2;
   int kernel_w=2;
   int imag_w=3;*/
   int sum_k=0;
   int count_i=0;
    //make input
    //    std::cout<<"data_in_fft"<<*(input+channel_size)<<std::endl;
    memset(output,0,out_size);
    //std::cout<<"channels="<<channels<<std::endl;
    //for(int k_l = 0;k_l < conv_out_channels;k_l++;)
    //{
        for(int c = 0;c < channels ; c++)
        {
          //memset(data_in_fft,0,channel_size);
          //memset(data_kernel,0,channel_size);
          Dtype data_in_fft[out_size*kernel_size] = {0};
          Dtype data_kernel[out_size*kernel_size] = {0};
        for(int count=0;count<out_size;count++)
        {   //memcpy(data_in_fft+count*kernel_size,input+(channel_size*c)+count,kernel_size*sizeof(Dtype));
          if(count_i==width-kernel_w+1)
                        {
                            count_i=0;
                          data_im+=kernel_w-1;
                        }
          count_i++;
          for(int kk=0;kk<kernel_h;kk++)
                    {
            
                      memcpy(data_in_fft+sum_k,data_im+width*kk,sizeof(Dtype)*kernel_w);
            sum_k+=kernel_w;
            /*for(int jj=0;jj<out_size*kernel_size;jj++)
                       std::cout<<count+kk<<"th"<<*(data_in_fft+jj)<<std::endl;
            std::cout<<"--------------------------"<<std::endl;*/
          }
          data_im++;
        }
      /*for(int jj=0;jj<out_size*kernel_size;jj++)
                std::cout<<*(data_in_fft+jj)<<std::endl;
      std::cout<<"--------------------------"<<std::endl;*/
            cv::Mat a(1, out_size*kernel_size, CV_32F, data_in_fft);
            //memcpy(data_kernel,this->blobs_[0]->mutable_cpu_data()+(kernel_size*c),kernel_size*sizeof(Dtype));
      memcpy(data_kernel,this->blobs_[0]->mutable_cpu_data()+(kernel_size*c),kernel_size*sizeof(Dtype));
      //std::cout<<"data_kernel"<<*data_kernel<<std::endl;

            const Dtype* weight = this->blobs_[0]->cpu_data();

            for(int j = 0;j < (kernel_size>>1);j++)
            {
                int tmp;
                tmp=data_kernel[j];
                data_kernel[j]=data_kernel[kernel_size-j-1];
                data_kernel[kernel_size-j-1]=tmp;

            }   
            /*for(int jj=0;jj<out_size*kernel_size;jj++)
            std::cout<<*(data_kernel+jj)<<std::endl;
            std::cout<<"--------------------------"<<std::endl; */                                                             
            cv::Mat b(1,out_size*kernel_size,CV_32F,data_kernel);
            cv::Mat padded(a);
           
            cv::Mat planes[] = {cv::Mat_<float>(padded), cv::Mat::zeros(padded.size(), CV_32F)};
            cv::Mat complexImg;
            cv::merge(planes, 2, complexImg);
            cv::dft(complexImg,complexImg);
            cv::split(complexImg, planes);

            cv::Mat padded1(b);
            cv::Mat planes2[] = {cv::Mat_<float>(padded1), cv::Mat::zeros(padded1.size(), CV_32F)};
            cv::Mat complexImg2;
            cv::merge(planes2, 2, complexImg2);
            cv::dft(complexImg2,complexImg2);
            cv::split(complexImg2, planes2);
    

            cv::Mat planes3[] = {cv::Mat::zeros(padded1.size(),CV_32F), cv::Mat::zeros(padded1.size(), CV_32F)};
            planes3[0]=planes[0].mul(planes2[0])-planes[1].mul(planes2[1]);
            planes3[1]=planes[0].mul(planes2[1])+planes[1].mul(planes2[0]);

            cv::merge(planes3, 2, complexImg);
            cv::idft(complexImg,complexImg,CV_DXT_INV_SCALE);
            cv::split(complexImg, planes);
            //std::cout<<planes[0]<<std::endl;
            //std::cout<<"+++++++++++++++++++++++++++++++"<<std::endl;
            //std::cout<<planes[0]<<std::endl;
            //std::cout<<"*******************************"<<std::endl;
            //std::cout<<"output"<<std::endl;
            for(int m = kernel_size-1,i=0; m < out_size*kernel_size ;m += kernel_size)
            {
                 *(output +i)=planes[0].at<float>(0,m);
                 i++;
                //std::cout<<*(output+i)<<std::endl;
            }
            //for(int iii=0;iii<4;iii++)
            //std::cout<<*(output+iii)<<std::endl;
    
        }
    //}
    
}

 

  • 1
    点赞
  • 6
    收藏
    觉得还不错? 一键收藏
  • 11
    评论
评论 11
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值