Can I Constrain ListInterpolation by Curvature?

Question

Using image analysis, I have found the positions of a circular ring and imported them as xx and yy coordinates. I am using ListInterpolation to interpolate the data:

xi = ListInterpolation[xx, {0, 1}, InterpolationOrder -> 4, PeriodicInterpolation -> True, Method -> "Spline"];
yi = ListInterpolation[yy, {0, 1}, InterpolationOrder -> 4, PeriodicInterpolation -> True, Method -> "Spline"];

I plot the results as:

splinePlot = ParametricPlot[{xi[s], yi[s]} , {s, 0, 1}, PlotStyle -> {Red}]

and the result looks like:

I am trying to study this shape as it deforms, and I will need to look at the derivatives of this interpretation (notably, second derivatives). I know that there are physical constraints that will not let the local curvature at any point be larger than, for example 1/10 in the units show (so, a radius of curvature 10). Is there a way that I can constrain the interpolation so that the local curvature never exceeds a given value?

Here is the data (Dropbox Link): https://www.dropbox.com/s/g9vajch0obbcplk/testShape.csv?dl=0

Answer 1

An approach to this problem may be found by searching for "person curve" on this site. In particular, I will adapt this answer by @SimonWoods and this answer by @J.M. to How to create a new "person curve"? @MikeY's approach is fundamentally the same, but the details (esp. polar vs. rectangular coordinates) differ. I might be able to add a curvature-based approach, if I have time.

ClearAll[fourierCoefficients, fourierFun, fourierInterpolate, 
  truncate];
(* Fourier coefficients of the coordinates of a list of points *)
fourierCoefficients[x_?MatrixQ] := 
  fourierCoefficients /@ Transpose@x;
fourierCoefficients[x_?VectorQ] := Module[{fc},
   fc = 2 Chop[
      Take[Fourier[x, FourierParameters -> {-1, 1}],
       Ceiling[Length[x]/2]]];
   fc[[1]] /= 2;
   fc
   ];
(* Function from Fourier coefficients
 *   fourierFun[fc][t] is converted to a symbolic-algebraic expression *)
fourierFun[fc_][t_] := fourierFun[fc, t];
fourierFun[fc_, t_] := 
  Abs[#].Cos[Pi (2 Range[0, Length[#] - 1] t - Arg[#]/Pi)] & /@ fc;
(* construct Fourier interpolation *)
fourierInterpolate[x_] := fourierFun[fourierCoefficients@x];
(* truncate a Fourier series equivalent to a least squares projection 
 * onto the low order subspace *)
truncate[fourierFun[fc_], {order_}] := fourierFun[
   With[{norms = Norm /@ Transpose@fc},
    Take[fc, All, Min[1 + order, Length@First@fc]]
    ]
   ];

Example:

xydata = Get["https://pastebin.com/raw/htNS50qt"];

fourierInterpolate[xydata];
{xFN[t_], yFN[t_]} = truncate[fourierInterpolate[xydata], {10}][t]
(*
{113.477 + 32.3787 Cos[π (0.530648 + 2 t)] + 
  0.495141 Cos[π (0.392001 + 4 t)] + 
  0.238509 Cos[π (-0.534649 + 6 t)] + 
  0.258527 Cos[π (-0.288496 + 8 t)] + 
  0.344561 Cos[π (-0.404418 + 10 t)] + 
  0.119943 Cos[π (-0.576067 + 12 t)] + 
  0.0223725 Cos[π (-0.852096 + 14 t)] + 
  0.0350636 Cos[π (0.117051 + 16 t)] + 
  0.0514668 Cos[π (0.516506 + 18 t)] + 
  0.0284801 Cos[π (0.9183 + 20 t)], 
 203.572 + 33.1083 Cos[π (-0.973853 + 2 t)] + 
  0.23916 Cos[π (0.619213 + 4 t)] + 
  0.361602 Cos[π (-0.905991 + 6 t)] + 
  0.135998 Cos[π (0.511184 + 8 t)] + 
  0.341424 Cos[π (0.136349 + 10 t)] + 
  0.0831821 Cos[π (-0.066907 + 12 t)] + 
  0.0387512 Cos[π (0.137154 + 14 t)] + 
  0.0488503 Cos[π (-0.0273644 + 16 t)] + 
  0.0234458 Cos[π (-0.377519 + 18 t)] + 
  0.0300021 Cos[π (0.601693 + 20 t)]}
*)

ParametricPlot[{xFN[t], yFN[t]}, {t, 0, 1}]

Answer 2

You can use a Fourier Transform to take it into the frequency space, then truncate coming back with the inverse transform, thus putting a limit on the max frequency component and therefore max curvature.

Load in the data

data = Import["testShape.csv"];

Center the data about the mean (or median or other measure of centered-ness)

cc = Mean[data];
dc = # - cc & /@ data;
ListPlot[dc]

Convert to polar coordinates, making the first column $\theta$ and second column $r$, and removing a repeated point. The Union also sorts it.

polar = Union@Transpose@Reverse@Transpose@(ToPolarCoordinates /@ dc);

ListPlot@polar

This is a repeating sequence with uneven $\theta$ values. Do an interpolation of it. Add points before and after so the interpolation covers greater than $\{-\pi,\pi\}$

polpolpol = Join[# - {2 π, 0} & /@ polar, polar, # + {2 π, 0} & /@ polar];
if = Interpolation[polpolpol];
Plot[if[θ], {θ, -π, π}];

Now create an evenly subdivided data set that we can do a Fourier transform on. We want $-\pi < \theta \leq \pi$. We can vary the number of points, here I chose 4096.

angles = Rest@Subdivide[-π, π, 2^12];

reg = if /@ angles;

Do a transform to go to the frequency domain space.

ft = Fourier[reg];

Come back to the $\theta$ space, but only carry the first 10 components, dropping the high frequency signal (and therefore high curvature signal).

regTrunc = Re@InverseFourier[PadRight[Take[ft, 10], 2^12]];

Compare

ListPlot[{reg, regTrunc}]

Going back to cartesian...

ListPlot[{FromPolarCoordinates /@ Transpose@{regTrunc, angles}, 
          FromPolarCoordinates /@ Transpose@{reg, angles}}, 
          Joined -> True]

At this point, you have a closed form equation of regTrunc in terms of the Fourier expansion so you can ultimately work through to the curvature using these relationships and bound it. Remainder left as an exercise!

Answer 3

After normalizing the data scaling conveniently we can proceed approximating the data points by a conic (ellipse). This can be done as

f[p_List] := c1 p[[1]]^2 + c2 p[[2]]^2 + c3 p[[1]] p[[2]] + c4 p[[1]] + c5 p[[2]] + c6

data0 = Import["testShape.csv"];
factor = 0.01;
data = data0*factor;
xmax = Max[Transpose[data][[1]]];
xmin = Min[Transpose[data][[1]]];
ymax = Max[Transpose[data][[2]]];
ymin = Min[Transpose[data][[2]]];

vars = {c1, c2, c3, c4, c5, c6};
restr = {c1 > 0.5, c2 > 0.5, c1 c2 > 0.25 c3^2};
obj = Sum[f[data[[k]]]^2, {k, 1, Length[data]}];
sol = NMinimize[{obj, restr}, vars]
conic = f[{x, y}] /. sol[[2]]
gr1 = ContourPlot[conic == 0, {x, xmin, xmax}, {y, ymin-0.1, ymax}, ContourStyle -> Black];
gr2 = ListPlot[data, PlotStyle -> Red];
Show[gr2, gr1, AspectRatio -> 1]

NOTE

The scaling can be done with the values

xmax = Max[Transpose[data][[1]]]
xmin = Min[Transpose[data][[1]]]
ymax = Max[Transpose[data][[2]]]
ymin = Min[Transpose[data][[2]]]

In the present case we choose a multiplicative factor of $0.01$ over data

Answer 4

Again I tend to duplicate, despite this is an interpolation and a real 2D curve under consideration.

What I did is using this solution. I like to meet the periodic boundary conditions, but did not really match them.

Use this

xxn = Mean[{#, RotateLeft@#}] &@(Transpose[xydata][[1, All]]);
yyn = Mean[{#, RotateLeft@#}] &@(Transpose[xydata][[2, All]]);

It do in this particular case not really enhance the problems.

Better make use of

xxn7 = MovingAverage[ArrayPad[#, {0, 1}, "Periodic"], 7] &@xxn;
xi1 = ListInterpolation[xxn7, {0, 1}, InterpolationOrder -> 4, 
  Method -> "Spline"]
Plot[xi1'[x], {x, 0, 1}]

The result is:

That compares to

.

The curvature is defined:

This is not so critical in this particular case because it is inversely dependent on the first derivative at the first look. The second derivative called a in the formular is even oscillating wilder at the problematic points. So smoothing the first derivative and then derive further is a nice selection for controlling and stay close to the real process of oscillations.

There is no general best approximation. For most of the points the Mean smoothing, neating, even, sleeking is to much, for others fair, for some not enough. The last ones govern the seconde derivative and make it even more jumpy.

On accident or selection, the second derivative from Your interpolation is already limited. 60000 is pretty much. MovingAverage of 7 neighbors reduces the maximum to 20000. There is a need for a function between neighbors for MovingAverage and the Maximum in x''.

Better is average different pairs of neighbors repetitive. The resulting curve and maximum value of the second derivative will get smaller and smaller. Stay on the original curve at the points fitting best, showing a little value in the second derivative to keep Your moving curve close to the original process. It is a matter of taste or mathematical authentic perception.

The path of calculation a parametric representation of the data points may be shorter but giving less control and insight into the behavior of the first and second derivative of the original curve.

If You decide for the MovingAverage approach, it might be sense full to transform the center of the curve to the origin of the coordinate system.

So after applying the MovingAverage over 7 neighbors, the second derivative is tamed:

The spikes a smaller about a third but there are more spikes. The distribution broadened a lot. This again contributes to the advice that a selection for better points might be better.

Keep in mind I did not alter the coordinate system nor did I add or remove points to or from the original curve. I worked on approximate curvature.

The resulting curve is a wavy as the given one:

The problem is it needs to be closed again.

A close comparison shows the power of averaging little:

In the upper half the averaging cuts of the spikes in the lower half it tends to move the curve outwards, enlarging the radius. So no symmetry at all. It seems there is more work ahead in dealing with the completion of the curve at the bottom. Again as with all steps it is up to the realism that applies.

Again I remembered not to make too much duplicate. For my introduce problem about the periodic option have a look at this article:

Higher order periodic interpolation (curve fitting). The answer by J. M. is in limbo should do the rest.

Obey the periodic option and stay away from evaluating the curvature but look deep and with care a the first and second derivative of the numerical interpolated representation of the given data.

I am just having problems with the part from 20 minutes past to close before half-past. The rest is comfortable and fine. There should be a rest of the curls in the curve. And the curve must not be as induced by changing to polar coordinate reside inside the very random speedy given one. There is a systematic error from the acquisition at the bottom that is in need to be dealt different with than numeric work.