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

curve fittingMATLAB

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.

Best Answer

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.
Related Question