Minimum enclosing circle/ellipse 2

In [1] where I discussed how to find the minimum enclosing circle and minimum enclosing ellipse around a set of points. This is a follow-up post where I extend this to sets of circles and ellipses. 

1. MINIMUM ENCLOSING CIRCLE

Here our data is a set of \(n\) circles (or disks) of different size. We want to find the smallest circle that contains all these circles. 

An example data set with random circles can look like:

----     30 PARAMETER circles  coordinates of center and radius

                   x           y           r

circle1        4.294      21.082       1.663
circle2       13.759       7.528       4.014
circle3        7.305       5.601       1.961
circle4        8.746      21.407       6.235
circle5        1.678      12.505       2.591
circle6       24.953      14.468       2.715
circle7       24.778      19.056       4.564
circle8        3.267      15.993       5.336
circle9        3.988       6.252       4.769
circle10      16.723      10.884       3.783
circle11       8.993       8.786       3.480
circle12       3.287       3.753       1.706
circle13      14.728      20.772       2.885
circle14       5.770      16.643       1.279
circle15      19.396       7.591       3.031

The task is to find the smallest circle around this:

The pink colored circles are completely inside another circle, so we can drop them by preprocessing the data. So, how do we detect if a circle is inside another one?

If we do this test for each pair of circles \((i,j)\) with \({\color{darkblue}r}_i \le {\color{darkblue}r}_j\), then we get:

----     49 PARAMETER circles  coordinates of center and radius

                   x           y           r     candrop

circle1        4.294      21.082       1.663       1.000
circle2       13.759       7.528       4.014
circle3        7.305       5.601       1.961
circle4        8.746      21.407       6.235
circle5        1.678      12.505       2.591
circle6       24.953      14.468       2.715
circle7       24.778      19.056       4.564
circle8        3.267      15.993       5.336
circle9        3.988       6.252       4.769
circle10      16.723      10.884       3.783
circle11       8.993       8.786       3.480
circle12       3.287       3.753       1.706       1.000
circle13      14.728      20.772       2.885
circle14       5.770      16.643       1.279       1.000
circle15      19.396       7.591       3.031

We can recognize the three red circles here. This is of course not as powerful as a convex hull, but it may be better than nothing. In addition, this step forces us to think about the conditions for a disk being completely inside another circle. We can use this as a constraint in the following models.

We can form the NLP model to find the smallest enclosing circle with center coordinates \({\color{darkred}C}\) and radius \({\color{darkred}R}\) using the same mechanism: 

NLP formulation 1
\[\begin{aligned}\min_{{\color{darkred}C},{\color{darkred}R}}\>&{\color{darkred}R}\\ & \sqrt{{\color{darkblue}\varepsilon} + \sum_k\left({\color{darkred}C}_k-{\color{darkblue}c}_{i,k}\right)^2} + {\color{darkblue}r}_i \le {\color{darkred}R} && \forall i \\ & k \in \{x,y\} \\& i \in \{1,\dots,n\} \\ & {\color{darkblue}\varepsilon} = 0.0001 \end{aligned}\]

I added a small constant \({\color{darkblue}\varepsilon}\gt 0\) to prevent the argument of the square root to become zero (the gradient is not defined there). A different formulation can look like:

NLP formulation 2

\[\begin{aligned}\min_{{\color{darkred}C},{\color{darkred}R}}\>&{\color{darkred}R}\\ &  \sum_k\left({\color{darkred}C}_k-{\color{darkblue}c}_{i,k}\right)^2  \le \left({\color{darkred}R}- {\color{darkblue}r}_i \right)^2 && \forall i \\ & {\color{darkred}R} \ge \max_i {\color{darkblue}r}_i \\& {\color{darkred}C}_k^0 := \frac{1}{n} \sum_i {\color{darkblue}c}_{i,k} \\ &{\color{darkred}R}^0 := \max_i \left[ \sqrt{ \sum_k\left({\color{darkred}C}_k^0-{\color{darkblue}c}_{i,k}\right)^2} +{\color{darkblue}r}_i \right]\\ & k \in \{x,y\}  \\& i \in \{1,\dots,n\}  \end{aligned}\]

I removed the square root. Now we need to protect against \({\color{darkred}R}\) becoming negative. Otherwise, a very large negative value would be very attractive. We could have used a lower bound of 0, but I used here that the enclosing circle must be at least as large as the largest inner circle. I also added the calculation of a good starting point \(({\color{darkred}C}_k^0, {\color{darkred}R}^0)\).

We can reformulate our NLP problem as a SOCP (Second-Order Conic Program) problem. It looks like:

SOCP formulation

\[\begin{aligned}\min\>&{\color{darkred}R} \\& {\color{darkred}\Gamma}_{i,k} = {\color{darkred}C}_k-{\color{darkblue}c}_{i,k} && \forall i,k \\ & {\color{darkred}\Delta}_{i} = {\color{darkred}R} - {\color{darkblue}r}_{i} && \forall i \\ & \sum_k {\color{darkred}\Gamma}_{i,k}^2 \le {\color{darkred}\Delta}_{i}^2 && \forall i \\ & {\color{darkred}\Gamma}_{i,k}\>{\bf{free}} \\ & {\color{darkred}\Delta}_{i} \ge 0 \\ & k \in \{x,y\} \\& i \in \{1,\dots,n\} \end{aligned}\]

We needed to introduce a number of additional variables in order to follow the strict rules for specifying SOCP constraints. Like for an LP, there is no need for an initial point. There is also no need to impose a lower bound on \({\color{darkred}R}\) (that is taken care of by \({\color{darkred}\Delta}_{i}\ge 0\)). 

All these models solve very easily and need very few iterations. We get identical solutions:

----    129 PARAMETER results  

               x           y           r

init      12.360      12.570      18.574
NLP       14.296      12.600      16.875
SOCP      14.296      12.600      16.875

In some cases, I found that NLP (using the active set solver CONOPT) was a little bit more precise than the Cplex SOCP solver. This may be the natural consequence of using an interior point code. Of course, we can always use the Cplex solution as starting point for CONOPT. This makes CONOPT essentially a crossover algorithm for Cplex.

 

GAMS Model

$ontext     Given n circles or disks (with different radii), find the smallest   enclosing circle.     We use a NLP model and a SOCP model   $offtext     *--------------------------------------------------------------------- * options *---------------------------------------------------------------------   * 0: no plot * 1: produce html plot $set htmlplot 1   * 0: don't remove circles contained in another circle * 1: do remove them $set drop 1   option nlp=conopt, qcp = cplex;   *--------------------------------------------------------------------- * data *---------------------------------------------------------------------   sets   i 'inner circles' /circle1*circle15/   a 'attributes for storing results' /x,y,r/   k(a) 'coordinate' /x,y/ ;   alias (i,j);   Parameters  circles(*,*)   'coordinates of center and radius' ;   circles(i,k) = uniform(0,25); circles(i,'r') = uniform(1,7); display circles;     *--------------------------------------------------------------------- * preprocessing * remove circles inside another one *--------------------------------------------------------------------- scalar d 'distance between centers';   set ne(i,j) ' not the same'; ne(i,j) = ord(i)<>ord(j);   circles(i,'candrop')=0; loop(ne(i,j)$(circles(i,'r')<=circles(j,'r') and circles(i,'candrop')=0),      d = sqrt(sum(k, sqr(circles(i,k)-circles(j,k))));    circles(i,'candrop')$(d + circles(i,'r') <= circles(j,'r')) = 1;   ); display circles;   set keep(i) 'circle is not dropped from data set'; keep(i) = yes; keep(i)$%drop% = circles(i,'candrop')=0; display keep;   *--------------------------------------------------------------------- * NLP Model *---------------------------------------------------------------------   variables    c(k)   'center'    r      'radius' ;   r.lo = smax(keep(i),circles(i,'r'));   * * initial point calculation * parameters    dist(i) 'distance'    init(a) 'initial values' ; init(k) = sum(keep(i),circles(i,k))/card(keep); dist(keep(i)) = sqrt(sum(k, sqr(init(k)-circles(i,k)))); init('r') = smax(keep(i), dist(i) + circles(i,'r')); c.l(k) = init(k); r.l = init('r');   Equations    inside(i) ;   inside(keep(i))..   sum(k,sqr(c(k)-circles(i,k))) =l= sqr(r-circles(i,'r'));   model m1 /all/; solve m1 minimizing r using nlp;   parameter results(*,*); results('init',a) = init(a); results('NLP',k) = c.l(k); results('NLP','r') = r.l;   display results;   *--------------------------------------------------------------------- * SOCP Model *---------------------------------------------------------------------   variables    diffc(i,k) 'difference between outer center and inner centers'    diffr(i) 'difference between outer radius and inner radii' ; positive variable diffr;   * no need for bound on r r.lo = -INF;   Equations    calcdiffc(i,k)    calcdiffr(i)    inside2(i) ;   calcdiffc(keep(i),k).. diffc(i,k) =e= c(k)-circles(i,k); calcdiffr(keep(i))..   diffr(i) =e= r-circles(i,'r'); inside2(keep(i))..     sum(k,sqr(diffc(i,k))) =l= sqr(diffr(i));   model m2 /all-m1/; solve m2 minimizing r using qcp;   results('SOCP',k) = c.l(k); results('SOCP','r') = r.l;   display results;     *--------------------------------------------------------------------- * visualization *---------------------------------------------------------------------   $set html  plot.html $set data  data.js     $if %htmlplot%==0 $goto skipplot   file fdata /%data%/; put fdata;     * plot first not dropped (large) circles * then plot dropped circles. set    sd /keep,drop,outer/    sdi(sd,*)   ; sdi('keep',i)$(not circles(i,'candrop')) = yes; sdi('drop',i)$circles(i,'candrop') = yes; sdi('outer','optimal') = yes; sdi('outer','initial') = yes;   circles('initial','r') = init('r'); circles('initial',k) = init(k); circles('optimal','r') = r.l; circles('optimal',k) = c.l(k); display sdi,circles;   Scalars    firstred /1/    firstblue /1/ ;   alias (*,item);   * circles put "circles = ["/; loop(sdi(sd,item),   put " {"/;   put "   type:'circle',"/;   put "   xref:'x',"/;   put "   yref:'y',"/;   put "   x0:",(circles(item,'x')-circles(item,'r')):0:4,","/;   put "   y0:",(circles(item,'y')-circles(item,'r')):0:4,","/;   put "   x1:",(circles(item,'x')+circles(item,'r')):0:4,","/;   put "   y1:",(circles(item,'y')+circles(item,'r')):0:4,","/;   if (sameas(item,'optimal'),      put "   line:{color:'Black'},"/;      put "   name:'optimal solution',"/;      put "   showlegend:true,"/;   elseif sameas(item,'initial'),      put "   line:{color:'Grey'},"/;      put "   name:'initial solution',"/;      put "   showlegend:true,"/;   elseif circles(item,'candrop'),      put "   line:{color:'Red'},"/;      put "   fillcolor:'Pink',"/;      put$firstred "   name:'dropped',"/;      put$firstred "   showlegend:true,"/;      firstred = 0;   else      put "   line:{color:'Blue'},"/;      put "   fillcolor:'LightSkyBlue',"/;      put$firstblue "   name:'data',"/;      put$firstblue "   showlegend:true,"/;      firstblue = 0;   );   put "   opacity:0.5,"/;    put " },"/; ); put "];"/;   scalars xmin,xmax,ymin,ymax,range0,range1; xmin = min(init('x')-init('r'),c.l('x')-r.l); xmax = max(init('x')+init('r'),c.l('x')+r.l); ymin = min(init('y')-init('r'),c.l('y')-r.l); ymax = max(init('y')+init('r'),c.l('y')+r.l); range0 = floor(min(xmin,ymin)); range1 = ceil(max(xmax,ymax));   put "range0=",range0:0:4,";"/; put "range1=",range1:0:4,";"/;   putclose;   $onecho > %html% <html> <script src="https://cdn.plot.ly/plotly-3.4.0.min.js" charset="utf-8"></script> <script src="%data%" charset="utf-8"></script> <h1>Smallest Encompassing Circle</h1> <div id="plotDiv"></div> <script>   var layout = {   autosize: false,   width: 700,   height: 650,   showlegend: true,   shapes: circles,   xaxis: { range: [range0,range1] },   yaxis: { range: [range0,range1] } } var options = {staticPlot: true, displayModeBar: false}; Plotly.newPlot('plotDiv', [], layout, options); </script> </html> $offecho   executetool 'win32.ShellExecute "%html%"';   $label skipplot

 

Conopt/Cplex log

NLP Model
=========

 
    C O N O P T   version 4.39.0
    Copyright (C) GAMS Software GmbH
                  GAMS Development Corporation
 
    Will use up to 4 threads.
 
 
    The user model has 15 constraints and 3 variables
    with 45 Jacobian elements, 45 of which are nonlinear.
    The Hessian of the Lagrangian has 3 elements on the diagonal,
    0 elements below the diagonal, and 3 nonlinear variables.
 
    Iter Phase   Ninf   Infeasibility   RGmax      NSB   Step  InItr MX OK
       0   0          0.0000000000E+00 (Input point)
 
    Preprocessed model has 15 constraints and 3 variables
    with 45 Jacobian elements, 45 of which are nonlinear.
 
    Iter Phase   Ninf   Infeasibility   RGmax      NSB   Step  InItr MX OK
                      0.0000000000E+00 (Full preprocessed model)
                      0.0000000000E+00 (After scaling)
 
 ** Feasible solution. Value of objective =    19.8898036129
 
    Iter Phase   Ninf     Objective     RGmax      NSB   Step  InItr MX OK
       1   3          1.6875101130E+01 7.3E+00       3 1.0E+00     2 T  T
       4   4          1.6875101115E+01 1.9E-11       1
 
 ** Optimal solution. Reduced gradient less than tolerance.


SOCP Model
==========

IBM ILOG CPLEX   53.1.0 6a03d3b9 Feb 16, 2026          WEI x86 64bit/MS Window

--- GAMS/CPLEX licensed for continuous and discrete problems.
--- GMO Q Extraction (ThreePass): 0.00s
--- GMO setup time: 0.00s
--- GMO memory 0.52 Mb (peak 0.52 Mb)
--- Dictionary memory 0.00 Mb
--- Cplex 22.1.2.0 link memory 0.00 Mb (peak 0.01 Mb)
--- Starting Cplex

Version identifier: 22.1.2.0 | 2024-11-25 | 0edbb82fd
CPXPARAM_Advance                                 0
CPXPARAM_MIP_Display                             4
CPXPARAM_MIP_Pool_Capacity                       0
CPXPARAM_MIP_Tolerances_AbsMIPGap                0
Tried aggregator 1 time.
QCP Presolve eliminated 0 rows and 1 columns.
Aggregator did 1 substitutions.
Reduced QCP has 44 rows, 49 columns, and 88 nonzeros.
Reduced QCP has 17 quadratic constraints.
Presolve time = 0.00 sec. (0.03 ticks)
Parallel mode: using up to 16 threads for barrier.
Number of nonzeros in lower triangle of A*A' = 372
Using Approximate Minimum Degree ordering
Total time for automatic ordering = 0.00 sec. (0.02 ticks)
Summary statistics for Cholesky factor:
  Threads                   = 16
  Rows in Factor            = 44
  Integer space required    = 98
  Total non-zeros in factor = 836
  Total FP ops to factor    = 19870
 Itn      Primal Obj        Dual Obj  Prim Inf Upper Inf  Dual Inf Inf Ratio
   0   3.0771673e+00   1.6629537e+00  3.80e+02  0.00e+00  2.30e+01  1.00e+00
   1   6.8680040e+01   1.1362489e+02  3.80e+02  0.00e+00  2.30e+01  2.13e-02
   2   3.6074399e+01   3.5719296e+01  3.29e+02  0.00e+00  2.00e+01  1.41e+01
   3   1.7305081e+01   1.7491949e+01  6.70e+01  0.00e+00  4.07e+00  3.88e+00
   4   1.6877161e+01   1.6946120e+01  1.12e+01  0.00e+00  6.77e-01  1.10e+01
   5   1.6932308e+01   1.6959923e+01  3.39e+00  0.00e+00  2.06e-01  2.89e+01
   6   1.6888164e+01   1.6893936e+01  1.10e+00  0.00e+00  6.68e-02  1.39e+02
   7   1.6877431e+01   1.6880210e+01  2.27e-01  0.00e+00  1.37e-02  2.94e+02
   8   1.6878549e+01   1.6880381e+01  9.79e-02  0.00e+00  5.94e-03  4.58e+02
   9   1.6875200e+01   1.6875394e+01  5.53e-02  0.00e+00  3.36e-03  4.33e+03
  10   1.6875281e+01   1.6875425e+01  5.75e-03  0.00e+00  3.49e-04  5.88e+03
  11   1.6875121e+01   1.6875137e+01  4.06e-03  0.00e+00  2.46e-04  5.48e+04
  12   1.6875102e+01   1.6875104e+01  4.25e-04  0.00e+00  2.58e-05  4.72e+05
  13   1.6875102e+01   1.6875102e+01  4.93e-05  0.00e+00  2.99e-06  2.21e+06

--- QCP status (1): optimal.
--- Cplex Time: 0.02sec (det. 0.76 ticks)

--- Computing dual multipliers for quadratic constraints...
--- Dual multiplier computation can be disabled using 'calcqcpduals = 0'.

Optimal solution found
Objective:           16.875102

For some other methods for this problem see [2]. They mention the SOCP model but not my proposed NLP approach. 

Next questions:

• Can we extend the principle of a convex hull to circles/disks.
• What about a model with inner ellipses and also an enclosing ellipse 

There are algorithms that extend the concept of a convex hull to circles (or disks) [3,4]. For our data set this convex hull can look like:

----     69 PARAMETER circles  coordinates of center and radius

                   x           y           r        hull    interior

circle1        4.294      21.082       1.663                   1.000
circle2       13.759       7.528       4.014                   1.000
circle3        7.305       5.601       1.961                   1.000
circle4        8.746      21.407       6.235       1.000
circle5        1.678      12.505       2.591                   1.000
circle6       24.953      14.468       2.715       1.000
circle7       24.778      19.056       4.564       1.000
circle8        3.267      15.993       5.336       1.000
circle9        3.988       6.252       4.769       1.000
circle10      16.723      10.884       3.783                   1.000
circle11       8.993       8.786       3.480                   1.000
circle12       3.287       3.753       1.706                   1.000
circle13      14.728      20.772       2.885                   1.000
circle14       5.770      16.643       1.279                   1.000
circle15      19.396       7.591       3.031       1.000

Obviously, this is much better than just dropping the circles that are completely inside another circle.

GAMS code

$onText     Find convex hull of a set of disks.   $offText   *--------------------------------------------------------------------- * data *---------------------------------------------------------------------   sets   i 'inner circles' /circle1*circle15/   a 'attributes for storing results' /x,y,r/   k(a) 'coordinate' /x,y/ ;   alias (i,j);   Parameters  circles(*,*)   'coordinates of center and radius' ;   circles(i,k) = uniform(0,25); circles(i,'r') = uniform(1,7); display circles;     *--------------------------------------------------------------------- * find convex hull *---------------------------------------------------------------------   set hull(i) 'circles forming the convex hull';   * assumes shapely is installed in gmspython   embeddedCode Python:   import numpy as np import gams.transfer as gt import pandas as pd   from shapely.geometry import Point from shapely.ops import unary_union   m = gt.Container(gams.db) df = m.data["circles"].records   cols = df.columns df2 = df.pivot(index=cols[0],columns=cols[1],values=cols[2]) print(df2)   # Build circle geometries geoms = {cid: Point(x, y).buffer(r) for cid, x, y, r in df2.itertuples()}   # Convex hull of union union = unary_union(list(geoms.values())) hull = union.convex_hull   # Find circles that touch the hull boundary hull_circles = [cid for cid, g in geoms.items() if g.intersects(hull.boundary)] print(hull_circles) gams.set("hull", hull_circles)   endEmbeddedCode hull     circles(i,'hull') = hull(i); circles(i,'interior') = not hull(i); display circles;     *--------------------------------------------------------------------- * visualization *---------------------------------------------------------------------   $set html  plot.html $set data  data.js     $if %htmlplot%==0 $goto skipplot   file fdata /%data%/; put fdata;     set    sd /inner,hull/    sdi(sd,*)   ; sdi('inner',i)$(not hull(i)) = yes; sdi('hull',i)$hull(i) = yes;   display sdi;   Scalars    firstred /1/    firstblue /1/ ;   alias (*,item);   * circles put "circles = ["/; loop(sdi(sd,item),   put " {"/;   put "   type:'circle',"/;   put "   xref:'x',"/;   put "   yref:'y',"/;   put "   x0:",(circles(item,'x')-circles(item,'r')):0:4,","/;   put "   y0:",(circles(item,'y')-circles(item,'r')):0:4,","/;   put "   x1:",(circles(item,'x')+circles(item,'r')):0:4,","/;   put "   y1:",(circles(item,'y')+circles(item,'r')):0:4,","/;   if (sameas(sd,'hull'),      put "   line:{color:'Red'},"/;      put "   fillcolor:'Pink',"/;      put$firstred "   name:'hull',"/;      put$firstred "   showlegend:true,"/;      firstred = 0;   else      put "   line:{color:'Blue'},"/;      put "   fillcolor:'LightSkyBlue',"/;      put$firstblue "   name:'interior',"/;      put$firstblue "   showlegend:true,"/;      firstblue = 0;   );   put "   opacity:0.5,"/;    put " },"/; ); put "];"/;   scalars xmin,xmax,ymin,ymax,range0,range1; xmin = smin(i,circles(i,'x')-circles(i,'r')); xmax = smax(i,circles(i,'x')+circles(i,'r')); ymin = smin(i,circles(i,'y')-circles(i,'r')); ymax = smax(i,circles(i,'y')+circles(i,'r')); range0 = floor(min(xmin,ymin)); range1 = ceil(max(xmax,ymax));   put "range0=",range0:0:4,";"/; put "range1=",range1:0:4,";"/;   putclose;   $onecho > %html% <html> <script src="https://cdn.plot.ly/plotly-3.4.0.min.js" charset="utf-8"></script> <script src="%data%" charset="utf-8"></script> <h1>Smallest Encompassing Circle</h1> <div id="plotDiv"></div> <script>   var layout = {   autosize: false,   width: 700,   height: 650,   showlegend: true,   shapes: circles,   xaxis: { range: [range0,range1] },   yaxis: { range: [range0,range1] } } var options = {staticPlot: true, displayModeBar: false}; Plotly.newPlot('plotDiv', [], layout, options); </script> </html> $offecho   executetool 'win32.ShellExecute "%html%"';   $label skipplot

 

References

1. Minimum enclosing ellipse, https://yetanothermathprogrammingconsultant.blogspot.com/2026/04/minimum-enclosing-ellipse.html
2. Xu, S., Freund, R. M., & Sun, J. (2003). Solution methodologies for the smallest enclosing circle problem. Computational Optimization and Applications, 25(1), 283-292.
3. David Rappaport, A convex hull algorithm for discs, and applications, Computational Geometry: Theory and Applications 1 (1992) 171-187.
4. Linh Kieu Nguyen, Chanyoung Song, Joonghyun Ryu, Phan Thanh An, Nam-Dũng Hoang, Deok-Soo Kim, QuickhullDisk: A faster convex hull algorithm for disks, Applied Mathematics and Computation, 363, 2019, 124626.