How do you fit a curve through 3D points, using splines, constrained within a volume.

Suppose I have an array of points in 3space,

[r1; r2; ... ; rn]

where

ri = [xi yi zi]

are it's coordinates. I want to fit a curve (a spline) in 3 dimensions between these points. However, I want the interpolated curve between points ri and rj to be constrained within a certain volume. How could I go about doing this?

EDIT: To elaborate, each pair of points exist on distinct faces of a shared tetrahedral. The interpolated curve between said points has to stay within this tetrahedral.

[Matt J]

You'll need to clarify what you mean by "between points ri and rj". This is 3D, so in what way can 2 points bound a 3D region? Also, please elaborate on the shape of the volume, e.g., whether convex, as Sean was asking.

[Laurence hutton-smith]

Thanks for your comment Matt, each pair of points exist on separate faces of a shared tetrahedral. The interpolated curve between them has to stay within this tetrahedral.

[Matt J]

The interpolated curve

Originally, you said you were fitting, but now you say you're interpolating. Originally, also, you said you were fitting a line, but now you say it's a curve. We need to know definitively which it's to be.

An interpolation of the points, for example, could be done just be connecting each pair of points sharing a given tetrahedron with a line segment. Do that for all points and you get a curve that definitely stays inside the union of all the tetahedra.

[Laurence hutton-smith]

Sorry for not being too specific, I said 'fitting through' which I thought implied the use of control points, also the line I was refering to was a curve. I do not have a background in numerical methods, so my terminology is a bit off.

I realise now I need to be generating splines. That is currently what is happening, i.e. a line connecting the points, I was hoping for a higher degree function to connect the two. For example be able to connect all the points with a cubic spline, but have each spline segment constrained within a volume.

[Matt J]

I see. Well, you definitely cannot do what you are attempting with a cubic spline if the spline has to actually pass through the points and satisfy all the usual constraints of spline fitting (e.g., continuous derivatives at the knots). All degrees of freedom on the spline coefficients are eaten up just by satisfying those usual constraints, and you are left with no degrees of freedom with which to satisfy your tetrahedral constraints.

[Laurence hutton-smith]

I see. So what sort of spline or method of interpolation might be able to produce such a result?

[Matt J]

You can liberate some degrees of freedom by relaxing requirements on the splines. For example, do you need the curve to have the usual twice differentiability of cubic splines or is once-differentiability enough? What is insufficient about connecting the points with line segments? Is it simply that you need something smooth?

I think it is a tough problem, regardless.

[John D'Errico]

Agreed, IF the interpolant is to be a smooth one.

[Laurence hutton-smith]

Once differentiability would be enough, anything above linear is an improvement on the current scheme. This is modelling a trajectory of sorts, linear approximations to paths always underestimates the length of the trajectory, a smooth curve would give more accurate statistics in later analysis.

[Matt J]

Note that even once-differentiability could be impossible, depending on the locations of your ri points. For example, if one of the points lies at the corner of one of the tetrahedra, it will be impossible to pass a differentiable curve through it that doesn't exit the constrained region of the tetrahedra.

Answer 1

[John D'Errico]

Ugh, enough to say about this that it seems worth an answer instead of just a comment on Sean's answer.

It sounds like you have a facet of a tetrahedron, and a set of points that lie inside that facet, and are in the plane of the facet. Now, you want to find a spline that passes through the points, and is ensured to lie entirely inside the facet, as well as staying entirely in-plane.

The obvious answer is you can project the points into a 2-dimensional subspace, one that is defined by the facet. Then one could do any fits needed in the planar subspace, and the problem is completely a 2-d problem. At the end, you can reverse the projection, returning the interpolated points back into 3-d.

In fact, I think this will not be necessary, since a spline interpolant will result in a linear combination of the points. If the points are all in-plane, then any linear combination of those points will also be in-plane too. (I'll think about this claim, but I'm pretty confident it will hold.)

The interpolation itself can be done using tools like pchip, my own interparc, etc. A problem is how to enforce that the points lie inside a polygonal domain as defined by that facet. The suggestion of vert2con will fail here however. Why?

vert2con will produce a set of linear constraints that could in theory be applied to the spline. But the spline is a nonlinear function, even though it is linear in the spline parameters. Those linear inequalities cannot be applied to the spline, as a FUNCTION. You could sample the spline at a gazillion (the technical term for a LOT) points, then apply those inequalities to each point. Essentially this will be a NECESSARY constraint on the spline, that it lies inside the volume, but it will not be a sufficient condition. This is the same as how monotonicity constraints can be applied to a spline. If you constraint the spline to be monotone at a finite set of points, then you will have a necessary, but not sufficient condition that the spline be monotonic as a function. So there could be a small region (or regions) where the spline actually fails to satisfy the constraints.

And of course, you need to be careful. For example, suppose the points happen to lie at the vertices of the triangular facet? For example, suppose you choose the points as rows of xy below:

xy = [1 0;1 1;0 1];

They define a triangular domain. Can you now find a spline that lies entirely inside that triangular region, AND also interpolate those points? Of course not, if the spline interpolant will be at least a C1 function. No matter how you define a smooth interpolant through those points, if it is differentiable at the vertices of the triangle, then it must go outside of the triangle at that second vertex.

xyt = [0 0;0 1;1 0];
xy = rand(100,2);
ind = sum(xy,2)<1;
xy = xy(ind,:);

plot(xyt(:,1),xyt(:,2),'ob-')
hold on
plot(xy(:,1),xy(:,2),'g+',xyint(:,1),xyint(:,2),'r:')

As you can see, the interpolated curve falls entirely inside the triangular domain, mainly because I used a pchip interpolant for the purpose. That helps, but it won't ensure the curve always does. I'm quite confident that with only a tiny amount of effort, I can come up with an example where that fails. For example...

xy = [0.1 0.1;0.3 0.69;0.69 0.3;0.2 0.2];
xyt = [0 0;0 1;1 0];
xyint = interparc(1000,xy(:,1),xy(:,2),'pchip');
plot(xyt(:,1),xyt(:,2),'ob-')
hold on
plot(xy(:,1),xy(:,2),'g+',xyint(:,1),xyint(:,2),'r:')

I could easily enough show you how to do the interpolation using pchip directly, without need to use interparc here. Regardless, the point is that ensuring the curve falls entirely inside a polygonal domain is difficult, UNLESS your spline is a piecewise linear one. Then it is entirely trivial. But few people ever say spline when they mean a piecewise linear (connect-the-dots) function.

xyint = interparc(1000,xy(:,1),xy(:,2),'linear');
plot(xyt(:,1),xyt(:,2),'ob-')
hold on
plot(xy(:,1),xy(:,2),'g+',xyint(:,1),xyint(:,2),'r:')

As you can see, a piecewise linear interpolant will never have a problem, but since it is just connect-the-dots, I'm not sure you would be happy with that result. Any smooth interpolant will potentially have problems with some set of points.

Again, the fact that I've done my examples in 2-d is not relevant, since one can always use projective means to turn a 3-d problem that lies entirely in a plane into a 2-d one. And given my claim that sine a spline interpolant will yield a linear combination of the points, then you don't really need to do the projection at all. Any linear combination of a set of points that lie in a plane will result in a point that also lies in the same plane.

If I've misunderstood the question, please explain and I'll see what I can do.

[John D'Errico]

One day, I'll write a short answer to some question. Really. I will. I promise. I think I will.

[Laurence hutton-smith]

Dear John, firstly I want to say how grateful I am for your answer, this is incredibly useful in improving my understanding of this problem.

I must say however I believe I have not been clear enough, each pair of points in 3space lie on two separate faces of a associated tetrahedral. Imagine a tetrahedral finite element mesh, where we have dotted the centres of many adjacent elements, that is the line we wish to connect. I have attached a picture which describes the problem in 2D (I could not draw it in 3D in an understandable way).

I hope that translates clearly to 3D. It is why I originally describes this as a collection of points, with each pair having an associated volume. I believe it may be possible to amalgamate the tetrahedrals together and simply describe the constraint on the spline to be enclosed by the the union of the tetrahedrals volumes. In which case it is only important the spline can stay inside a well defined volume.

I very much like your example with the triangular face however. The problem you raise with differentiability at the vertices of the constraints might not be a problem here, as the control points are strictly defined to be in the centroids of the tetrahedrals, therefore will ALWAYS lie away from the vertices of the tetrahedral.

In your first image would you be able to explain where the constraint that the points fall within the triangle is implemented? Your interparc function looks very useful, how did you constrain it in this case?

Many thanks again for your assist with this.

[John D'Errico]

Must run out to play bridge soon, so I'll answer in more depth later. That first example I posed had no constraints at all. I used pchip as the interpolant, chosen because pchip is much more robust against creating spurious extrema in the interpolated function.

[John D'Errico]

Ok, so reading your comment, it sounds like you have a series of adjacent tetrahedra, and a set of points on the shared facets? Then you need a smooth curve that passes through those points (interpolates that set), as well as falls entirely inside the volume of this simplicial complex?

So if this is correct, I assume the goal would be to find the spline that is "smoothest" yet still satisfies the desired properties. Here, by smoothest, I'd assume a minimum integral of the square of the second derivative would make sense.

[Laurence hutton-smith]

Yes, you are correct, each tetrahedral has two faces with 1 defined control point on each face, exactly like the drawing, except in 3D.

So yes, we need a curve that passes through all of these points, yet does not meander outside the tetrahedral. In 2D this translates to the curve segment between two triangle faces not escaping outside of their respective triangle.

Yes exactly, however I am unsure of the context of the statement 'minimum integral of the square of the second derivative', as I say, I have never studied splines or numerical methods.

[John D'Errico]

A cubic spline is a mathematical model for the shape of a thin, flexible beam that passes through a set of points. The idea is a thin beam will take the shape of the minimum potential energy stored in that beam, where the energy I refer to is an energy of bending. So a curvy, wiggly beam is not a minimal potential energy state. The beam "wants" to be as straight as possible, and if no constraints were applied, it would be a straight line. We can measure that energy by computing the square of the second derivative of the function, integrating that over the length of the spline.

I think the same idea might apply here, that one would wish to find the smoothest curve that passes through the set of points and satisfies the constraint that the curve lies in a volume. (One might even use the idea of a tension spline, wherein some tension can be placed along the axis of the beam, to help control the shape in that way too.)

The problem with all of this is it is not a simple task. The constraints that an entire curve (composed of infinitely many points) lies inside a polyhedral volume effectively mean that you will need infinitely many linear constraints to ensure your result. Alternatively, it would probably require a set of rather nasty looking nonlinear inequality constraints to deal with.

To solve this problem in any automatic sense really would be very non-trivial. So perhaps a better solution might be the idea of a tension spline, as I mentioned above. Here one formulates a spline with axial tension in the curve. In fact, one can create such a spline where the tension is chosen to vary along the curve. (It has been a while since I played with them, but I recall this being possible.) So the idea would be to use an unconstrained spline fit through the points (i.e., zero tension.) Then in any individual simplex, IF the curve is found to fall outside the bounds of that simplex, then you increase the local tension in that curve until it no longer fails the constraint.

I claim the tension spline idea might be tractable.

[Laurence hutton-smith]

Dear John, thank you so much for the time and thought you've put into answering this question, I believe looking into using tension splines would be a very fruitful course of action for my work. Many thanks, I very much appreciate the guidance in a field I am not familiar with.

Answer 2

[Sean de Wolski]

Here are my thoughts:

If the curve you're trying to fit is linear, lsqlin can handle the inequality constraints you'll need to describe the volume (assuming it's convex). If it's nonlinear, you'll need to frame it for fmincon which can handle those constraints.

If the volume is convex to build the constraints, use vert2con or vert2lcon on the File Exchange. If it is not, I defer to Alan, Matt, or John who will hopefully chime in.

[Laurence hutton-smith]

Dear Sean, thank you very much for your answer and link to this function. So using that function I could constrain the answer through my data points using A_eq x < b_eq, then constrain it's location space with Ax<b?

[Sean de Wolski]

No, not quite. You can't use the equality constraints except for say an intercept.

You would still be minimizing the sum of squared errors but would have to frame the SSE calculation for fmincon.

Is the function non-linear? And is the volume convex? Those are probably good first things to know.