Last active
August 15, 2022 15:00
-
-
Save finia2NA/82fb32e9fd58f23121bd8cf576732b31 to your computer and use it in GitHub Desktop.
Rational decasteljau
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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; | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment