//C is the spline points container, x,y

void draw_spline(std::vector<std::vector<double>> C, double red, double green, double blue, double alpha){
    std::vector<POINT>::iterator it;
    //the forward and backward function has only to be recalculated when spline controlpoints are edited by user.
    std::vector<POINT> forward=draw_forward_spline(C,red,green,blue,alpha);
    std::reverse(C.begin(),C.end());
    std::vector<POINT> backward=draw_backward_spline(C,red,green,blue,alpha);
    std::reverse(backward.begin(),backward.end());
    std::vector<POINT> differce;
    POINT p;

    for(int i=0; i<forward.size(); i++){
        p.x=(forward.at(i).x+backward.at(i).x)/2;
        p.y=(forward.at(i).y+backward.at(i).y)/2;
        differce.push_back(p);
    }

    glLineWidth(1);
    glColor4d(red,green,blue,alpha);
    glBegin(GL_LINE_STRIP);
    for(it=differce.begin(); it<differce.end(); it++){
            glVertex2d(it->x,it->y);
    }
    glEnd();
}

std::vector<POINT> draw_forward_spline(std::vector<std::vector<double>> C, double red, double green, double blue, double alpha){
    std::vector<POINT> points;
    double S[C.size()-1][2][4];
    double w, b[C.size()], d[C.size()], x[C.size()];
    int i, dim;
    int n = C.size()-1; // number of splines

    for(dim=0; dim<2; dim++)
    {
        // define d[]:
        d[0]=3.0*(C[1][dim]-C[0][dim]);
        for(i=1; i<n; i++)
            d[i]=3.0*(C[i+1][dim]-C[i-1][dim]);
        d[n]=3.0*(C[n][dim]-C[n-1][dim]);
        // forward sweep: (simplified tridiagonal solution: all a[i]=1 and all c[i]=1)
        b[0]=2.0;
        for(i=1; i<C.size(); i++)
        {
            w = 1.0/b[i-1];
            b[i] = 4.0 - w;
            d[i] -= (w*d[i-1]);
        }
        // calculate solution vector x[i] = D[i]:
        // (Wikipedia x[] = Wolfram D[])
        x[n]=d[n]/b[n];
        for(i=n-1; i>=0; i--)
            x[i]=(d[i]-x[i+1])/b[i];
        // calculate spline S[i][dim] a, b, c and d:
        for(i=0; i<n; i++)
        {
            S[i][dim][0]=C[i][dim];
            S[i][dim][1]=x[i];
            S[i][dim][2]=3.0*(C[i+1][dim]-C[i][dim]) - 2.0*x[i] - x[i+1];
            S[i][dim][3]=2.0*(C[i][dim]-C[i+1][dim]) + x[i] + x[i+1];
        }
    }


    for(int p=0; p<C.size()-1; p++)  //spline points
    {
        for(double t=0; t<1; t+=0.1)  //time 0-1 for each spline
        {
            POINT point;
            point.x = ((S[p][0][3]*t+S[p][0][2])*t+S[p][0][1])*t+S[p][0][0];
            point.y = ((S[p][1][3]*t+S[p][1][2])*t+S[p][1][1])*t+S[p][1][0];
            //std::cout<< "x:" << px << " y:" << py << std::endl;
            points.push_back(point);
        }
    }
    return points;
}

std::vector<POINT> draw_backward_spline(std::vector<std::vector<double>> C, double red, double green, double blue, double alpha){
    std::vector<POINT> points;
    double S[C.size()-1][2][4];
    double w, b[C.size()], d[C.size()], x[C.size()];
    int i, dim;
    int n = C.size()-1; // number of splines

    for(dim=0; dim<2; dim++)
    {
        // define d[]:
        d[0]=3.0*(C[1][dim]-C[0][dim]);
        for(i=1; i<n; i++)
            d[i]=3.0*(C[i+1][dim]-C[i-1][dim]);
        d[n]=3.0*(C[n][dim]-C[n-1][dim]);
        // forward sweep: (simplified tridiagonal solution: all a[i]=1 and all c[i]=1)
        b[0]=2.0;
        for(i=1; i<C.size(); i++)
        {
            w = 1.0/b[i-1];
            b[i] = 4.0 - w;
            d[i] -= (w*d[i-1]);
        }
        // calculate solution vector x[i] = D[i]:
        // (Wikipedia x[] = Wolfram D[])
        x[n]=d[n]/b[n];
        for(i=n-1; i>=0; i--)
            x[i]=(d[i]-x[i+1])/b[i];
        // calculate spline S[i][dim] a, b, c and d:
        for(i=0; i<n; i++)
        {
            S[i][dim][0]=C[i][dim];
            S[i][dim][1]=x[i];
            S[i][dim][2]=3.0*(C[i+1][dim]-C[i][dim]) - 2.0*x[i] - x[i+1];
            S[i][dim][3]=2.0*(C[i][dim]-C[i+1][dim]) + x[i] + x[i+1];
        }
    }

    for(int p=0; p<C.size()-1; p++)  //spline points
    {
        for(double t=0; t<1; t+=0.1)  //time 0-1 for each spline
        {
            POINT point;
            point.x = ((S[p][0][3]*t+S[p][0][2])*t+S[p][0][1])*t+S[p][0][0];
            point.y = ((S[p][1][3]*t+S[p][1][2])*t+S[p][1][1])*t+S[p][1][0];
            //std::cout<< "x:" << px << " y:" << py << std::endl;
            points.push_back(point);
        }
    }
    return points;
}