3
$\begingroup$

This my question is concerned with previous about Tube cross-section. Let's use the same pts as in that question.

pts = {{0.`, 0.`, 0.`}, {7.497665500157259`*^-14, 
   0.00036555534709634656`, 
   0.041887873232098444`}, {2.9988378071419946`*^-13, 
   0.0014621100333863195`, 
   0.08376298663421493`}, {6.746528662189823`*^-13, 
   0.0032893300277940163`, 
   0.12561258426325045`}, {1.1991697498398906`*^-12, 
   0.005846658724918654`, 
   0.16742391794868797`}, {1.8732746539195436`*^-12, 
   0.009133317114586991`, 
   0.20918425117592362`}, {2.696762233490476`*^-12, 
   0.013148304019154402`, 
   0.25088086296604745`}, {3.669381638827055`*^-12, 
   0.017890396398482316`, 
   0.2925010517508915`}, {4.790836591859004`*^-12, 
   0.02335814972249914`, 
   0.33403213924216535`}, {6.0607854764232425`*^-12, 
   0.029549898411231134`, 
   0.37546147429349974`}, {7.478841442326563`*^-12, 
   0.03646375634216925`, 
   0.41677643675422266`}, {9.044572523187444`*^-12, 
   0.04409761742481735`, 
   0.45796444131369324`}, {1.0757501768021084`*^-11, 
   0.052449156242246726`, 
   0.499012941335023`}, {1.26171073865276`*^-11, 0.06151582875946169`,
    0.5399094326770161`}, {1.4622822908039118`*^-11, 
   0.07129487309836018`, 
   0.5806414575031648`}, {1.6774037354077335`*^-11, 
   0.08178331037905351`, 0.6211966080765395`}, {1.9070095424469`*^-11,
    0.09297794562728895`, 
   0.6615625305394178`}, {2.1510297696962602`*^-11, 
   0.10487536874769866`, 
   0.7017269286765003`}, {2.4093900840285483`*^-11, 
   0.11747195556257856`, 
   0.7416775676605682`}, {2.6820117840576448`*^-11, 
   0.13076386891588063`, 
   0.7814022777794407`}, {2.9688118241124915`*^-11, 
   0.14474705984208247`, 0.8208889581430966`}}

Now I need to construct parametric form of the curve. One method is using Interpolation:

path[u_] = Interpolation[#, u] & /@ Transpose[pts]

The main advantage of Interpolation is that path is going through all the points of pts exactly. Then I use code from @SquareOne:

frenet[u_] = FrenetSerretSystem[path[u], u][[2]]; // AbsoluteTiming
transform[u_] := 
  Composition[TranslationTransform[pts[[u]]], 
   FindGeometricTransform[
     frenet[u], {{0, 0, 1}, {1, 0, 0}, {0, 1, 0}}][[2]]];

I included AbsoluteTiming here to measure time of operation. And I see that frenet[u_] takes significant time, about 0.27--0.3 sec on my machine.

Another approach is using BSplineFunction:

path1[u_] := BSplineFunction[#][u] & /@ (Transpose@pts);
frenet1[u_] = FrenetSerretSystem[path1[u], u][[2]]; // AbsoluteTiming

This is much faster, it takes about 0.09--0.1 sec on my machine, i.e. BSplineFunction is 2.5--3 times faster than Interpolation. Why is this difference in speed? Can I speed-up Interpolation if I want to use it because of (as I think, probably wrong) it gives more accurate approximation.

Thank you.

$\endgroup$
  • $\begingroup$ What happens if you use Method -> "Spline" in Interpolation[]? $\endgroup$ – J. M. will be back soon Dec 19 '16 at 10:12
  • $\begingroup$ If I run path[u_] = Interpolation[#, u, Method -> "Spline"] & /@ Transpose[pts], then frenet[u_] takes about 0.32--0.35 sec, more or the same time or even longer. $\endgroup$ – Alx Dec 19 '16 at 10:28
  • $\begingroup$ Just to be clear, you want to speed up the construction of the code generated by FrenetSerretSystem, which is done only once for each set of pts? Interpolation itself is fast, but it is FrenetSerretSystem that takes the extra time when analyzing (I surmise) InterpolatingFunction compared to analyzing BSplineFunction. $\endgroup$ – Michael E2 Dec 19 '16 at 18:33
  • $\begingroup$ @MichaelE2, right. I have many such sets of points, and I had a hope it is possibe to decrease time of these operations. $\endgroup$ – Alx Dec 19 '16 at 23:21
0
$\begingroup$

Perhaps this will help: FrenetSerretSystem seems to take a symbolic approach, which slows it down. A more straightforward numeric approach, which will be much faster, may be more satisfactory.

The following function works for all (twice differentiable) paths. The savings are two-fold. One doesn't have to compute a new Frenet frame function for each set of points; only an interpolating path needs to be computed. Second, as seen below, this function executes faster than the one created with FrenetSerretFrame[] in the OP.

ClearAll[frenet2];
fsfC = Compile[{{v, _Real, 1}, {a, _Real, 1}},
   Block[{nv, na},
    If[v == {0, 0, 0},                        (* singular point *)
     {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
     nv = v * Internal`ReciprocalSqrt[v.v];   (* use v / Sqrt[v.v] pre V10 *)
     If[a == {0, 0, 0},                       (* no turning *)
      {nv, {0, 0, 0}, {0, 0, 0}},
      na = # * Internal`ReciprocalSqrt[#.#] &[a - (a.nv) nv];
      {nv, na, Cross[nv, na]}
      ]]
    ]];
frenet2[path_, {u_, u0_?NumericQ}] :=
  With[{v = D[path, u]}, With[{a = D[v, u]}, Block[{u = u0},
     fsfC[v, a]
     ]]];

Using a single interpolating function instead of three seems to speed things up a little, too:

path2[u_] = Interpolation[MapIndexed[{#2, #1} &, pts]][u]

Timings. Note that frenet2 works on the OP's paths. The frame is slightly different for path1 because the B-spline path is slightly different from the interpolated path.

uu0 = 12.5;
frenet2[path[u], {u, uu0}] // RepeatedTiming
frenet2[path1[u], {u, Rescale[uu0, {1, 21}]}] // RepeatedTiming
frenet2[path2[u], {u, uu0}] // RepeatedTiming
(*
  {0.00013,
    {{4.08916*10^-11, 0.199371, 0.979924},
     {2.00986*10^-10, 0.979924, -0.199371},
     {-1., 2.05103*10^-10, 1.9387*10^-26}}}

  {0.00016,
    {{4.03653*10^-11, 0.196805, 0.980443},
     {2.01092*10^-10, 0.980443, -0.196805},
     {-1., 2.05103*10^-10, 9.69352*10^-26}}}

  {0.000077,
    {{4.08916*10^-11, 0.199371, 0.979924},
     {2.00986*10^-10, 0.979924, -0.199371},
     {-1., 2.05103*10^-10, 1.9387*10^-26}}}
*)

The first and third agree with the frame computed by FrenetSerretFrame used in the OP's frenet; they're also much faster than frenet:

frenet[uu0] // RepeatedTiming
(*
  {0.0011,
    {{4.08916*10^-11, 0.199371, 0.979924},
     {2.00986*10^-10, 0.979924, -0.199371},
     {-1., 2.05103*10^-10, 1.28793*10^-26}}
*)
$\endgroup$
  • $\begingroup$ many thanks! This works very fast! $\endgroup$ – Alx Dec 20 '16 at 2:57
  • $\begingroup$ sorry for late question, had no time. As I undestand, in your fsfC parameters nv and na are normalized. So we can write nv = Normalize[v], and na = Normalize[a - (a.nv) nv]. This works. Next step is replace (a.nv) nv with Projection[a, nv], and this gives me errors that parameters have different rank and uncompiled function will be used. What is wrong? I read many times that one shoud preferably use built-in functions, they are extremely optimized. $\endgroup$ – Alx Dec 22 '16 at 9:17
  • $\begingroup$ @Alx: Yes, built-in are often (should be always but in practice not) to be preferred. If we're using Compile, then we should use functions built into the compiler. Execute the following: Needs["CompiledFunctionTools`"]; then CompilePrint@Compile[{{v, _Real, 1}, {w, _Real, 1}}, Projection[v, w]] and CompilePrint@Compile[{{v, _Real, 1}}, Normalize[v]]. They show MainEvaluate[..] called on Projection and Normalize, which means the functions are not "compilable"..... $\endgroup$ – Michael E2 Dec 22 '16 at 14:28
  • $\begingroup$ @Alx ...Further Compile has no (internal) declared return type for Normalize[], so it assumes it returns a real number. You can declare the return type by added a declaration at the end: fsfC = Compile[{{v, _Real, 1}, {a, _Real, 1}}, Block[{nv, na},,,], {{Normalize[_], _Real, 1}}]. See "Compiled Function Operation" in the compiler manual. Just to be clear, the use of Projection and Normalize inside Compile should be slower. $\endgroup$ – Michael E2 Dec 22 '16 at 14:37
  • $\begingroup$ thank you for clarification. I've read topic on List of compilable functions and the main conclusion (as I understood) is that (in most cases) simple ("low-level") functions such as Sin or Dot can be used in Compile, but not "high-level" already optimized ones as e.g. NIntegrate or Projection. $\endgroup$ – Alx Dec 22 '16 at 15:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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