Rational decasteljau

bezier_math.cpp

#include "bezier_math.h"

// References:
// https://en.wikipedia.org/wiki/B%C3%A9zier_curve#Rational_B%C3%A9zier_curves
// http://resources.mpi-inf.mpg.de/departments/d4/teaching/ss2012/geomod/slides_public/12_Rational_Splines.pdf
// https://doi.org/10.1016/B978-155860737-8/50013-2

namespace Bezier_Math{


float interpolate(float x, float y, float t){
    return t * y + (1.0f - t) * x;
}


// The modified DeCasteljau algorithm calculates points and normals on a rational bezier curve.
// The idea is as follows: As we know, the normal DeCasteljau linearly interpolates between points
// until only one remains. Now, to give more weight to one of these points, we can just
// change the interpolation method to one which favours that point.

//                          Diagrams: Weight of the 2nd point
//     |      .                                   |      .
//     |    .                                     |   .
//     |  .                                       | .
//     |.                                         |.
//     |-------                                   |-------
// Normal Decasteljau    Our modified Algorithm when the 2nd point has more weight

// The point with more weight influences the resulting interpolation point for longer.
// Crucially however, at t=0 the first point still holds all the influence, and at t=1 the second point does.
// This is regardless of weights.

// The normals are then calculated just like described in the hint.

void de_casteljau(const std::vector<Control_Point> &control_points, float t, glm::vec2 *point, glm::vec2 *normal){
    // We define n to be the length of the list of control points for efficiency
    unsigned long n = control_points.size();


    // If we only have one control point remaining, the coordinates of this point are the coordinates of the point we wanted to calculate.
    // We can this set the point and are done with the recursion.
    if (n == 1){
        *point = glm::vec2(control_points.at(0).x_, control_points.at(0).y_);
        return;
    }


    // If we only got 2 control points left, it's time to calculate the normal, like described in the hint.
    // For this, we take the vector between the 2 points, rotate and normalize it.
    if (n == 2){
        // vector coords
        float x = control_points.at(1).x_ - control_points.at(0).x_;
        float y = control_points.at(1).y_ - control_points.at(0).y_;

        // rotation
        float rotated_x = y;
        float rotated_y = -x;

        // determining normalization factor
        float length_before_normlize = (float)sqrt(rotated_x * rotated_x + rotated_y * rotated_y);

        // normalizing
        float normalized_x = rotated_x / length_before_normlize;
        float normalized_y = rotated_y / length_before_normlize;

        // set normal variable
        *normal = glm::vec2(normalized_x, normalized_y);
    }


    // If we made it to this point, we have more than 2 points remaining in the control_points list.
    // So, according to the deCasteljau Algorithm, we will interpolate between consecutive control points.
    // This will yield a list of new control points that is by 1 shorter than the list we got as an argument.
    // On this list, we can again run this method, until we run with just 1 point.

    std::vector<Control_Point> interpolated_points = std::vector<Control_Point>();

    // iterating through the list of n-1 consecutive control points
    for (unsigned long i = 1; i < n; i++){
        // We want to interpolate the control points at control_points[i-1] and control_points[i] and.
        // For convenience, we write the x,y position and weight of these points into variables.
        float x_0 = control_points.at(i - 1).x_;
        float y_0 = control_points.at(i - 1).y_;
        float x_1 = control_points.at(i).x_;
        float y_1 = control_points.at(i).y_;
        float w_0 = control_points.at(i - 1).w_;
        float w_1 = control_points.at(i).w_;

        // The weight of the new control point is just the linear interpolation of
        // the weights of the constituating control points.
        float new_w = interpolate(w_0, w_1, t);

        // The coordinates are interpolated, but not linearly, but in a way that takes their weights into account.
        //   Through (w_0 / new_w) * (1-t), we get the amount of influence p_0 should have on the point at position t,
        //   (w_1 / new_w) * t gives us the correct influence of p_1.
        //   Multiply that with the coordinate values, and you get the correct new coordinates at that point.
        //   This is what is being done here, but expressed in terms of linear interpolation.
        float new_x = interpolate((w_0 / new_w) * x_0, (w_1 / new_w) * x_1, t);
        float new_y = interpolate((w_0 / new_w) * y_0, (w_1 / new_w) * y_1, t);

        // Construct a new Control_Point from these values and add it to the list
        Control_Point newPoint = Control_Point(new_x, new_y);
        newPoint.setWeight(new_w);
        interpolated_points.push_back(newPoint);
    }

    // Run de_casteljau on the new list of interpolated points, which is shorter by 1.
    de_casteljau(interpolated_points, t, point, normal);

    return;
}


// The method uses a modified DeCasteljau type algorithm to calulate the currect correct curve normals
// and coordinates at num_points positions.
// For this, we first calculate the t interval between our sampling points.
// We then execute the modified DeCasteljau at multiples of this interval to get the point coordinates
//and normals, which we then write to their respective lists.
bool calculateBezierCurve(const std::vector<Control_Point> &control_points, const int &num_points,
                                       std::vector<glm::vec2> *curve_points, std::vector<glm::vec2> *normals){

    // Some simple sanity checks:
    if(control_points.size() < 2 || num_points < 2)
        return false;


    // Calculate the t-interval between sample points
    float interval = 1.0f / ((float) (num_points - 1));


    // The de_casteljau algorithm expects two pointers in which to write the point and normal results, thus we define those here
    glm::vec2* point = (glm::vec2*) malloc(sizeof(glm::vec2));
    glm::vec2* normal = (glm::vec2*) malloc(sizeof(glm::vec2));


    // Then we calculate point and normal positions at num_points multiples of the interval.
    for (int i = 0; i < num_points; i++){

        // The de_casteljau algorithm calculates the point and normal coordinates at the i'th multiple of the interval
        de_casteljau(control_points, interval * (float) i, point, normal);

        // The results are then pushed onto the list of results.
        curve_points->push_back(*point);
        normals->push_back(*normal);
    }

    // free the result vars before returning
    free(point);
    free(normal);

    // We are done here!
    return true;
}
}