bicubic interpolation in C

bicubic.c

float cubic_interpolate(
    float p0, float p1, float p2, float p3, float ratio, float cubic_coeff_a) {

    float cca = cubic_coeff_a;

    float coeffects_0 = ((cca * (ratio + 1) - 5 * cca) * (ratio + 1) + 8 * cca) * (ratio + 1) - 4 * cca;
    float coeffects_1 = ((cca + 2) * ratio - (cca + 3)) * ratio * ratio + 1;
    float coeffects_2 = ((cca + 2) * (1 - ratio) - (cca + 3)) * (1 - ratio) * (1 - ratio) + 1;
    float coeffects_3 = ((cca * ((1 - ratio) + 1) - 5 * cca) * ((1 - ratio) + 1) + 8 * cca) * ((1 - ratio) + 1) - 4 * cca;

    return (p0 * coeffects_0 +
            p1 * coeffects_1 +
            p2 * coeffects_2 +
            p3 * coeffects_3 );
}



float cubic_interpolate_help(float *input, float ratio, float cubic_coeff_a, int idx_base, int row_offset, int highest_index) {
    int j0 = row_offset + CLAMP(0, highest_index, (int)idx_base + 0);
    int j1 = row_offset + CLAMP(0, highest_index, (int)idx_base + 1);
    int j2 = row_offset + CLAMP(0, highest_index, (int)idx_base + 2);
    int j3 = row_offset + CLAMP(0, highest_index, (int)idx_base + 3);

    float v0 = input[j0];
    float v1 = input[j1];
    float v2 = input[j2];
    float v3 = input[j3];

    float dxyz = cubic_interpolate(v0, v1,v2,v3, ratio, cubic_coeff_a);

    return dxyz;
}


void bicubic_interpolation(int32_t input_depth, int32_t input_height,
                             int32_t input_width, int32_t output_depth,
                             int32_t output_height, int32_t output_width,
                             float cubic_coeff_a,
                             float *scales,
                             float *input, float *output) {

    int length = input_width;
    int coeffects_length = 4;

    float scale_x = scales[2];
    float scale_y = scales[1];
    float scale_z = scales[0];

    int32_t z = 0;

    for (int32_t y = 0; y < output_height; ++y) {
        float origin_idx_y = ((float)y + 0.5) / scale_y - 0.5;
        float ratio_y = calc_ratio(origin_idx_y);
        int32_t idx_base_y  = calc_idx_base(origin_idx_y, 4);

    for (int32_t x = 0; x < output_width; ++x) {
        float output_length = scale_x * input_width;

        float origin_idx_x = ((float)x + 0.5) / scale_x - 0.5;
        float ratio_x = calc_ratio(origin_idx_x);
        int32_t idx_base_x  = calc_idx_base(origin_idx_x, 4);



        // offset for the "row"
        int oh0 = input_width * input_height * z + input_width * CLAMP(0, input_height - 1, idx_base_y + 0);
        int oh1 = input_width * input_height * z + input_width * CLAMP(0, input_height - 1, idx_base_y + 1);
        int oh2 = input_width * input_height * z + input_width * CLAMP(0, input_height - 1, idx_base_y + 2);
        int oh3 = input_width * input_height * z + input_width * CLAMP(0, input_height - 1, idx_base_y + 3);

        // TODO should this be input height or width?
        // final index of the current row
        int highest_index = input_width - 1;

        float d0 = cubic_interpolate_help(input, ratio_x, cubic_coeff_a, idx_base_x, oh0, highest_index);
        float d1 = cubic_interpolate_help(input, ratio_x, cubic_coeff_a, idx_base_x, oh1, highest_index);
        float d2 = cubic_interpolate_help(input, ratio_x, cubic_coeff_a, idx_base_x, oh2, highest_index);
        float d3 = cubic_interpolate_help(input, ratio_x, cubic_coeff_a, idx_base_x, oh3, highest_index);

        float dxyz = cubic_interpolate(d0, d1,d2,d3, ratio_y, cubic_coeff_a);

        // the first index on the current row
        int row_offset = output_width * output_height * z + output_width * y;
        int index = row_offset + x;
        output[index] = dxyz;
    }
    }
}