模拟3维卷积的计算

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <ctime>
#include <time.h>
using namespace std;
ofstream out("output.txt");
// Initilize a 3D array
double*** initialize_3DArray( int x, int y, int z)
{
    double*** arr;
    arr = new double**[z];

    for(int i=0; i<z; i++)
    {
        arr[i] = new double*[x];
        for(int ii=0; ii<x; ii++)
        {
            arr[i][ii] = new double[y];
            for(int iii=0; iii<y; iii++)
            {
                arr[i][ii][iii] = (double)rand() / RAND_MAX;
            }
        }
    }
    return arr;
}

void free_3DArray(double*** arr, int x, int y, int z)
{
    for(int i = 0; i < z; i++)
    {
        for(int ii =0; ii < x; ii++)
        {
            delete [] arr[i][ii];
            arr[i][ii]=NULL;
        }
        delete [] arr[i];
        arr[i] = NULL;
    }
    delete [] arr;
    arr = NULL;
}

double** initialize_2DArray(int x, int y)
{
    double** arr = new double*[x];
    for(int i = 0; i < x; i++)
    {
        arr[i] = new double[y];
    }
    return arr;
}

void free_2DArray(double** arr, int x, int y)
{
    for(int i = 0; i < x; i++)
    {
        delete [] arr[i];
        arr[i] = NULL;
    }
    delete [] arr;
    arr = NULL;
}

// Purpose: perform the convolution between convoluted matrix and convolutor kernel, both them are 2D
// convoluted: the matrix beconvoluted, x0: size of its X-dimension, y0: size of its Y-dimension
// convolutor: the convolution kernel, x1: size of its X-dimension, y1: size of its Y-dimension
// the function return a convolution resulting matrix, whose x-Dimension: x0-x1+1, Y-Dimension:
// y0-y1+1.
double** convolution_2D(double** convoluted, int x0, int y0, double** convolutor, int x1, int y1)
{
    if(x0 < x1 || y0 < y1)
    {
        cout<<"Error: the kernel lager than convoluted matrix!"<<endl;
        exit(1);
    }

    int stepsize = 1;
    double** result;
    int result_x = x0-x1+1;
    int result_y = y0-y1+1;
    result = initialize_2DArray(result_x, result_y);

    for(int row = 0; row < result_x; row += stepsize)
    {
        for(int col = 0; col < result_y; col += stepsize)
        {
            double temp_sum = 0;

            // k_x: index of the X-dimesnion of the kernel, k_y: index of the Y-dimension of the kernel
            // row0: index of the X-dimesnion of the convoluted matrix, col0:index of the Y-dimension of the convoluted matrix
            int k_x, k_y, row0, col0;
            for(k_x = 0, row0 = row; k_x < x1, row0 < row+x1; k_x++, row0++)
            {
                for(k_y = 0, col0 = col; k_y < y1, col0 < col+y1; k_y++, col0++)
                {
                    temp_sum += convoluted[row0][col0] * convolutor[k_x][k_y];
                }
            }
            result[row][col] = temp_sum;
        }
    }

    return result;
}

// pupose: add two 2D matrix together, both them have size of x by y
double** add_2D(double** matrix1, double** matrix2, int x, int y)
{
    //double** tmp = initialize_2DArray(x,y);
    for(int i = 0; i < x; i++)
    {
        for (int ii = 0; ii < y; ii++)
        {
            matrix1[i][ii] = matrix1[i][ii] + matrix2[i][ii];
        }
    }
    return matrix1;
}

// In this version, we assume z0 = z1, but in general cases, the z0 can larger than z1.
double** convolution_3D(double*** convoluted, int x0, int y0, int z0, double*** convolutor, int x1, int y1, int z1)
{
    // store the temp result of convolution
    double** tmp_conv;
    int conv_result_x = x0-x1+1;
    int conv_result_y = y0-y1+1;
    double** conv_result = initialize_2DArray(conv_result_x, conv_result_y);
    // this loop summary the results of cube_z times convolution
    for(int recep = 0; recep < z0; recep++)
    {
        tmp_conv = convolution_2D(convoluted[recep],x0,y0,convolutor[recep],x1,y1);
        conv_result = add_2D(conv_result, tmp_conv, conv_result_x, conv_result_y);
    }
    return conv_result;
}

void show_2D(double** matrix, int x, int y)
{
    for(int i = 0; i < x; i++)
    {
        for (int ii = 0; ii < y; ii++)
        {
            cout<<matrix[i][ii]<<" ";
        }
        cout<<endl;
    }
}

void show_3D(double*** matrix, int x, int y, int z)
{
    for(int i = 0; i < x; i++)
    {
        for(int ii = 0; ii < x; ii++)
        {
            for (int iii = 0; iii < y; iii++)
            {
                cout<<matrix[i][ii][ii]<<" ";
            }
            cout<<endl;
        }
        cout<<endl;
    }
}

int main()
{
    //设置随即数生成器的种子
    srand( (unsigned)time(NULL) );

    // length of receptive field
    int recep_field = 3;
    // cube: it formed by stacked multiple frames(feature maps)
    double*** cube;
    /*
        sizes of three dimension of the cube
    */
    int cube_x;
    int cube_y;
    int cube_z;
    cube_x = 30;
    cube_y = 10;
    //cube_z = 7;
    cube_z = 3;
    cube = initialize_3DArray(cube_x,cube_y,cube_z);
    //free_3DArray(cube,cube_x,cube_y,cube_z);

    // Design the 3D-kernel
    double*** kernel;
    int kernel_x;
    int kernel_y;
    int kernel_z;
    kernel_z = 3;
    kernel_x = 7;
    kernel_y = 7;
    kernel = initialize_3DArray(kernel_x, kernel_y, kernel_z);
    //free_3DArray(kernel, kernel_x, kernel_y, kernel_z);

    // allocate space for convolutional result
    double** conv_result;
    int conv_result_x = cube_x-kernel_x+1;
    int conv_result_y = cube_y-kernel_y+1;
//    int conv_result_z = cube_z-kernel_z+1;
    //conv_result = initialize_2DArray(conv_result_x, conv_result_y);
    //free_2DArray(conv_result, conv_result_x, conv_result_y);

    /*****************************************************

                        3D convolution

    *****************************************************/
    conv_result = convolution_3D(cube,cube_x,cube_y,cube_z,kernel,kernel_x,kernel_y,kernel_z);
    show_2D(conv_result,conv_result_x,conv_result_y);
//    // store the temp result of convolution
//    double** tmp_conv;
//    // this loop summary the results of cube_z times convolution
//    for(int recep = 0; recep < cube_z; recep++)
//    {
//        tmp_conv = convolution_2D(cube[recep],cube_x,cube_y,kernel[recep],kernel_x,kernel_y);
//        conv_result = add_2D(conv_result, tmp_conv, conv_result_x, conv_result_y);
//    }

    free_3DArray(cube,cube_x,cube_y,cube_z);
    free_3DArray(kernel, kernel_x, kernel_y, kernel_z);

    delete [] conv_result;
    //free_2DArray(conv_result, conv_result_x, conv_result_y);
    return 0;
}