#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define BICUBIC_KERNEL_SIZE 4

double cubic(double x) {
    double absx = fabs(x);
    if (absx <= 1) {
        return 1 - 2 * pow(absx, 2) + pow(absx, 3);
    } else if (absx < 2) {
        return 4 - 8 * absx + 5 * pow(absx, 2) - pow(absx, 3);
    } else {
        return 0;
    }
}

double bicubic(double coeffs[BICUBIC_KERNEL_SIZE][BICUBIC_KERNEL_SIZE], double x, double y) {
    double result = 0;
    for (int i = 0; i < BICUBIC_KERNEL_SIZE; i++) {
        for (int j = 0; j < BICUBIC_KERNEL_SIZE; j++) {
            double cx = cubic(x - (double)i);
            double cy = cubic(y - (double)j);
            result += cx * cy * coeffs[i][j];
        }
    }
    return result;
}

void rescale_image(unsigned char *input_image, int input_width, int input_height, unsigned char *output_image, int output_width, int output_height) {
    double x_ratio = (double)input_width / (double)output_width;
    double y_ratio = (double)input_height / (double)output_height;
    double coeffs[BICUBIC_KERNEL_SIZE][BICUBIC_KERNEL_SIZE];
    int x, y;
    double x_diff, y_diff;
    int pixel_value;

    // Iterate over each pixel in the output image
    for (int i = 0; i < output_height; i++) {
        for (int j = 0; j < output_width; j++) {
            // Compute the coordinates in the input image
            double x_coord = ((double)j + 0.5) * x_ratio - 0.5;
            double y_coord = ((double)i + 0.5) * y_ratio - 0.5;

            // Compute the integer and fractional parts of the coordinates
            x = (int)x_coord;
            y = (int)y_coord;
            x_diff = x_coord - x;
            y_diff = y_coord - y;

            // Compute the bicubic coefficients
            for (int k = -1; k <= 2; k++) {
                for (int l = -1; l <= 2; l++) {
                    int row = y + k;
                    int col = x + l;

                    // Handle edge cases
                    if (row < 0) {
                        row = 0;
                    } else if (row >= input_height) {
                        row = input_height - 1;
                    }
                    if (col < 0) {
                        col = 0;
                    } else if (col >= input_width) {
                        col = input_width - 1;
                    }

                    // Get the pixel value at the current position
                    pixel_value = *(input_image + row * input_width + col);

                    // Compute the bicubic coefficient
                    coeffs[k+1][l+1] = (double)pixel_value;
                }
            }

            // Compute the new pixel value using bicubic interpolation
            pixel_value = (int)round(bic