Skip to main content

graphics - Why BoundingRegion does not work well?


Bug introduced in 8 or earlier and fixed in 12.0.0


The underlying bug is unstable and wrong rendering of Disk[] and Circle[] primitives after applying GeometricTransformation, see answer below.



I have a set of data points describing an ellipse in the plane. I want to obtain the best ellipse that fits them.



As a first attempt, I use this Q&A, and went well. However, I used in the past the following:


param = NArgMin[{Norm[
Function[{x, y}, ((x - h)*Cos[\[Alpha]] - (y - k)*Sin[\[Alpha]])^2/a^2 + 
((x - h)*Sin[\[Alpha]] + (y - k)*Cos[\[Alpha]])^2/b^2 - 1] @@@ elip]}, 
{a, b, h, k, \[Alpha]}]

However, now this code does not work, and I do not know why.


Then, I went over the function BoundingRegion, and run:


elipse = BoundingRegion[elip, "FastEllipse"]
Graphics[{{LightBlue, elipse}, Point[elip]}, ImageSize -> Medium, 

Axes -> True, PlotRange -> {{1, 578}, {1, 724}}]

and I got


enter image description here


Why do my two last attempts fail? Further, I do not understand why of the result by means of BoundingRegion as I should get the best ellipse that contains the points, shouldn't I?


My dataset is:


elip={{238., 277.}, {238., 278.}, {238., 279.}, {238., 280.}, {238., 
281.}, {238., 282.}, {238., 283.}, {238., 284.}, {238., 
285.}, {238., 286.}, {238., 287.}, {238., 288.}, {238., 
289.}, {238., 290.}, {238., 291.}, {238., 292.}, {238., 

293.}, {238., 294.}, {238., 295.}, {238., 296.}, {238., 
297.}, {238., 298.}, {238., 299.}, {238., 300.}, {238., 
301.}, {238., 302.}, {238., 303.}, {238., 304.}, {238., 
305.}, {238., 306.}, {238., 307.}, {238., 308.}, {238., 
309.}, {238., 310.}, {239., 271.}, {239., 272.}, {239., 
273.}, {239., 274.}, {239., 275.}, {239., 313.}, {239., 
314.}, {239., 315.}, {239., 316.}, {239., 317.}, {239., 
318.}, {239., 319.}, {240., 266.}, {240., 267.}, {240., 
268.}, {240., 269.}, {240., 321.}, {240., 322.}, {240., 
323.}, {240., 324.}, {240., 325.}, {240., 326.}, {241., 

263.}, {241., 264.}, {241., 265.}, {241., 327.}, {241., 
328.}, {241., 329.}, {241., 330.}, {241., 331.}, {242., 
260.}, {242., 261.}, {242., 262.}, {242., 333.}, {242., 
334.}, {242., 335.}, {242., 336.}, {243., 258.}, {243., 
259.}, {243., 338.}, {243., 339.}, {243., 340.}, {243., 
341.}, {244., 256.}, {244., 257.}, {244., 342.}, {244., 
343.}, {244., 344.}, {244., 345.}, {245., 254.}, {245., 
255.}, {245., 346.}, {245., 347.}, {245., 348.}, {246., 
253.}, {246., 254.}, {246., 350.}, {246., 351.}, {246., 
352.}, {247., 251.}, {247., 252.}, {247., 353.}, {247., 

354.}, {247., 355.}, {248., 250.}, {248., 251.}, {248., 
356.}, {248., 357.}, {248., 358.}, {249., 249.}, {249., 
250.}, {249., 359.}, {249., 360.}, {249., 361.}, {250., 
248.}, {250., 249.}, {250., 362.}, {250., 363.}, {250., 
364.}, {251., 247.}, {251., 248.}, {251., 365.}, {251., 
366.}, {251., 367.}, {252., 247.}, {252., 368.}, {252., 
369.}, {253., 246.}, {253., 370.}, {253., 371.}, {253., 
372.}, {254., 245.}, {254., 373.}, {254., 374.}, {255., 
245.}, {255., 375.}, {255., 376.}, {256., 244.}, {256., 
377.}, {256., 378.}, {256., 379.}, {257., 244.}, {257., 

380.}, {257., 381.}, {258., 243.}, {258., 382.}, {258., 
383.}, {259., 243.}, {259., 384.}, {259., 385.}, {260., 
243.}, {260., 386.}, {260., 387.}, {261., 243.}, {261., 
388.}, {261., 389.}, {262., 242.}, {262., 390.}, {262., 
391.}, {263., 242.}, {263., 392.}, {263., 393.}, {264., 
242.}, {264., 394.}, {264., 395.}, {265., 242.}, {265., 
396.}, {265., 397.}, {266., 242.}, {266., 397.}, {266., 
398.}, {267., 242.}, {267., 399.}, {267., 400.}, {268., 
242.}, {268., 401.}, {268., 402.}, {269., 242.}, {269., 
403.}, {269., 404.}, {270., 242.}, {270., 404.}, {270., 

405.}, {271., 242.}, {271., 406.}, {271., 407.}, {272., 
242.}, {272., 408.}, {272., 409.}, {273., 242.}, {273., 
409.}, {273., 410.}, {274., 411.}, {274., 412.}, {275., 
243.}, {275., 412.}, {275., 413.}, {276., 243.}, {276., 
414.}, {276., 415.}, {277., 243.}, {277., 416.}, {278., 
243.}, {278., 417.}, {278., 418.}, {279., 244.}, {279., 
419.}, {280., 244.}, {280., 420.}, {280., 421.}, {281., 
244.}, {281., 421.}, {281., 422.}, {282., 245.}, {282., 
423.}, {282., 424.}, {283., 245.}, {283., 424.}, {283., 
425.}, {284., 246.}, {284., 426.}, {285., 246.}, {285., 

427.}, {285., 428.}, {286., 247.}, {286., 428.}, {286., 
429.}, {287., 247.}, {287., 429.}, {287., 430.}, {288., 
248.}, {288., 431.}, {289., 248.}, {289., 249.}, {289., 
432.}, {289., 433.}, {290., 249.}, {290., 433.}, {290., 
434.}, {291., 250.}, {291., 434.}, {291., 435.}, {292., 
250.}, {292., 251.}, {292., 435.}, {292., 436.}, {293., 
251.}, {293., 437.}, {294., 252.}, {294., 438.}, {295., 
253.}, {295., 439.}, {296., 253.}, {296., 254.}, {296., 
440.}, {297., 254.}, {297., 255.}, {297., 441.}, {298., 
255.}, {298., 442.}, {299., 256.}, {299., 443.}, {300., 

257.}, {300., 444.}, {301., 258.}, {301., 444.}, {301., 
445.}, {302., 259.}, {302., 445.}, {302., 446.}, {303., 
260.}, {303., 446.}, {304., 261.}, {304., 447.}, {305., 
262.}, {305., 448.}, {306., 263.}, {306., 448.}, {306., 
449.}, {307., 264.}, {307., 265.}, {307., 449.}, {308., 
265.}, {308., 266.}, {308., 450.}, {309., 266.}, {309., 
267.}, {309., 450.}, {309., 451.}, {310., 268.}, {310., 
451.}, {311., 269.}, {311., 452.}, {312., 270.}, {312., 
271.}, {312., 452.}, {313., 271.}, {313., 272.}, {313., 
453.}, {314., 273.}, {314., 453.}, {315., 274.}, {315., 

275.}, {315., 454.}, {316., 275.}, {316., 276.}, {316., 
454.}, {317., 277.}, {317., 278.}, {317., 455.}, {318., 
278.}, {318., 279.}, {318., 455.}, {319., 280.}, {319., 
281.}, {319., 456.}, {320., 281.}, {320., 282.}, {320., 
456.}, {321., 283.}, {321., 284.}, {321., 456.}, {322., 
284.}, {322., 285.}, {322., 456.}, {323., 286.}, {323., 
287.}, {323., 457.}, {324., 288.}, {324., 289.}, {324., 
457.}, {325., 289.}, {325., 290.}, {325., 457.}, {326., 
291.}, {326., 292.}, {326., 457.}, {327., 293.}, {327., 
294.}, {327., 457.}, {328., 294.}, {328., 295.}, {328., 

296.}, {328., 457.}, {329., 296.}, {329., 297.}, {329., 
457.}, {330., 298.}, {330., 299.}, {330., 457.}, {331., 
300.}, {331., 301.}, {331., 457.}, {332., 302.}, {332., 
303.}, {332., 457.}, {333., 304.}, {333., 305.}, {333., 
457.}, {334., 306.}, {334., 307.}, {334., 457.}, {335., 
308.}, {335., 309.}, {335., 457.}, {336., 310.}, {336., 
311.}, {336., 457.}, {337., 312.}, {337., 313.}, {337., 
457.}, {338., 314.}, {338., 315.}, {338., 456.}, {339., 
316.}, {339., 317.}, {339., 456.}, {340., 318.}, {340., 
319.}, {340., 456.}, {341., 320.}, {341., 321.}, {341., 

322.}, {341., 456.}, {342., 322.}, {342., 323.}, {342., 
324.}, {342., 455.}, {343., 325.}, {343., 326.}, {343., 
455.}, {344., 327.}, {344., 328.}, {344., 329.}, {344., 
454.}, {345., 330.}, {345., 331.}, {345., 454.}, {346., 
332.}, {346., 333.}, {346., 334.}, {346., 453.}, {347., 
335.}, {347., 336.}, {347., 452.}, {347., 453.}, {348., 
337.}, {348., 338.}, {348., 339.}, {348., 452.}, {349., 
340.}, {349., 341.}, {349., 342.}, {349., 451.}, {350., 
343.}, {350., 344.}, {350., 345.}, {350., 450.}, {351., 
346.}, {351., 347.}, {351., 348.}, {351., 449.}, {352., 

349.}, {352., 350.}, {352., 351.}, {352., 447.}, {352., 
448.}, {353., 352.}, {353., 353.}, {353., 354.}, {353., 
446.}, {353., 447.}, {354., 356.}, {354., 357.}, {354., 
358.}, {354., 445.}, {354., 446.}, {355., 359.}, {355., 
360.}, {355., 361.}, {355., 362.}, {355., 443.}, {355., 
444.}, {356., 363.}, {356., 364.}, {356., 365.}, {356., 
366.}, {356., 441.}, {356., 442.}, {357., 367.}, {357., 
368.}, {357., 369.}, {357., 370.}, {357., 371.}, {357., 
439.}, {357., 440.}, {358., 372.}, {358., 373.}, {358., 
374.}, {358., 375.}, {358., 376.}, {358., 436.}, {358., 

437.}, {358., 438.}, {359., 377.}, {359., 378.}, {359., 
379.}, {359., 380.}, {359., 381.}, {359., 432.}, {359., 
433.}, {359., 434.}, {359., 435.}, {360., 383.}, {360., 
384.}, {360., 385.}, {360., 386.}, {360., 387.}, {360., 
388.}, {360., 389.}, {360., 427.}, {360., 428.}, {360., 
429.}, {360., 430.}, {360., 431.}, {361., 391.}, {361., 
392.}, {361., 393.}, {361., 394.}, {361., 395.}, {361., 
396.}, {361., 397.}, {361., 398.}, {361., 399.}, {361., 
400.}, {361., 401.}, {361., 419.}, {361., 420.}, {361., 
421.}, {361., 422.}, {361., 423.}, {361., 424.}, {361., 

425.}, {361., 426.}, {362., 408.}, {362., 409.}, {362., 
410.}, {362., 411.}, {362., 412.}, {362., 413.}}

Thanks for your time.



Answer



Here is an answer to your first question: Using Norm[] you'll get Abs[]-terms in the functional, which are sometimes problematical. Using the sum of squares


opt = FindMinimum[{#.# &[
Apply[Function[{x, 
   y}, ((x - h)*Cos[\[Alpha]] - (y - k)*Sin[\[Alpha]])^2/
    a^2 + ( (x - h)*Sin[\[Alpha]] + (y - k)*Cos[\[Alpha]])^2/

    b^2 - 1], elip, 1]]
, -Pi <= \[Alpha] <= Pi, a > b   },
{ a , b , {h, Mean[elip][[1]]}, {k, Mean[elip][[2]]} , \[Alpha] },MaxIterations -> 1000, AccuracyGoal -> 4, PrecisionGoal -> 5]
(* {0.131336, {a -> 114.631, b -> 50.1975, h -> 299.194,k -> 350.063, \[Alpha] -> 1.93331}}*)

gives this result


Show[{ContourPlot[(((x - h)*Cos[\[Alpha]] - (y - k)*Sin[\[Alpha]])^2/
    a^2 + ( (x - h)*Sin[\[Alpha]] + (y - k)*Cos[\[Alpha]])^2/
    b^2 - 1 /. opt[[2]]) == 0, {x, 200, 400}, {y, 200, 500}]
,Graphics[{Red,Point[elip] }]},PlotRange -> All]


enter image description here


for the approximation.


Comments

Post a Comment

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more