RGB三通道图像二维卷积运算

方法一：分离通道，二维卷积 （处理大矩阵效率低）

1. 对图像进行padding，保证有效数据宽高和缩放后图像的宽高一致(卷积核中心点对应于实际坐标点）

int kernel_height = 7;
int kernel_width = 7;
int h_padding = kernel_height / 2;
int w_padding = kernel_width / 2;
int resize_height = 224;
int resize_width = 224;
int padding_height = resize_height + h_padding  * 2;
int padding_width = resize_width + w_padding * 2;
cv::Mat src = cv::imread(filename);
cv::Mat dst;
cv::Mat padding_dst(padding_height, padding_width, CV_32FC3, Scalar(0,0,0));
cv::resize(src, dst, cv::Size(resize_height, resize_width);
dst.convertTo(dst, CV_32FC3);
dst.copyTo(padding_dst(Rect(w_padding, h_padding, dst.cols, dst.rows)));

2. 基于opencv的二维卷积函数filter2D计算卷积

int kernel_shape[4] = {3,7,7,64};
vector<vector<Mat> > kernels;
//导入卷积核数据，构建卷积核矩阵，通道顺序为BGR
vector<Mat> channels;
split(padding_dst, channels);//BGR
int channel_size = padding_height * padding_width * kernel_shape[3];
float *filter_data = new float[3 * channel_size];
//对于每一个channel，针对每一个卷积核构建卷积结果矩阵
vector<vector<Mat> > result;
for(int i=0; i<kernel_shape[0]; i++){
    vector<Mat> temp;
    int start_offset = i*channel_size;
    for(int j=0; j<kernel_shape[3]; j++){
        int offset = start_offset + j*padding_height*padding_width;
        Mat mat(padding_height, padding_width, CV_32FC1, filter_data+offset);
    	temp.push_back(mat);
	}
	result.push_back(temp);
}
for(int i=0; i<kernel_shape[0]; i++){
    for(int j=0; j<kernel_shape[3]; j++){
        filter2D(channels[i], result[i][j], -1, kernels[i][j]);
    }    
}
for(int j=0; j<kernel_shape[3]; j++){
    for(int i=1; i<kernel_shape[0]; i++){
        result[0][j] += result[i][j];         
    }   
}

方法二：转换RGB图像和卷积核成二维矩阵，进行矩阵乘法（处理大矩阵效率高）

1. 利用im2col_cpu函数将RGB图像转换为二维矩阵（3*7*7 x 224*224），每一列对应一个像素点的RGB三通道卷积向量（3*7*7），通道顺序为RGB

int resize_height = 224;
int resize_width = 224;
cv::Mat src = cv::imread(filename);
cv::Mat dst;
cv::resize(src, dst, cv::Size(resize_height, resize_width);
dst.convertTo(dst, CV_32FC3);

int kernel_shape[4] = {3,7,7,64};
int stride_h = 1;
int stride_w = 1;
int rows = kernel_shape[0] * kernel_shape[1] * kernel_shape[2];
int cols = resize_height / stride_h * resize_width / stride_w;
float *convert_data = new float[rows * cols];
im2col_cpu((float *)dst.data, resize_height, resize_width, kernel_shape[1], kernel_shape[2], 
    kernel_shape[1]/2, kernel_shape[2]/2, stride_h, stride_w, 1, 1, convert_data);

template <typename Dtype>
void im2col_cpu(const Dtype* data_im,
                const int height, const int width,
                const int kernel_h, const int kernel_w,
                const int pad_h, const int pad_w,
                const int stride_h, const int stride_w,
                const int dilation_h, const int dilation_w,
                Dtype* data_col) 
{
    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 col_size = height * width * kernel_h * kernel_w / stride_h / stride_w;
    float* data_col_r = data_col;
    float* data_col_g = data_col + col_size;
    float* data_col_b = data_col + col_size * 2;
    for (int kernel_row = 0; kernel_row < kernel_h; kernel_row++) {
        for (int kernel_col = 0; kernel_col < kernel_w; kernel_col++) {
            int input_row = -pad_h + kernel_row * dilation_h;
            for (int output_rows = output_h; output_rows; output_rows--) {
                if (!is_a_ge_zero_and_a_lt_b(input_row, height)) {
                    for (int output_cols = output_w; output_cols; output_cols--) {
                        *(data_col_r++) = 0;
                        *(data_col_g++) = 0;
                        *(data_col_b++) = 0;
                    }
                }
                else {
                    int input_col = -pad_w + kernel_col * dilation_w;
                    for (int output_col = output_w; output_col; output_col--) {
                        if (is_a_ge_zero_and_a_lt_b(input_col, width)) {
                            int index = (input_row * width + input_col) * 3;
                            *(data_col_r++) = data_im[index + 0];
                            *(data_col_g++) = data_im[index + 1];
                            *(data_col_b++) = data_im[index + 2];
                        }
                        else {
                            *(data_col_r++) = 0;
                            *(data_col_g++) = 0;
                            *(data_col_b++) = 0;
                        }
                        input_col += stride_w;
                    }
                }
                input_row += stride_h;
            }
        }
    }
}

2. 初始化卷积核矩阵（64 x 3*7*7）的每一行为一个卷积核（3*7*7），通道顺序为RGB；利用openblas进行矩阵乘法，计算卷积结果

const int M = kernel_shape[3];//A的行数，C的行数
const int K = rows;//A的列数，B的行数
const int N = cols;//B的列数，C的列数
float *kernels = new float[M * K];
//初始化卷积核矩阵，通道顺序为RGB
float* filter_data = new float[rows * cols];
// MxK * KxN = MxN
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, M, N, K, 
     1, kernels, K, convert_data, N, 0, filter_data, N);