tangent and curvature

[pietroblandizzi]

Hello i am trying to use the library to fit splines on some data. It looks like the derivatives are not working properly. Maybe i am using the library in the wrong way. From what i can see the time is not exactly related to the x of each point. I compute the first derivative at time t : auto d1 = ucbs.getTangent(u).tangent; I compute the second derivative at time t : auto d2 = ucbs.getCurvature(u).curvature;

then i extract x and y and plot them. The data are really simple:
(x,y); (1, 30); (2, 29); (3, 27); (4, 24); (5, 20); (6, 16); (7, 9);

here the curve fitted with an uniform bspline and the two derivatives from the library:

can you help me? Am i using the lib in the wrong way?

[ejmahler]

Without digging too deep into the numbers, I think this is working as intended. The tangent is the rate of change of the spline -- but it's the rate of change with respect to T, not with respect to X.

Your X values increase by 1 for each point, and this library gives one unit of T for each input point, so I would expect the X component of the tangent to always be 1. Eyeballing the graph, it looks like all the points have X=1

The curvature is the rate of change of the tangent, again with respect to T. Because the tangent is 1 at every data point, it is not changing, and its rate of change is 0. And eyeballing the curvature graph, all the points have X=0

[ejmahler]

It sounds like you have an independent variable and a dependent variable, and you want to use a spline to relate them? Right now the library is only supporting cases where the independent variable is a fixed interval, but it wouldn't be too hard to add a constructor to the spline type you want to use, which takes T values and their corresponding data points

[pietroblandizzi]

Yes actually you are completely right. I am in the case you mentioned. Thank you for the answer. How could i do this for an uniform B-Spline? Could you help me? I would like to obtain d'(x) and d''(x) Thank you again.

[pietroblandizzi]

Hello, can you please tell me which of the classes i should modify in order to handle my case of independent variable not at fixed interval? thank you.

[ejmahler]

You can use UniformCubicBSpline<float>

Note that the T values will start from 0, instead of from 1

[pietroblandizzi]

But if i have point in x y and my x data are not equally distributed and not equally spaced how can i relate those points with t? because then i would like to evaluate the spline as f(x) f'(x) and f''(x) Is there an easy way to implement it in the library? Thank you

[ejmahler]

Ah, I missed that your X values weren't evenly spaced.

The is is an interesting problem. There isn't really a mathematically correct way to solve this.

One possible idea is having two parallel UniformCubicBSplines - one containing your X values and one containing your Y values.

UniformCubicBSpline<float> xValues = ....; UniformCubicBSpline<float> yValues = ....;

Then, in order to look up the y value at X=1.25f, create a SplineInverter

SplineInverter<float, float, 1> xInverter(xValues);

float t = xInverter.findClosestT(1.25f);

float yResult = yValues.getCurvature(t);

[ejmahler]

Another option is to add a constructor to GenericBSpline: https://github.com/ejmahler/SplineLibrary/blob/master/spline_library/splines/generic_b_spline.h#L189

The existing constructor takes a list of points and auto-computes the knots. The new constructor could take both the points and knots. It could look something like this:

GenericBSpline(std::vector<floating_t> knots, const std::vector<InterpolationType> &points, size_t degree)
        :SplineImpl<GenericBSplineCommon, InterpolationType,floating_t>(points, points.size() - degree)
    {
        assert(points.size() > degree);
        assert(knots.size() == points.size() + degree - 1);
        assert(knots[degree - 1] == 0.0f) // Make this a fuzzy compare, and put "knots[degree - 1] = 0.0f" on the next line?

        for(size_t i = 1; i < knots.size(); i++)
        {
              assert(knots[i] >= knots[i - 1]);
        }

        this->common = GenericBSplineCommon<InterpolationType, floating_t>(points, std::move(knots), degree);
}

Note that you actually need more knots than you do points. Specificially, you need an extra degree - 1 knots, added to the beginning so that knots[degree - 1] == 0. This is just a consequence of the Deboor algorithm used to compute the Generic B Spline.

Since we're talking about cubic splines, I would start by setting degree = 3. So you need to insert 3 - 1 = 2 "dummy knots" at the beginning of the knot array. There's no mathematically correct way to choose these, but here are some options:

• Repeat 0.0f 3 times. IE, knots[0] = 0 and knots[1] = 0 and knots[2] = 0.
• Increment by 1, so knots[0] = -2 and knots[1] = -1andknots[2] = 0`.
• Compute the average difference between X values, and increment by that. If the average difference is 0.6, then knots[0] = -1.2 and knots[1] = -0.6andknots[2] = 0`

Once you do all of this, then you can construct a GenericBSpline with your X values as knots, and your Y values as points:

GenericBSpline<float> mySpline(xValues_WithDummyValues, yValues, 3);
float y = mySpline.getPosition(1.5f);

If I had to guess, I'd say this approach is goingto be easier than the double-spline suggestion I made above, so start with this. Good luck :)

[pietroblandizzi]

I am going to try. There are problems if you try to use float or double as InterpolationType due to some functions that assume arrays like this:

auto segmentFunction = [this, index](floating_t t) -> floating_t { auto tangent = computeTangent(index, t); tangent.length(); };

I tried to use those formula: http://www.qc.edu.hk/math/Certificate%20Level/Parametric%20differentiation.htm

but i always get the second derivative scaled of some factor. I dont know the reason. Do you know id those are suitable for the problem? I mean computing the t for a certain x then using the getTangent and the get Curvature and applying those formula?

[pietroblandizzi]

Maybe i am getting something wrong but i think you cant use as InterpolationType something that is not an array unless changing the templates and also the SplineInverter plus some other utils classes. Am i wrong? Do you have an easy fix to use your proposed solution of the 2 splines?

[pietroblandizzi]

I created the Generic spline with your idea. The position work but if i plot the derivatives it is completely wrong. Do you have an idea of the reason?

Here the plot of the arc fitted first and second derivative

each plot is representing getPosition(x) x,y getTangent(x).tangent x,y getCurvature(x).curvature x,y

If you ask for a certain x the getPosition(x),x() is different. Moreover is not interpolating from the beginning. This knot[2] = 0 is making an huge shift in the all fitting. Do you have an idea?

[ejmahler]

• You're right, you can't use float. The Vector class has a template parameter for number of dimensions. So you can use GenericBSpline<Vector<1>>

• Can you post your code? I tried a setup of this and it seemed to produce meaningful values. I used the exact constructor I posted, and I used this code to set up my data:

template<size_t N>
 static std::vector<Vector<N>> generateTriangleNumberData(size_t size) {
     std::vector<Vector<N>> result(size);
     size_t currentTriangleNumber = 0;
     for(size_t i = 0; i < size; i++)
     {
         currentTriangleNumber += i;

         Vector<N> nextPosition{};
         for(size_t c = 0; c < N; c++) {
             nextPosition[c] = currentTriangleNumber;
         }

         result[i] = nextPosition;
     }
     return result;
 }

void makeSpline(void)
{
    float degree = 3;
    auto knots = std::vector<float>{-2,-1,0,0.75,1,5,3,4,5,6,7};
    auto points = generateTriangleNumberData<1>(knots.size() - (degree - 1));

    GenericBSpline<Vector<1>> spline(knots, points, 3);

    float interval = 0.5f;

    for(size_t i = 0; i < 10; i++) {
        auto computedData = spline.getCurvature(i * interval);
    }
}

• For the offset, you could try adding one padding value at the beginning, and another padding value at the end. So knots[1] is 0, knots[0] is -1, and then there's an extra knot on the end.

[pietroblandizzi]

` void GenericSpline::Fit(const std::vector& X, const std::vector& Y) { assert(X.size() == Y.size()); assert(X.size() > 2);

std::vector knots; knots.push_back(-2.0); knots.push_back(-1.0); knots.push_back(0.0);

data_.push_back(QVector2D(X.at(0), Y.at(0))); for(auto i = 0; i < X.size(); ++i) { knots.pushback(X.at(i)); data.push_back(QVector2D(X.at(i), Y.at(i))); } ptoSplineT.reset(new GenericBSpline<QVector2D, double>(knots, data_, size_t(3))); auto last = X.at(X.size()-1); for (auto i = 0.0; i < last; i += 0.25 ) auto p = ptoSplineT->getPosition(x); auto dx1 = ptoSplineT->getTangent(x).tangent; auto dx2 = ptoSplineT->getCurvature(x).curvature; }` So the two vectors contain the X and Y data Maybe i create the knots in the wrong way. But when i plot p.x(), p.y()
d1x.x(), d1x.y() d2x.x(), d2x.y(); I get the plot i have shown you above

From where it come from this auto knots = std::vector{-2,-1,0,0.75,1,5,3,4,5,6,7};

and why do you use a variable interval?

Can you please help me? Thank you

[ejmahler]

The problem could be that you're treating your X values as dependent variables.

Instead of plotting these values:

auto p = p_toSplineT_->getPosition(i);
auto dx1 = p_toSplineT_->getTangent(i).tangent;
auto dx2 = p_toSplineT_->getCurvature(i).curvature;

//plot:
p.x(), p.y()
d1x.x(), d1x.y()
d2x.x(), d2x.y();

What happens if you plot these:

auto p = p_toSplineT_->getPosition(i);
auto dx1 = p_toSplineT_->getTangent(i).tangent;
auto dx2 = p_toSplineT_->getCurvature(i).curvature;

//plot:
i, p.y()
i, d1x.y()
i, d2x.y();

Think of the input into the spline as your X value (in this case, the variable i is the input that we pass to getPosition() etc), and the output of the spline as the Y value. Does that improve your results?

[ejmahler]

From where it come from this auto knots = std::vector{-2,-1,0,0.75,1,5,3,4,5,6,7};

They're just arbitrary numbers I chose for my x values. Let me rename some variablesto make it more similar to your setup and see if that helps:

void makeSpline(void)
{
    float degree = 3;
    std::vector<float> X = std::vector<float>{-2,-1,0,0.75,1,5,3,4,5,6,7};
    std::vector<Vector<1>> Y = generateTriangleNumberData<1>(knots.size() - (degree - 1));

    GenericBSpline<Vector<1>> spline(X, Y, 3);

    float interval = 0.5f;

    for(size_t i = 0; i < 10; i++) {
        float x = i * interval;
        auto computedData = spline.getCurvature(x);
        float y = computedData.position[0];
        float dy = computedData.tangent[0];
        float d2y = computedData.curvature[0];
    }
}

[ejmahler]

and why do you use a variable interval?

It's just a different way of writing the same thing that you wrote.

This:

for(float x = 0.0; x < max; x += 0.5) {
    //loop code
}

is the same as this:

float interval = 0.5;
for(int i = 0; i * interval < max; i++) {
    float x = i * interval;

    //loop code
}

But I prefer the second one because there are fewer floating point errors

[pietroblandizzi]

Thank you. I will try but just a question: The X and Y size should be the same or is it ok that Y[0] will correspond to X[0] which is -2?? Same for the other values. This create a shift