Find optimal BezierFunction

Question

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 08: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 08:49
Your bez is polluted. Check DownValues@bez, then Clear@bez and retry. — xzczd, Oct 31 '19 at 08: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 09: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
Thanks, there was an error in the definition of ui. I modified my question. — Ulrich Neumann, Oct 31 '19 at 10:42
Now try Clear@point; Norm[point - {1, 2}] and think about what's wrong here. — xzczd, Oct 31 '19 at 10:51
In my case , if bez[...] is evaluated to a list of size 2, the argument of Norm[bez[...],p[[i]]] is ok I think. — Ulrich Neumann, Oct 31 '19 at 10:59
Norm won't wait, NMinimize doesn't have Hold* Attributes. — xzczd, Oct 31 '19 at 11:02

Sjoerd Smit · Accepted Answer · 2019-11-01T11:18:15.027

4

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}]}]
 ]

edited Nov 01 '19 at 11:18

answered Oct 31 '19 at 17:02

Sjoerd Smit

23,370
46
75

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

score 2 · Answer 2 · answered Nov 01 '19 at 05:48

There're 2 issues here:

The evaluation order should be properly controlled.
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 01 '19 at 07:19

Find optimal BezierFunction

2 Answers2

Linked