blur

#include "blur.h"

#include <valarray>
#include <chrono>
#include <iostream>


namespace {
    double CalculateGaussian(size_t delta_x, size_t delta_y, double sigma) {
        double coefficient = sigma * sigma * 2 * std::numbers::pi_v<double>;
        double degree = (delta_x * delta_x + delta_y * delta_y) / (2 * sigma * sigma);
        double result = exp(-1 * degree) / coefficient;
        return result;
    }

    NormalChannel CalculateGaussianMatrix(double radius, double sigma) {
        auto start = std::chrono::high_resolution_clock::now();
        NormalChannel matrix(2 * radius, std::vector<double>(2 * radius, 0));
        for (size_t x = 0; x < 2 * radius; ++x) {
            for (size_t y = 0; y < 2 * radius; ++y) {
                matrix[y][x] = CalculateGaussian(x, y, sigma);
            }
        }
        auto stop = std::chrono::high_resolution_clock::now();
        auto execution_time = duration_cast<std::chrono::milliseconds>(stop - start).count();

        std::cout << "Matrix built in " << std::round(execution_time / 10) / 100 << " seconds" << std::endl;
        return matrix;
    }

    NormalRGBMatrix ApplyGaussian(NormalRGBMatrix &matrix, size_t height, size_t width, double sigma) {
        auto new_matrix = matrix;
        static const double RADIUS = 20;
        const auto gauss_matrix = CalculateGaussianMatrix(RADIUS, sigma);

        double accumulation_red;
        double accumulation_green;
        double accumulation_blue;

        for (size_t x0 = 0; x0 < width; ++x0) {
            for (size_t y0 = 0; y0 < height; ++y0) {
                accumulation_red = 0;
                accumulation_green = 0;
                accumulation_blue = 0;

                size_t left_border = x0 < RADIUS ? 0 : x0 - RADIUS;
                size_t top_border = y0 < RADIUS ? 0 : y0 - RADIUS;

                if (x0 + RADIUS >= width) {
                    left_border = std::max<double>(width - 2 * RADIUS, 0);
                }

                if (y0 + RADIUS >= height) {
                    top_border = std::max<double>(width - 2 * RADIUS, 0);
                }

                size_t right_border = std::min<double>(left_border + 2 * RADIUS, width);
                size_t bottom_border = std::min<double>(top_border + 2 * RADIUS, height);

                for (size_t x = left_border; x < right_border; ++x) {
                    for (size_t y = top_border; y < bottom_border; ++y) {
                        double delta_x = x0 > x ? x0 - x : x - x0;
                        double delta_y = y0 > y ? y0 - y : y - y0;

                        if (delta_x > RADIUS || delta_y > RADIUS) {
                            continue;
                        }

                        auto gaussian = gauss_matrix[delta_x][delta_y];

                        accumulation_red += matrix.red_chanel[y][x] * gaussian;
                        accumulation_green += matrix.green_chanel[y][x] * gaussian;
                        accumulation_blue += matrix.blue_chanel[y][x] * gaussian;
                    }
                }
                new_matrix.red_chanel[y0][x0] = std::max<double>(0, std::min<double>(accumulation_red, 1));
                new_matrix.green_chanel[y0][x0] = std::max<double>(0, std::min<double>(accumulation_green, 1));
                new_matrix.blue_chanel[y0][x0] = std::max<double>(0, std::min<double>(accumulation_blue, 1));
            }
        }
        return new_matrix;
    }
}


void Blur::Process(Image &input) {
    auto rgb = input.GetNormalChannels();

    auto new_rgb = ApplyGaussian(rgb, input.height, input.width, sigma_);

    input = Image(new_rgb);
}