0
$\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$
  • $\begingroup$ I Believe that you should provide the functions. $\endgroup$ – Öskå Jul 1 '14 at 13:48
  • $\begingroup$ There are the functions! @Öskå $\endgroup$ – Kaisey Jul 1 '14 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$ – rhomboidRhipper Jul 1 '14 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 '14 at 13:34
2
$\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}]]

enter image description here

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 ans

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

{2.3914, 1.4344}

1.4344

$\endgroup$
  • $\begingroup$ Thank you! @Chenminqi Although, with the shorthand, I'm having trouble understanding the dichotomy function. Would you mind explaining? $\endgroup$ – Kaisey Jul 1 '14 at 16:46
  • $\begingroup$ @Kaisey Does this answer reply your concerns? $\endgroup$ – Öskå Jul 1 '14 at 18: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.