Find optimal BezierFunction

For a given set of points (examplary quadrant)

n = 5; (* number of inner points*)
p = Table[{Cos[φ], Sin[φ]}, {φ, Join[{0}, RandomReal[{0, Pi/2}, n], {Pi/2}]}];


I want to calculate the optimal controlpoint {kx,ky} of a 3point bezier function(using BezierFunction).

An optimal controlpoint seems to exist for this problem:

Manipulate[
Show[{
ParametricPlot[
BezierFunction[{p[[1]], {kx, ky }, p[[-1]]}][u] // Evaluate, {u, 0, 1}],
ListPlot[p],
Graphics[{Red, Point[{kx, ky}],Text["{kx,ky}", {kx, ky}, {0, -1}]}]}],
{{kx, .7}, 0, 1, Appearance -> "Labeled"}, {{ky, .7}, 0, 1, Appearance -> "Labeled"}]


I tried

some variable bezier parameters(to be optimized too):

ui = Table[u[i], {i, 2, Length[p] - 1}];


general bezier curve with variable control point :

bez[k_ ] := BezierFunction[k]


But optimization gives several errors

NMinimize[Sum[Norm[ bez[{p[[1]], {kx, ky}, p[[-1]]}][u[i]] - p[[i]]],
{i, 2, Length[p] - 1}],
Join[{kx, ky}, ui] ]
(*BezierFunction::invcpts: {{1,0},{kx,ky},{0,1}} should be a rectangular array of
machine-sized real numbers of any depth,
whose dimensions are greater than 1., ...*)


The argument of bez[...] is rectangular, though I don't understand the message.

What's wrong with my attempt? Thanks!

=> SOLUTION: Thanks a lot @SjoerdSmid

With userdefined bezier function bezCurve

bezCurve[{pt : {_, _}}] := pt &;
bezCurve[ctrlPts_?MatrixQ] := bezCurve[ctrlPts] = With[{
b1 = bezCurve[Most[ctrlPts]],b2 = bezCurve[Rest[ctrlPts]]}
, (1 - #)*b1[#] + #*b2[#] &];
bezCurve[ctrlPts_?MatrixQ, t_] := Simplify[bezCurve[ctrlPts][t]];


the optimization is done in ~7seconds

opt = NMinimize[{Sum[Norm[bezCurve[{p[[1]], {kx, ky}, p[[-1]]}][u[i]] - p[[i]]],
{i, 2,Length[p] - 1}],
Join[{0}, ui, {1}] /. List -> Less, 0 <= kx <= 1 , 0 <= ky <= 1},
Join[{kx, ky}, ui ]]
(*{1.22189, {kx -> 0.899875, ky -> 0.999976,
u[2] -> 0.152739, u[3] -> 0.152739,
u[4] -> 0.152739, u[5] -> 0.152739,u[6] -> 0.244299}}*)

• Definition of bez is wrong. Just test it with any valid argument. – xzczd Oct 31 '19 at 8:44
• @xzczd Thanks, I tried Show[{ParametricPlot[ bez[{p[[1]], {0.75, 0.75 }, p[[-1]]}][u], {u, 0, 1}], ListPlot[p]}]  and it looks fine! – Ulrich Neumann Oct 31 '19 at 8:49
• Your bez is polluted. Check DownValues@bez, then Clear@bez and retry. – xzczd Oct 31 '19 at 8:52
• That means the conditional argument check is wrong? If I define Clear[bez]; bez[k_] := BezierFunction[k] the function can be plotted but the NMinimze-error is still present. – Ulrich Neumann Oct 31 '19 at 9:03
• Yes, it's wrong, try NumericQ /@ {{1, 2}, {3, 4}}. Then, execute Sum[Norm[ bez[{p[[1]], {kx, ky}, p[[-1]]}][u[i]] - p[[i]]], {i, 2, Length[p] - 1}] outside of NMinimize, you'll see u[5], u[6] therein. – xzczd Oct 31 '19 at 10:34

The easiest way is, I think, to use the definition of a Bézier curve to define a ParametricRegion and then use RegionDistance to find out how well the curve approximates the points. Here's my suggestion. Define the points:

n = 5;
p = Table[{Cos[\[CurlyPhi]], Sin[\[CurlyPhi]]}, {\[CurlyPhi],
Join[{0}, RandomReal[{0, Pi/2}, n], {Pi/2}]}];


Definition of the curve:

bezCurve[{pt : {_, _}}] := pt &;
bezCurve[ctrlPts_?MatrixQ] := bezCurve[ctrlPts] =
With[{b1 = bezCurve[Most[ctrlPts]],
b2 = bezCurve[Rest[ctrlPts]]}, (1 - #)*b1[#] + #*b2[#] &];
bezCurve[ctrlPts_?MatrixQ, t_] := Simplify[bezCurve[ctrlPts][t]];


Define the loss function to be minimized:

ClearAll[loss];
loss[pt : {__?NumericQ}] := loss[pt] = With[{
regDist = RegionDistance[
ParametricRegion[
bezCurve[{p[[1]], pt, p[[-1]]}, t],
{{t, 0, 1}}
]
]
},
Total@regDist[p]
]


Use FindMinimum to find the solution (will take a while):

FindMinimum[{loss[{x, y}],
0 < x < 3 && 0 < y < 3}, {{x, 0.937}, {y, 0.883}},
EvaluationMonitor :> Print[{{x, y}, loss[{x, y}]}]
]

• Thanks a lot, I'll try to understand your code. I stumble first at the definitions bezCurve[{pt : {_, _}}] (?means: default argument is an empty list of length 2?) and loss[pt : {__?NumericQ}] (?means: only numerical lists?) . – Ulrich Neumann Oct 31 '19 at 19:11
• bezCurve[{pt : {_, _}}] means that it will match something of the form bezCurve[{{x, y}}] and then substitute pt on the r.h.s. with {x, y} (i.e., with the outer list stripped off). There is no default argument here; that would involve a second colon. The confusion is understandable, though, since the colon is a bit weird in Mathematica. You're correct about loss[pt : {__?NumericQ}]: it will only match things like loss[{1.2, 3.5}]. A slightly better pattern would've been loss[pt : {_?NumericQ, _?NumericQ}], which will only match a list of exactly 2 numbers. – Sjoerd Smit Oct 31 '19 at 19:27
• The loss function is really tricky!!! I'm still wondering why the use of the original BezierFunction doesn't work inside FindMinimum – Ulrich Neumann Oct 31 '19 at 19:39
• It's because of the order of the evaluation. That's why I have this loss function that won't evaluate symbolically. That ensures correct evaluation order when numerical values are substituted for x and y. – Sjoerd Smit Oct 31 '19 at 22:25

There're 2 issues here:

1. The evaluation order should be properly controlled.

2. The argument of BezierFunction should be between $$0$$ and $$1$$, so we need to add constraints to NMinimize.

The following is the fixed code, I've also adjust Method option of NMinimize a bit to obtain better result:

SeedRandom[1];
n = 5; p = Table[{Cos[φ], Sin[φ]}, {φ,
Join[{0}, RandomReal[{0, Pi/2}, n], {Pi/2}]}];
ui = Table[u[i], {i, 2, Length[p] - 1}];

Clear@bez
bez[k : {{_?NumericQ, _?NumericQ} ..}] := BezierFunction[k]

Clear@norm;
norm[point_List, point2_] := Norm[point - point2]

NMinimize[{Sum[
norm[bez[{p[[1]], {kx, ky}, p[[-1]]}][u[i]], p[[i]]], {i, 2, Length[p] - 1}],
0 <= # <= 1 & /@ ui}, Join[{kx, ky}, ui], Method -> "RandomSearch"] // AbsoluteTiming
(*
{3.00299, {0.00958796, {kx -> 0.946537, ky -> 0.938663, u[2] -> 0.838486,
u[3] -> 0.0969485, u[4] -> 0.811814, u[5] -> 0.16807, u[6] -> 0.220674}}}
*)

• Thanks for this nice solution. The two modifications of BezierFunction and Norm eliminate the problems I had. Furthermore the restrictions 0<u[i]<1 work much better than 0 u[1]<u[2]<...<1 – Ulrich Neumann Nov 1 '19 at 7:19