#include <stdio.h>
#include <time.h>
#include <arm_neon.h>

float inputImage[5][5] = {
    {1, 2, 3, 4, 5},
    {6, 7, 8, 9, 10},
    {11, 12, 13, 14, 15},
    {16, 17, 18, 19, 20},
    {21, 22, 23, 24, 25}
};

float kx[3] = { 1.0f / 4, 2.0f / 4, 1.0f / 4 };
float ky[3] = { 1.0f / 4, 2.0f / 4, 1.0f / 4 };
float buf[3][3] = { 0 };

void applySeparableGaussianBlurNEON(float src[][5], float dst[][5], float kx[], float ky[]) {
    float32x4_t kx_vec = vld1q_f32(kx);
    float32x4_t ky_vec = vld1q_f32(ky);

    for (int i = 1; i < 3; ++i) {
        for (int j = 0; j < 3; j += 4) {
            float32x4_t left = vld1q_f32(&src[i - 1][j + 1]);
            float32x4_t center = vld1q_f32(&src[i][j + 1]);
            float32x4_t right = vld1q_f32(&src[i + 1][j + 1]);

            float32x4_t result_left = vmulq_f32(left, kx_vec);
            float32x4_t result_center = vmulq_f32(center, kx_vec);
            float32x4_t result_right = vmulq_f32(right, kx_vec);

            float32x4_t total_result = vaddq_f32(result_left, result_center);
            total_result = vaddq_f32(total_result, result_right);

            vst1q_f32(&buf[i][j], total_result);
        }
    }

    for (int i = 0; i < 3; ++i) {
        for (int j = 1; j < 3; j++) {
            float32x4_t buf_data = vld1q_f32(&buf[i][j - 1]);

            float32x4_t result = vmulq_f32(buf_data, ky_vec);

            vst1q_f32(&dst[i + 1][j + 1], result);
        }
    }
}

int main() {
    float outputImage[5][5] = {0}; 
    clock_t start, end;

    start = clock();
    applySeparableGaussianBlurNEON(inputImage, outputImage, kx, ky);
    end = clock();

    double timeSpent = (double)(end - start) / CLOCKS_PER_SEC;

    for (int i = 0; i < 5; i++) {
        for (int j = 0; j < 5; j++) {
            printf("%.6f ", outputImage[i][j]);
        }
        printf("\n");
    }

    printf("Time spent: %.6f seconds\n", timeSpent);
    return 0;
}