3
$\begingroup$

I have coordinates of points that represent a curve (for example, helix or centerline of moebius strip) in $ \mathbb{R^3}$. I want to fit a Bspline through these points. How to do this in Mathematica?

Example of curve: see at the end

I get noisy values for curvature and torsion of this curve if I use the pts as control points of BSplineCurve function.

Also, I have used Interpolation function, with that too, computed curvature and torsion are noisy. The real data I want to fit:

{{-0.127862, 7.23797, 0.080385}, {-0.386039, 7.22523, 
  0.238487}, {-0.643736, 7.19972, 0.395825}, {-0.900618, 7.16136, 
  0.551879}, {-1.15633, 7.11002, 0.706108}, {-1.4105, 7.04555, 
  0.857945}, {-1.66271, 6.96778, 1.00679}, {-1.9125, 6.8765, 
  1.15199}, {-2.15936, 6.77151, 1.29287}, {-2.40271, 6.6526, 
  1.42869}, {-2.64192, 6.51958, 1.55865}, {-2.87627, 6.37229, 
  1.68191}, {-3.10494, 6.21063, 1.79756}, {-3.32704, 6.03455, 
  1.90464}, {-3.54155, 5.84411, 2.00217}, {-3.74738, 5.63948, 
  2.08911}, {-3.9433, 5.42099, 2.1644}, {-4.12798, 5.18913, 
  2.227}, {-4.30001, 4.94459, 2.27587}, {-4.45787, 4.68829, 
  2.31007}, {-4.59997, 4.42139, 2.3287}, {-4.72468, 4.14532, 
  2.33102}, {-4.83036, 3.8618, 2.31646}, {-4.91539, 3.57279, 
  2.28468}, {-4.9782, 3.28055, 2.23561}, {-5.01738, 2.98756, 
  2.16947}, {-5.03167, 2.69648, 2.08687}, {-5.02007, 2.41015, 
  1.98876}, {-4.98187, 2.13143, 1.87651}, {-4.9167, 1.86318, 
  1.75186}, {-4.82458, 1.60813, 1.61691}, {-4.70594, 1.36881, 
  1.47412}, {-4.56163, 1.14739, 1.32615}, {-4.39291, 0.945657, 
  1.17587}, {-4.2014, 0.764887, 1.02622}, {-3.98905, 0.605816, 
  0.880086}, {-3.75804, 0.4686, 0.740214}, {-3.51073, 0.352823, 
  0.609092}, {-3.24952, 0.257522, 0.488851}, {-2.97682, 0.181251, 
  0.381183}, {-2.6949, 0.122163, 0.287282}, {-2.40588, 0.0781093, 
  0.207801}, {-2.11161, 0.0467557, 0.142845}, {-1.81367, 0.0256998, 
  0.0919741}, {-1.51331, 0.0125879, 0.0542331}, {-1.21151, 0.00522497,
   0.0282002}, {-0.908944, 0.00167311, 0.0120463}, {-0.60604, 
  0.00033535, 0.00360597}, {-0.303029, 0.0000222426, 0.000456223}, {0,
   0, 0}, {0, 0, 0}, {0.303029, 0.000022046, -0.000449164}, {0.606041,
   0.000333794, -0.00357979}, {0.908945, 
  0.00166793, -0.0119886}, {1.21151, 0.00521293, -0.028099}, {1.51332,
   0.012565, -0.0540773}, {1.81368, 0.0256614, -0.0917535}, {2.11165, 
  0.0466967, -0.142551}, {2.40594, 0.0780248, -0.207425}, {2.69499, 
  0.122048, -0.286817}, {2.97694, 0.181102, -0.380626}, {3.24969, 
  0.257335, -0.488198}, {3.51095, 0.352596, -0.608342}, {3.75834, 
  0.468332, -0.739366}, {3.98943, 0.605506, -0.879142}, {4.20188, 
  0.764539, -1.02518}, {4.39349, 0.945273, -1.17475}, {4.56232, 
  1.14698, -1.32494}, {4.70675, 1.36837, -1.47283}, {4.82551, 
  1.60768, -1.61555}, {4.91775, 1.86271, -1.75043}, {4.98305, 
  2.13096, -1.87502}, {5.02137, 2.40969, -1.98722}, {5.0331, 
  2.69604, -2.08527}, {5.01892, 2.98713, -2.16782}, {4.97985, 
  3.28015, -2.23391}, {4.91715, 3.57242, -2.28293}, {4.83222, 
  3.86146, -2.31467}, {4.72663, 4.14502, -2.32918}, {4.602, 
  4.42112, -2.32682}, {4.45998, 4.68806, -2.30815}, {4.30219, 
  4.9444, -2.27392}, {4.13022, 5.18897, -2.225}, {3.94559, 
  5.42087, -2.16236}, {3.74972, 5.63939, -2.08703}, {3.54394, 
  5.84404, -2.00006}, {3.32946, 6.03451, -1.90249}, {3.1074, 
  6.21061, -1.79537}, {2.87875, 6.37229, -1.67969}, {2.64443, 
  6.51959, -1.5564}, {2.40524, 6.65262, -1.42641}, {2.1619, 
  6.77154, -1.29057}, {1.91506, 6.87653, -1.14966}, {1.66529, 
  6.96781, -1.00444}, {1.41309, 7.04558, -0.855578}, {1.15893, 
  7.11005, -0.703726}, {0.903221, 7.16138, -0.549484}, {0.646343, 
  7.19974, -0.393421}, {0.38865, 7.22524, -0.236076}, {0.130474, 
  7.23797, -0.077971}, {-0.127862, 7.23797, 0.080385}}
Improve this question
$\endgroup$
3
  • $\begingroup$ Graphics3D[ BSplineCurve[pts] ] only gives a Graphics3D object. How did you calculate curvature and torsion? $\endgroup$
    – Ulrich Neumann
    Commented Oct 30, 2019 at 11:39
  • $\begingroup$ These are standard formulae for curvature and torsion compuation (ref: A Persseley- "Elementary Differential Geometry"), where t is any parameter. In the code below Ifn is the interpolated function. Curvature: curvK = (Dot[Cross[Ifn''[t], Ifn'[t]] , Cross[Ifn'[t], Ifn'[t]]]/ (Dot[Ifn'[t], Ifn'[t]])^3 )^(1/2) Torsion of curve: torsionT = Dot[ Cross[Ifn'[t] , Ifn''[t]], Ifn'''[t]] / Dot[ Cross[Ifn'[t], Ifn''[t]] , Cross[Ifn'[t], Ifn'[t]] ] $\endgroup$
    – akr
    Commented Oct 30, 2019 at 11:51
  • $\begingroup$ What is the definition of Ifn[t]? $\endgroup$
    – Ulrich Neumann
    Commented Oct 30, 2019 at 11:55

3 Answers 3

Reset to default
6
$\begingroup$

BezierCurve is only a Graphics3D-object. Use BezierFunction to describe a function.

Examplary here I give the calculation of the curvature kappa of the BezierFunction

pts = Table[{2 Cos[t], 2 Sin[t], 3*t}, {t, 0, 6, 0.01}];
x = BezierFunction[pts]

kappa = Sqrt[#.#] &[ D[x'[u]/Sqrt[x'[u].x'[u]], u]/Sqrt[x'[u].x'[u]]];
Plot[kappa, {u, 0, 1}, PlotRange -> {0, Automatic}]

enter image description here

which shows a smooth curvature function.

answer to the modified question

Let's call the data you gave p, then removing duplicate points

p=p// DeleteDuplicates;

Calculate arclength 

si = Accumulate[Join[{{0}}, Map[Sqrt[#.#] &, Rest[pi] - Most[pi]]] ] // Flatten

Interpolate the points with piecewise splines

ip = Interpolation[MapThread[{#1, #2} &, {si, pi}], Method -> "Spline"]
ParametricPlot3D[ip[s], {s, 0, Max[si]}]

enter image description here

The curvature can be calculated as in the first part of my answer.

Improve this answer
$\endgroup$
8
  • $\begingroup$ Thank you for your reply. But it is still taking pts as control points, not fitting a curve through pts. Though the curvature and torsion are smooth now, the curve I get after selecting with BezierCurve far away from the pts. $\endgroup$
    – akr
    Commented Oct 30, 2019 at 12:50
  • $\begingroup$ I have edited the question and have included the real data points to fit. $\endgroup$
    – akr
    Commented Oct 30, 2019 at 12:57
  • $\begingroup$ Ok , what is the curve parameter of your new curve? $\endgroup$
    – Ulrich Neumann
    Commented Oct 30, 2019 at 13:09
  • $\begingroup$ It is arc legth which I intend to compute later. But you may use any equally spaced parameter, including the one in range [0,1]. $\endgroup$
    – akr
    Commented Oct 30, 2019 at 15:19
  • 1
    $\begingroup$ This demonstation approximates the points, you asked for interpolation ! $\endgroup$
    – Ulrich Neumann
    Commented Oct 30, 2019 at 16:21
4
$\begingroup$

To eliminate noise in the data, it is best to perform a smoothing and not an interpolation. With the procedure that follows, three splines are adjusted one for each dimension, being that the degree of adjustment and smoothing is associated with the number of intermediate nodes. We will choose a grade 3 spline to make the adjustments. Given the data set raw we proceed as follows:

si = Accumulate[Join[{{0}}, Map[Sqrt[#.#] &, Rest[raw] - Most[raw]]]] // Flatten
X = Table[{si[[k]], raw[[k, 1]]}, {k, 1, Length[si]}];
Y = Table[{si[[k]], raw[[k, 2]]}, {k, 1, Length[si]}];
Z = Table[{si[[k]], raw[[k, 3]]}, {k, 1, Length[si]}];

SplineModel[data_, deg_, knots_] := Block[{basis, allKnots},
  basis = Array[\[FormalX]^# &, deg + 1, 0]~Join~Table[BSplineBasis[{deg, knots}, i, \[FormalX]], {i, 0, Length[knots] - deg - 2}];
  LinearModelFit[data, basis, \[FormalX]]];

knots = Range[0, 30, 3];
modX = SplineModel[X, 3, knots];
modY = SplineModel[Y, 3, knots];
modZ = SplineModel[Z, 3, knots];
ParametricPlot3D[{modX[s], modY[s], modZ[s]}, {s, 0, 30}]

So in knots = Range[0, 30, 3]; we choose the number of knots. The present values are choose as an adequate minimum.

Improve this answer
$\endgroup$
1
$\begingroup$ 
searchSpan[knots_, u0_] :=
 With[{max = Max[knots]},
  If[u0 == max,
   Position[knots, max][[1, 1]] - 2,
   Ordering[UnitStep[u0 - knots], 1][[1]] - 2]
  ]

calcFitKnots[{m_, n_}, deg_, paras_] :=

 With[{d = (m + 1)/(n - deg + 1)},
  Join[ConstantArray[0, deg + 1],
   (*calculate the interior knots*)
   Function[{j},
     Module[{i, alpha}, i = Floor[j*d];
      alpha = j*d - i;
      (1 - alpha) paras[[i]] + alpha paras[[i + 1]]]
     ] /@ Range[1, n - deg],
   ConstantArray[1, deg + 1]]
  ]

calcParas[pts_, type_] :=
 Which[
  type === Automatic || type === "ChordLength",
  FoldList[
    Plus, 0, Normalize[(Norm /@ Differences[pts]), Total]] // N, 
  type === "Centripetal",
  FoldList[
    Plus, 0, Normalize[(Norm /@ Differences[pts])^(1/2), Total]] // N,
   type === "EqualSpaced",
  Range[0, 1, 1/(Length@pts - 1)] // N
  ]

(* Spline Fit *)
(* Data from simulations *)

pts = data3;
m = Length@pts - 1;
ControlPointsNumber = IntegerPart[Length[pts]/2];
splineDegree = 5;
Parametrization = "ChordLength";
(*achieve the value of options*)
cpn = ControlPointsNumber
pz = Parametrization
sd = splineDegree



paras = calcParas[pts, pz];
(*calculate the knots*)

knots = calcFitKnots[{m, cpn - 1}, sd, paras];
(*calculate the coefficients of matrix*)
coeffMat =
  (Function[{u0},
      With[{i = searchSpan[knots, u0]},
       Join[ConstantArray[0, i - sd],
        BSplineBasis[{sd, knots}, #, u0] & /@ Range[i - sd, i],
        ConstantArray[0, cpn - 1 - i]]]] /@ ArrayPad[paras, -1])[[All,
    2 ;; -2]];

(*solve the control points of the B-Spline curve*)

R = Transpose[coeffMat].(pts[[2 ;; -2]] -
     With[{Q0 = First@pts, Qm = Last@pts},
      BSplineBasis[{sd, knots}, 0, #] Q0 +
         BSplineBasis[{sd, knots}, cpn - 1, #] Qm & /@ 
       ArrayPad[paras, -1]]);
ctrlpts =
  Join[{First@pts},
   LinearSolve[[email protected], R], {Last@pts}];
(*visualize the fitting result*)
p1 = Graphics3D[
   BSplineCurve[ctrlpts, SplineDegree -> sd, SplineKnots -> knots]
   ];
p2 = ListPointPlot3D[pts];
Show[p1, p2]


Ifn3 = BSplineFunction[ctrlpts, SplineDegree -> sd, 
  SplineKnots -> knots, SplineClosed -> True]
curvK3I = (Dot[Cross[Ifn3''[t], Ifn3'[t]] , 
      Cross[Ifn3''[t], Ifn3'[t]]]/ (Dot[Ifn3'[t], Ifn3'[t]])^3 )^(1/2);
torsion3I = 
  Dot[ Cross[Ifn3'[t] , Ifn3''[t]], Ifn3'''[t]] / 
   Dot[ Cross[Ifn3'[t], Ifn3''[t]] , Cross[Ifn3'[t], Ifn3''[t]] ];
plotK3I = 
 Plot[Evaluate[{curvK3, curvK3I}], {t, 0., 1}, 
  PlotRange -> {{0., 01}, {-0.1, 0.6}}]
plotT3I = Plot[torsion3I, {t, 0., 1}]
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.