#include <iostream>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/core.hpp>
using namespace std;
using namespace cv;
/*
* 矩阵的掩码操作通俗解释：
*     给定一个小矩阵，把这个矩阵放置到要操作的大矩阵上
*     把两个矩阵重叠对应的像素值相乘，再加起来。通过移动
*     小矩阵，直到要做掩码操作的大矩阵所有像素位置都计算
*     完成。
* 
* 小矩阵叫做核kernel, 大小取奇数
* kernel = [ k00, k01, k02,
*            k10, k11, k12,
*            k20, k21, k22]
* 
* 被操作的矩阵为n行m列
* mat   = [ e(00), e(01), e(02), ..., e(0m),
*           e(10), e(11), e(12), ..., e(1m),
*           .
*           .
*           .
*          e(n0), e(n1), e(n2), ..., e(nm)]
*
* 进行掩码操作时, 一般会先按照核的大小扩展被操作的矩阵的大小
* 保证进行掩码操作前后被操作的矩阵的大小不变
* mat被做掩码操作之后第i行第j列元素为：
* a = i - 1, b = i + 1, c = j - 1, d = j + 1
* res(i, j) = mat_e(ac) * k00 + mat_e(aj) * k01 + mat_e(ad) * k02 + 
              mat_e(ic) * k10 + mat_e(ij) * k11 + mat_e(id) * k12 +
              mat_e(bc) * k20 + mat_e(bj) * k21 + mat_e(bd) * k22
* 如果有多个通道，每个通道分别计算
*/

static Mat getOpenCVGaussianKernel2D(int kSize = 3, double sigma = 0)
{
    if (!sigma)  // 标准差为0, 由核的大小来计算
    {
        int t = kSize - 1;
        sigma = 0.3 * (t * 0.5 - 1) + 0.8;
    }
    Mat kernel_x = getGaussianKernel(kSize, sigma, CV_64F);
    Mat kernel_y = getGaussianKernel(kSize, sigma, CV_64F).t(); // 一维的列向量是(3 * 1)的向量, 要把它转置成(1 * 3)的行向量
    Mat kernel = kernel_x * kernel_y;   
    return kernel;
}

static Mat getManualGaussianKernel2D(int kSize = 3, double sigma = 0)
{
    if (!sigma)  // 标准差为0, 由核的大小来计算
    {
        int t = kSize - 1;
        sigma = 0.3 * (t * 0.5 - 1) + 0.8;
    }

    // 二维高斯函数: G(x, y) = exp(-(x * x + y * y) / 2 * sigma * sigma)
    Mat kernel(kSize, kSize, CV_64FC1, Scalar::all(0));
    double denom = 2.0 * sigma * sigma;
    int center = kSize >> 1;
    double* ptr = nullptr;
    for (int row = 0; row < kSize; ++row)
    {
        ptr = kernel.ptr<double>(row);
        for (int col = 0; col < kSize; ++col)
        {
            double x = (double)col - center;
            double y = (double)row - center;
            ptr[col] = exp(-(x * x + y * y) / denom);
        }
    }

    // 归一化, 即每一项都除以所有项的和
    Scalar sum = cv::sum(kernel);
    kernel /= sum;
    return kernel;
}

// opencv提供了filter2D函数来进行掩码操作, 这里自己实现一下
static Mat applyGaussianKernel2D(const Mat& mat, const Mat& kernel)
{
    CV_Assert(mat.depth() == CV_8U && mat.channels() == 1);
    CV_Assert(kernel.rows == 3 && kernel.cols == 3 
              && kernel.depth() == CV_64F && kernel.channels() == 1);

    int rows = mat.rows, cols = mat.cols, type = mat.type();
    Mat mat_;
    // 对图像的边界进行扩展, 保证掩码操作前后被操作的矩阵的大小不变
    // 核的大小为3 * 3向四周扩展一个像素的宽度即可, 使用 BORDER_REFLECT_101 扩展方式
    copyMakeBorder(mat, mat_, 1, 1, 1, 1, BORDER_REFLECT_101);

    Mat res(rows, cols, type);
    const Mat_<double>& kernel_ = kernel;

    for (int row = 1; row < rows + 1; ++row)
    {
        const uchar* pPrevRow = mat_.ptr<uchar>(row - 1);
        const uchar* pCurrRow = mat_.ptr<uchar>(row);
        const uchar* pNextRow = mat_.ptr<uchar>(row + 1);
        uchar* ptr = res.ptr<uchar>(row - 1);

        for (int col = 1; col < cols + 1; ++col)
        {
            ptr[col - 1] = saturate_cast<uchar>(
                    kernel_(0, 0) * pPrevRow[col - 1] + kernel_(0, 1) * pPrevRow[col] + kernel_(0, 2) * pPrevRow[col + 1] +
                    kernel_(1, 0) * pCurrRow[col - 1] + kernel_(1, 1) * pCurrRow[col] + kernel_(1, 2) * pCurrRow[col + 1] +
                    kernel_(2, 0) * pNextRow[col - 1] + kernel_(2, 1) * pNextRow[col] + kernel_(2, 2) * pNextRow[col + 1]);
        }
    }

    return res;
}

int main()
{
    Mat mat(16, 16, CV_8UC1);   
    cv::randu(mat, Scalar::all(0), Scalar(255));
    cout << "src mat = \n" << format(mat, Formatter::FMT_PYTHON) << endl;

    int kSize = 3;
    double sigma = 0;
    Mat kernel1 = getOpenCVGaussianKernel2D(kSize, sigma);
    cout << "2D Gaussian kernel using getGaussianKernel() = " << endl
         << format(kernel1, Formatter::FMT_PYTHON) << endl;

    Mat kernel2 = getManualGaussianKernel2D(kSize, sigma);
    cout << "2D Gaussian kernel by manual calculatiton = " << endl
         << format(kernel2, Formatter::FMT_PYTHON) << endl;

    Mat kernelDiff = kernel1 - kernel2;
    if (kernelDiff.empty())
        return -1;
    cout << "kernel difference by two methods = " << endl
         << format(kernelDiff, Formatter::FMT_PYTHON) << endl;

    Mat matFilter2D;
    filter2D(mat, matFilter2D, mat.depth(), kernel1, Point(-1, -1), 0, BORDER_REFLECT_101);
    cout << "result mat by filter2D() = \n" << format(matFilter2D, Formatter::FMT_PYTHON) << endl;

    Mat matManualMask = applyGaussianKernel2D(mat, kernel2);
    cout << "result mat by manual = \n" << format(matManualMask, Formatter::FMT_PYTHON) << endl;

    Mat matDiff = matManualMask - matFilter2D;
    if (matDiff.empty())
        return -1;
    cout << "result matrixs difference by two methods = " << endl
        << format(matDiff, Formatter::FMT_PYTHON) << endl;

    return 0;
}