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

#define MAX_POINTS 10

FILE* set_plot(int n,double p[][2])
{
    FILE* gp=popen("gnuplot -persistent","w");
    fprintf(gp,"set title 'bezier curve'\n");
    fprintf(gp,"set xrange[0:150]\n");
    fprintf(gp,"set yrange[0:150]\n");
    fprintf(gp,"set grid\n");

    fprintf(gp, "plot '-' with linespoints dt 2 pointsize 1.5 linewidth 2 linecolor rgb 'red' title 'Control Points', "
                    "'-' with linespoints pointtype 7 pointsize 1.5 linewidth 3 linecolor rgb 'blue' title 'Bezier Curve'\n");

    // Control points
    for(int i = 0; i <= n; i++) {
        fprintf(gp, "%f %f\n", p[i][0], p[i][1]);
    }
    fprintf(gp, "e\n");

    return gp;
}

void plot_curve(FILE* gp,double result[2])
{
    fprintf(gp, "%f %f\n", result[0], result[1]);
}

// De Casteljau Algorithm
void de_casteljau(double points[][2],int n, double t,double result[2])
{
    double temp[MAX_POINTS][2];

    // Copy original points
    for(int i = 0; i <= n; i++)
    {
        temp[i][0] = points[i][0];
        temp[i][1] = points[i][1];
    }

    // Iteratively interpolate
    for(int k = 1; k <= n; k++)
    {
        for(int i = 0; i <= n - k; i++)
        {
            temp[i][0] = (1 - t) * temp[i][0] + t * temp[i+1][0];
            temp[i][1] = (1 - t) * temp[i][1] + t * temp[i+1][1];
        }
    }

    // Final point
    result[0] = temp[0][0];
    result[1] = temp[0][1];
}

int main()
{
    int n = 3; // cubic Bezier (4 control points)

    double points[MAX_POINTS][2] = {
        {10, 10},
        {30, 120},
        {100, 120},
        {140, 20}
    };

    FILE* gp = set_plot(n, points);

    double result[2];

    // Generate curve points
    for(double t = 0.0; t <= 1.0; t += 0.01)
    {
        de_casteljau(points, n, t, result);
        plot_curve(gp, result);
    }

    fprintf(gp, "e\n"); // end Bezier curve data
    pclose(gp);

    return 0;
}