2
$\begingroup$

I generated a spline function called f using the BSplineFunction with domain {t, 0, 1}.

f = BSplineFunction[{{0., 5.81152}, {-0.909122, 5.73997}, {-1.86805, 5.5031}, 
  {-1.79586, 5.52709}, {-2.40811, 5.28912}, {-3.4109, 4.70527}, {-4.68131,   3.501}, 
  {-4.93131, 2.99958}, {-5.09697, 2.6673}, {-5.34697, 2.16588}, {-5.59697, 1.66446}, 
  {-5.84697, 1.16304}}, SplineWeights -> {10, 10, 10, 1, 1, 1, 1, 1, 1, 10, 10, 10}]

I then wanted to find the points at which the function intersects the piecewise linear function joining the spline's control points. Graphically, I know the intersection occurs within the domain of f.

k[t_] = {-4.93131 - 0.915666 t, 2.99958 - 1.83654 t}

However, running the code:

FindRoot[First[k[0]] == First[f[t]], {t, .5, 0, 1}]

where k[0] is the end of the line, produces a negative number. Is there a way to solve for the intersection points of those two functions?

$\endgroup$
4
  • $\begingroup$ I Believe that you should provide the functions. $\endgroup$
    – Öskå
    Jul 1, 2014 at 13:48
  • $\begingroup$ There are the functions! @Öskå $\endgroup$
    – Kaisey
    Jul 1, 2014 at 16:15
  • 1
    $\begingroup$ What's f? Did you mean g? Plus there are more problems (in addition to the extra square bracket in your FindRoot). For instance, First[g[t]] just yields t. So your FindRoot statements reduces to FindRoot[-4.93131==t,{t,0.5}] which will return t->-4.9313. $\endgroup$ Jul 1, 2014 at 16:34
  • $\begingroup$ @rhomboidRhipper The syntax problems should be fixed. f[t] does not return a value of t...it yields x and y coordinates for the domain 0<t<1. FindRoot[First[k[0]] == First[f[t]], {t, .5, 0, 1}] should yield the value of t at which the x coordinate is equal to -4.93131. Is that correct or am I completely missing something? $\endgroup$
    – Kaisey
    Jul 2, 2014 at 13:34

4 Answers 4

4
$\begingroup$

This is yet another problem that can be solved by explicitly representing the B-spline curve with BSplineBasis[], as was done here:

pts = {{1, 1}, {2, 3}, {3, -1}, {4, 1}, {5, 0}};

(* piecewise linear interpolant *)
lin = Interpolation[pts, InterpolationOrder -> 1];

n = Min[Length[pts] - 1, 3]; (* B-spline degree *)
m = Length[pts];
(* clamped uniform knots for B-spline *)
knots = {ConstantArray[0, n + 1], Range[m - n - 1]/(m - n),
         ConstantArray[1, n + 1]} // Flatten;

{xu, yu} = Transpose[pts];

(* B-spline component functions *)
bf[t_] = xu.Table[BSplineBasis[{n, knots}, i - 1, t], {i, Length[pts]}];
bg[t_] = yu.Table[BSplineBasis[{n, knots}, i - 1, t], {i, Length[pts]}];

From here we use a trick similar to what was done in this answer; use the MeshFunctions option of Plot[] to find initial estimates for the intersection, and then polish those estimates with FindRoot[]:

crs = Cases[Normal[Plot[bg[t] - lin[bf[t]], {t, 0, 1},
                        Mesh -> {{0}}, MeshFunctions -> {#2 &}]], 
            Point[{x_, _}] :> (Through[{bf, bg}[\[FormalT]]] /. 
                               FindRoot[bg[\[FormalT]] - lin[bf[\[FormalT]]],
                                        {\[FormalT], x}]), ∞]
   {{2.3914022076831603, 1.4343911692673563}, {3.707060176019142, 0.4141203520382841}}

Visualize the geometry:

ParametricPlot[{bf[t], bg[t]}, {t, 0, 1}, 
               Epilog -> {{Directive[AbsoluteThickness[1.6], ColorData[97, 2]], 
                           Line[pts]},
                          {Directive[AbsolutePointSize[6], ColorData[97, 3]], 
                           Point[pts]},
                          {Directive[AbsolutePointSize[8], ColorData[97, 4]], 
                           Point[crs]}},
               PlotRange -> {Automatic, {-1, 3}}, PlotRangePadding -> Scaled[.05]]

B-spline and piecewise linear function, with intersections

$\endgroup$
1
$\begingroup$
Clear["Global`*"]
pts = {{1, 1}, {2, 3}, {3, -1}, {4, 1}, {5, 0}};
f1 = Interpolation[pts, InterpolationOrder -> 1];
f2 = BSplineFunction[pts];
Show[Graphics[{Red, Point[pts], PointSize[0.02], 
   Point[f2 /@ {0.2, 0.4}], Green, Line[pts]}, Axes -> True], 
   ParametricPlot[f2[t], {t, 0, 1}]]

figure

First,determine approximate range of t, in my example tmin = 0.2, tmax = 0.4.

Then,use method of dichotomy to find the point we need.

dichotomy[{t1_, t2_}] :=
 Module[{x1, x2, y1, y2, tmid, x3, y3},
  {x1, x2} = First@f2[#] & /@ {t1, t2};
  {y1, y2} = f1 /@ {x1, x2};
  tmid = (t1 + t2)/2;
  {x3, y3} = f2[tmid];
  If[(y1 - Last@f2@t1)*(f1[x3] - Last@f2@tmid) < 0, {t1, tmid}, {tmid, t2}]
  ]
tmin = 0.2; tmax = 0.4; error = 10^-5.;
ans = First@NestWhile[dichotomy, {0.2, 0.4}, 
     With[{p = f2@First@#}, f1@First@p - Last@p > error] &]

0.306755

Last, check the answer:

f2[ans]
f1@First@f2[ans]

{2.3914, 1.4344}

1.4344

$\endgroup$
2
  • $\begingroup$ Thank you! @Chenminqi Although, with the shorthand, I'm having trouble understanding the dichotomy function. Would you mind explaining? $\endgroup$
    – Kaisey
    Jul 1, 2014 at 16:46
  • $\begingroup$ @Kaisey Does this answer reply your concerns? $\endgroup$
    – Öskå
    Jul 1, 2014 at 18:07
1
$\begingroup$

Define the spline for numerical values only and solve for x and y coordinate simultaneously.

g[t_?NumericQ] = f[t]

fr = FindRoot[{First[k[0]], a} == g[t], {t, .7}, {a, 1}]

(*   {t -> 0.674177, a -> 2.99958}   *)

{f[t], k[0]} /. fr

(*   {{-4.93131, 2.99958}, {-4.93131, 2.99958}}   *)
$\endgroup$
1
$\begingroup$

Since V10, one can use Indexed[.., 1] instead of First[..] or Part:

FindRoot[Indexed[k[0], 1] == Indexed[f[t], 1], {t, .5, 0, 1}]
(*  {t -> 0.674177}  *)
$\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.