Largest ellipse inscribed in a polygon

[markz]

Hello all. I am looking for a way to inscribe the largest possible ellipse in a convex polygon.
With some effort looking, I have found a possible approach, as:
Biggest ellipse included in a convex polygon
which has code, written in matlab, that appears to work well, based off of approach suggested by chapter 8 of Boyd & Vandenberghe’s Convex Optimization
Unfortunately I never have learned matlab, and am unsure how to convert this code to C or Python.

Has anyone coded up an approach for this, either mathematical or iterative, that will solve for the largest area ellipse, or a close approximation?

I am dealing with cases with N sides to the polygon, but typically 4 to 7 sides. I thought of the Steiner inellipse (3 sides) and Laurent’s work solving for largest rectangle, but do not see how to apply to the case of 5 or more sides

To clarify, looking for the largest area ellipsi

[Joseph_Oster]

@markz

I am looking for a way to inscribe the largest possible ellipse in a convex polygon.
[…]
Has anyone coded up an approach for this, either mathematical or iterative, that will solve for the largest area ellipse, or a close approximation?

Is that ellipse a circle?

poly_radius_2021_Aug1a.gh (9.2 KB)

[markz]

Thanks Joe, if only!

Perhaps I should have stated, 'polygon with unequal sides"

for example

the matlab code does produce a solution, but again I am limited by my ignorance of matlab

[Joseph_Oster]

@markz

Perhaps I should have stated, 'polygon with unequal sides"

Yeah, I figured that was too easy. Here’s a Galapagos model. The white group interactively defines an arbitrary “Convex Hull” polygon, the cyan group is to “Locate, Rotate, R1 and R2 for Ellipse”, all controlled by Galapagos.

poly_radius_2021_Aug1b.gh (21.9 KB)

Simple.

Not bad! Galapagos stopped with Area = 51.281883 for the given shape.

Of course, you can supply your own polygon instead of using the white group to create it.

[markz]

Thank you Joe, pretty good actually, for a quick attempt. eyeballing the results, looks like it is within less than 1% of theoretical maximum. Pretty good. Thank you.

Of course, typical design, 1,200+ polygons to run through, then I will decide that I don’t like it, so I will tweak it then run it again, and so on.

I can use your approach, but would prefer to find an algorithmic approach for easier inclusion into design flow. here is the matlib code, if anyone is fluent in its’ syntax:

% Section 8.4.1, Boyd & Vandenberghe “Convex Optimization”
% Original version by Lieven Vandenberghe
% Updated for CVX by Almir Mutapcic - Jan 2006
% (a figure is generated)
%
% We find the ellipsoid E of maximum volume that lies inside of
% a polyhedra C described by a set of linear inequalities.
%
% C = { x | a_i^T x <= b_i, i = 1,…,m } (polyhedra)
% E = { Bu + d | || u || <= 1 } (ellipsoid)
%
% This problem can be formulated as a log det maximization
% which can then be computed using the det_rootn function, ie,
% maximize log det B
% subject to || B a_i || + a_i^T d <= b, for i = 1,…,m

% problem data
n = 2;
px = [0 .5 2 3 1];
py = [0 1 1.5 .5 -.5];
m = size(px,2);
pxint = sum(px)/m; pyint = sum(py)/m;
px = [px px(1)];
py = [py py(1)];

% generate A,b
A = zeros(m,n); b = zeros(m,1);
for i=1:m
A(i, = null([px(i+1)-px(i) py(i+1)-py(i)])’;
b(i) = A(i,.5[px(i+1)+px(i); py(i+1)+py(i)];
if A(i,:)*[pxint; pyint]-b(i)>0
A(i, = -A(i,:);
b(i) = -b(i);
end
end

% formulate and solve the problem
cvx_begin
variable B(n,n) symmetric
variable d(n)
maximize( det_rootn( B ) )
subject to
for i = 1:m
norm( B*A(i,:)’, 2 ) + A(i,:)*d <= b(i);
end
cvx_end

% make the plots
noangles = 200;
angles = linspace( 0, 2 * pi, noangles );
ellipse_inner = B * [ cos(angles) ; sin(angles) ] + d * ones( 1, noangles );
ellipse_outer = 2*B * [ cos(angles) ; sin(angles) ] + d * ones( 1, noangles );

clf
plot(px,py)
hold on
plot( ellipse_inner(1,:), ellipse_inner(2,:), ‘r–’ );
plot( ellipse_outer(1,:), ellipse_outer(2,:), ‘r–’ );
axis square
axis off
hold off

[Joseph_Oster]

I modified the code a bit, dropped ‘Condition B’ because it was redundant, limited location of the ellipse plane to the untrimmed surface with X and Y slider values in the 0 to 1 range. Got a larger area value of 52.9482 vs. 51.2818, though had to restart Galapagos a couple of times. It must have gotten stuck?

poly_radius_2021_Aug1c.gh (22.4 KB)

[TomTom (Tom)]

Could you wrap the source code in ```? This doesn’t mess up the code formatting. Furthermore, the script itself is simple. Rather, the implementation of the in build-matrix operations is a problem. You can translate that by using https://numerics.mathdotnet.com/ , I’m sure there is equivalent functionality implemented. In the end, it’s a “simple” LLS problem. Minimize or maximize the delta until the error is below a threshold. In the end, what Galapagos does overcomplicated and far less efficient than by using matrix operations. Rhinocommon has a small Matrix library with the ability to inverse a matrix. Maybe this is already sufficient.
Matrix Class

[markz]

Thank you Tom. My apologies, here is source code again:

'''
 Section 8.4.1, Boyd & Vandenberghe "Convex Optimization"
% Original version by Lieven Vandenberghe
% Updated for CVX by Almir Mutapcic - Jan 2006
% (a figure is generated)
%
% We find the ellipsoid E of maximum volume that lies inside of
% a polyhedra C described by a set of linear inequalities.
%
% C = { x | a_i^T x <= b_i, i = 1,...,m } (polyhedra)
% E = { Bu + d | || u || <= 1 } (ellipsoid)
%
% This problem can be formulated as a log det maximization
% which can then be computed using the det_rootn function, ie,
%     maximize     log det B
%     subject to   || B a_i || + a_i^T d <= b,  for i = 1,...,m

% problem data
n = 2;
px = [0 .5 2 3 1];
py = [0 1 1.5 .5 -.5];
m = size(px,2);
pxint = sum(px)/m; pyint = sum(py)/m;
px = [px px(1)];
py = [py py(1)];

% generate A,b
A = zeros(m,n); b = zeros(m,1);
for i=1:m
  A(i,:) = null([px(i+1)-px(i) py(i+1)-py(i)])';
  b(i) = A(i,:)*.5*[px(i+1)+px(i); py(i+1)+py(i)];
  if A(i,:)*[pxint; pyint]-b(i)>0
    A(i,:) = -A(i,:);
    b(i) = -b(i);
  end
end

% formulate and solve the problem
cvx_begin
    variable B(n,n) symmetric
    variable d(n)
    maximize( det_rootn( B ) )
    subject to
       for i = 1:m
           norm( B*A(i,:)', 2 ) + A(i,:)*d <= b(i);
       end
cvx_end

% make the plots
noangles = 200;
angles   = linspace( 0, 2 * pi, noangles );
ellipse_inner  = B * [ cos(angles) ; sin(angles) ] + d * ones( 1, noangles );
ellipse_outer  = 2*B * [ cos(angles) ; sin(angles) ] + d * ones( 1, noangles );

clf
plot(px,py)
hold on
plot( ellipse_inner(1,:), ellipse_inner(2,:), 'r--' );
plot( ellipse_outer(1,:), ellipse_outer(2,:), 'r--' );
axis square
axis off
hold off
'''

[markz]

That’s great Joseph, thank you. Example looks like it is very close to theoretical maximum, again by eyeballing, probably within 0.2% of maximum. I have not used galapagos very much. This will work for a single case, although I have not tried it, will it crash if I feed it a tree structure? Another idea might be to extract optimal case, then adjust ellipse, as is, so it is equidistant from polygon, i.e., distributing the error, then do small offsets until the polygon is intersected again. If the error is small, it would work for my purposes. Not good enough for engineering work where exact contact/tangency is needed, but still pretty good.
Thank you for your input! I will be tinkering with this throughout the day

[Joseph_Oster]

@markz

A(i, = null([px(i+1)-px(i) py(i+1)-py(i)])’;
b(i) = A(i, .5 [px(i+1)+px(i); py(i+1)+py(i)];
if A(i,:)*[pxint; pyint]-b(i)>0
A(i, = -A(i,:);

Try wrapping that code as ‘Preformatted text’ which will avoid the emoticons. Instead of ‘Blockquote’.

[markz]

Blockquote

% generate A,b
A = zeros(m,n); b = zeros(m,1);
for i=1:m
  A(i,:) = null([px(i+1)-px(i) py(i+1)-py(i)])';
  b(i) = A(i,:)*.5*[px(i+1)+px(i); py(i+1)+py(i)];
  if A(i,:)*[pxint; pyint]-b(i)>0
    A(i,:) = -A(i,:);
    b(i) = -b(i);
  end
end

Blockquote

[markz]

thankyou Joe

[Joseph_Oster]

Section 8.4.1, Boyd & Vandenberghe “Convex Optimization”
% Original version by Lieven Vandenberghe
% Updated for CVX by Almir Mutapcic - Jan 2006
% (a figure is generated)
%
% We find the ellipsoid E of maximum volume that lies inside of
% a polyhedra C described by a set of linear inequalities.
%
% C = { x | a_i^T x <= b_i, i = 1,…,m } (polyhedra)
% E = { Bu + d | || u || <= 1 } (ellipsoid)
%
% This problem can be formulated as a log det maximization
% which can then be computed using the det_rootn function, ie,
% maximize log det B
% subject to || B a_i || + a_i^T d <= b, for i = 1,…,m

% problem data
n = 2;
px = [0 .5 2 3 1];
py = [0 1 1.5 .5 -.5];
m = size(px,2);
pxint = sum(px)/m; pyint = sum(py)/m;
px = [px px(1)];
py = [py py(1)];

% generate A,b
A = zeros(m,n); b = zeros(m,1);
for i=1:m
A(i,:slight_smile: = null([px(i+1)-px(i) py(i+1)-py(i)])’;
b(i) = A(i,:slight_smile: *.5* [px(i+1)+px(i); py(i+1)+py(i)];
if A(i,:)*[pxint; pyint]-b(i)>0
A(i,:slight_smile: = -A(i,:);
b(i) = -b(i);
end
end

% formulate and solve the problem
cvx_begin
variable B(n,n) symmetric
variable d(n)
maximize( det_rootn( B ) )
subject to
for i = 1:m
norm( B*A(i,:)’, 2 ) + A(i,:)*d <= b(i);
end
cvx_end

% make the plots
noangles = 200;
angles = linspace( 0, 2 * pi, noangles );
ellipse_inner = B * [ cos(angles) ; sin(angles) ] + d * ones( 1, noangles );
ellipse_outer = 2*B * [ cos(angles) ; sin(angles) ] + d * ones( 1, noangles );

clf
plot(px,py)
hold on
plot( ellipse_inner(1,:), ellipse_inner(2,:), ‘r–’ );
plot( ellipse_outer(1,:), ellipse_outer(2,:), ‘r–’ );
axis square
axis off
hold off

[AndersDeleuran]

See this post for how to define proper syntax highlighting:

Category definition for Scripting

Topics related to RhinoScript and Python scripting. Syntax Highlighting: Format and highlight your script code by fencing your code block within three back ticks. If you name the language your script is in, you get syntax highlighting, too: Python: ```python # Your python code here print “Hello Pythonic World” ``` becomes #Your python code here print "Hello Pythonic World" VBScript: ```vbscript ’ Your VBScript code here print “Hello VBScript” ``` becomes ' Your VBScript code he…

[markz]

Thankyou Anders. I know, must be irritating seeing the same basic mistakes repeatedly, at least we can hope the same mistake is made by different people rather that one person doing all of the repetition, I will try to stay out of that category.
Now, onto making this approach work. Porting it to GH
Picking the code apart, '%problem data, I think that is the set of input points, which the GH component will get from its input
the ‘clf’ block at the end will be replaced by the GH component outputs
now, I need to convert from matlib to python
I am thinking that, in the comments, although the author is calling it an ellipsoid, I think this is still a 2 dimensional case
for those reading this in the future, I am attaching again what I hope is a clean view of the source code

% Section 8.4.1, Boyd & Vandenberghe "Convex Optimization"
% Original version by Lieven Vandenberghe
% Updated for CVX by Almir Mutapcic - Jan 2006
% (a figure is generated)
%
% We find the ellipsoid E of maximum volume that lies inside of
% a polyhedra C described by a set of linear inequalities.
%
% C = { x | a_i^T x <= b_i, i = 1,...,m } (polyhedra)
% E = { Bu + d | || u || <= 1 } (ellipsoid)
%
% This problem can be formulated as a log det maximization
% which can then be computed using the det_rootn function, ie,
%     maximize     log det B
%     subject to   || B a_i || + a_i^T d <= b,  for i = 1,...,m

% problem data
n = 2;
px = [0 .5 2 3 1];
py = [0 1 1.5 .5 -.5];
m = size(px,2);
pxint = sum(px)/m; pyint = sum(py)/m;
px = [px px(1)];
py = [py py(1)];

% generate A,b
A = zeros(m,n); b = zeros(m,1);
for i=1:m
  A(i,:) = null([px(i+1)-px(i) py(i+1)-py(i)])';
  b(i) = A(i,:)*.5*[px(i+1)+px(i); py(i+1)+py(i)];
  if A(i,:)*[pxint; pyint]-b(i)>0
    A(i,:) = -A(i,:);
    b(i) = -b(i);
  end
end

% formulate and solve the problem
cvx_begin
    variable B(n,n) symmetric
    variable d(n)
    maximize( det_rootn( B ) )
    subject to
       for i = 1:m
           norm( B*A(i,:)', 2 ) + A(i,:)*d <= b(i);
       end
cvx_end

% make the plots
noangles = 200;
angles   = linspace( 0, 2 * pi, noangles );
ellipse_inner  = B * [ cos(angles) ; sin(angles) ] + d * ones( 1, noangles );
ellipse_outer  = 2*B * [ cos(angles) ; sin(angles) ] + d * ones( 1, noangles );

clf
plot(px,py)
hold on
plot( ellipse_inner(1,:), ellipse_inner(2,:), 'r--' );
plot( ellipse_outer(1,:), ellipse_outer(2,:), 'r--' );
axis square
axis off
hold off

[Joseph_Oster]

Additional code changes:

• Restored ‘Condition B’ to make sure the center of the ellipse is inside the polygon. It is required after all, otherwise the ellipse can be valid outside the polygon.

• Extended the range of the rotation angle slider from 90 degrees to 180 degrees.

• Made the R1 and R2 sliders (ellipse radii) parametric (0 to 1) instead of absolute. They are multiplied by the diagonal length of the polygon’s bounding box.

Galapagos has a tendency to hang in certain starting conditions, requiring the GH window to be closed to regain control. I don’t know Galapagos very well, maybe it’s possible to adjust options to avoid that?

P.S. The R1 and R2 sliders should be multiplied by half the diagonal length of the polygon’s bounding box, but even so, I can’t seem to avoid Galapagos quickly hanging when the starting condition is a failure (the ellipse intersects the polygon or is entirely outside of it).

Even worse, solutions are not consistent and again seem to depend too much on the starting condition. Maybe if my test conditions (B, C & D) were fractional instead of Boolean? This reminds me why I gave up on Galapagos long ago,

poly_radius_2021_Aug2b.gh (28.0 KB)

[Joseph_Oster]

Extraordinary Conics: The Most Difficult Math Problem I Ever Solved

[markz]

hmmm. or perhaps if the test conditions adjust for the error in between the ellipse and the polygon, i.e., calc distances between each polygon line segment to the ellipse curve, average them, and then move center point of ellipse to correspond, and run it again.
I haven’t played much with Galapagos, I always like to see a calculated solution, but your initial approach seemed to work pretty well Joe, brute force combined with reasoned guesses.

[Joseph_Oster]

@markz

perhaps if the test conditions adjust for the error in between the ellipse and the polygon, i.e., calc distances between each polygon line segment to the ellipse curve, average them, and then move center point of ellipse to correspond, and run it again.

I don’t know how to make that work, given the “Fitness” criteria I’m using, but you are welcome to try.

I did manage to avoid the Boolean failure condition and freeze-ups I was seeing. The “Fitness” value being used by Galapagos is the area of the ellipse, which Galapagos tries to maximize. If any of three “violation” conditions were encountered, the area (“Fitness”) value was set to zero in earlier code.

Instead, I now sum up the violations as follows:

• whether or not the centroid of the ellipse is inside the polygon (‘Condition B’, 20 points if not)

• number of polygon vertices inside the ellipse (‘Condition C’, 1 point each)

• number of intersection points between the polygon and the ellipse (‘Condition D’, 1 point each)

If the total of these violation points is zero, all is well and the area of the ellipse is the fitness value. Otherwise the area is multiplied by the negative total points. So far, this appears to work better.

The ‘Evolutionary Solver’ is still not reliably repeatable and it sometimes terminates too soon (without hanging, it just stops looping). The ‘Annealing Solver’ appears to never stop?

poly_radius_2021_Aug2c.gh (28.2 KB)

[markz]

On a related subject, how to find the largest rectangle in a polygon.
There is a new paper published July 5, 2021, so I have yet to digest it fully, but here is the link:
Finding Largest Rectangles in Convex Polygons∗

1405.1223.pdf

Laurent Delrieu has done some excellent work in this area.