математика — Найти эллипс минимальной площади, содержащий множество точек в переполнении стека

Они не идут так далеко, как дать вам код на C ++, но включают в себя подробное обсуждение эффективных алгоритмов того, что вам нужно сделать.

https://www.cs.cornell.edu/cv/OtherPdf/Ellipse.pdf

http://www.stsci.edu/~RAB/Backup%20Oct%2022%202011/f_3_CalculationForWFIRSTML/Gaertner%20&% 20Schoenherr.pdf

как предполагает Рори Донтон, вам нужно четко указать ограничения решения, и устранение любой воли значительно усложнит ситуацию. Для начала предположим это сейчас:

• это 2D проблема
• эллипс выровнен по оси
• центр произвольный вместо (0,0)

Я бы атаковал это как стандарт жанр и тест проблема с поиск приближения (что является гибридом между бинарным и линейным поиском), чтобы ускорить его (но вы также можете попробовать грубую силу с самого начала, чтобы увидеть, работает ли она).

1. вычислить ограничения решения

   Чтобы ограничить поиск, вам нужно найти приблизительное положение размещения и размер эллипса. Для этого вы можете использовать закрашенный круг для ваших очков. Понятно, что площадь эллипса будет меньше или равна кругу, а расположение будет рядом. Круг не обязательно должен быть наименьшим из возможных, поэтому мы можем использовать, например, это:

   1. найти ограничивающий прямоугольник
   2. пусть окружность будет центрирована в этой ограничительной рамке, а радиус будет максимальным расстоянием от его центра до любой из точек.

   Это будет O(n) сложность где n это количество ваших очков.

2. искать «все» возможные эллипсы и помнить лучшее решение

   поэтому нам нужно найти центр эллипса (x0,y0) и полуоси rx,ry в то время как area = M_PI*rx*ry минимален При приближенном поиске каждая переменная имеет коэффициент O(log(m)) и каждая итерация должна проверить правильность, которая O(n) итоговая сложность будет O(n.log^4(m)) где m среднее количество возможных вариаций каждого параметра поиска (зависит от точности и ограничений поиска). С простым грубым поиском это было бы O(n.m^4) что действительно страшно, особенно для плавающей запятой, где m может быть действительно большим.

   Чтобы ускорить это, мы знаем, что площадь эллипса будет меньше или равна площади найденного круга, поэтому мы можем игнорировать все большие эллипсы. Ограничения к rx,ry может быть получено из соотношения сторон ограничительной рамки +/- некоторый резерв.

Здесь простой маленький C ++ пример использования этого approx класс по ссылке выше:

//---------------------------------------------------------------------------
// input points
const int n=15; // number of random points to test
float pnt[n][2];
// debug bounding box
float box_x0,box_y0,box_x1,box_y1;
// debug outscribed circle
float circle_x,circle_y,circle_r;
// solution ellipse
float ellipse_x,ellipse_y,ellipse_rx,ellipse_ry;
//---------------------------------------------------------------------------
void compute(float x0,float y0,float x1,float y1) // cal with bounding box where you want your points will be generated
{
int i;
float x,y;
// generate n random 2D points inside defined area
Randomize();
for (i=0;i<n;i++)
{
pnt[i][0]=x0+(x1-x0)*Random();
pnt[i][1]=y0+(y1-y0)*Random();
}
// compute bounding box
x0=pnt[0][0]; x1=x0;
y0=pnt[0][1]; y1=y0;
for (i=0;i<n;i++)
{
x=pnt[i][0]; if (x0>x) x0=x; if (x1<x) x1=x;
y=pnt[i][1]; if (y0>y) y0=y; if (y1<y) y1=y;
}
box_x0=x0; box_x1=x1;
box_y0=y0; box_y1=y1;
// "outscribed" circle
circle_x=0.5*(x0+x1);
circle_y=0.5*(y0+y1);
circle_r=0.0;
for (i=0;i<n;i++)
{
x=pnt[i][0]-circle_x; x*=x;
y=pnt[i][1]-circle_y; y*=y; x+=y;
if (circle_r<x) circle_r=x;
}
circle_r=sqrt(circle_r);
// smallest area ellipse

int N;
double m,e,step,area;
approx ax,ay,aa,ab;
N=3; // number of recursions each one improves accuracy with factor 10
area=circle_r*circle_r; // solution will not be bigger that this
step=((x1-x0)+(y1-y0))*0.05; // initial position/size step for the search as 1/10 of avg bounding box size
for (ax.init(         x0,          x1,step,N,&e);!ax.done;ax.step())    // search x0
for (ay.init(         y0,          y1,step,N,&e);!ay.done;ay.step())    // search y0
for (aa.init(0.5*(x1-x0),2.0*circle_r,step,N,&e);!aa.done;aa.step())    // search rx
for (ab.init(0.5*(y1-y0),2.0*circle_r,step,N,&e);!ab.done;ab.step())    // search ry
{
e=aa.a*ab.a;
// is ellipse outscribed?
if (aa.a>=ab.a)
{
m=aa.a/ab.a; // convert to circle of radius rx
for (i=0;i<n;i++)
{
x=(pnt[i][0]-ax.a);   x*=x;
y=(pnt[i][1]-ay.a)*m; y*=y;
// throw away this ellipse if not
if (x+y>aa.a*aa.a) { e=2.0*area; break; }
}
}
else{
m=ab.a/aa.a; // convert to circle of radius ry
for (i=0;i<n;i++)
{
x=(pnt[i][0]-ax.a)*m; x*=x;
y=(pnt[i][1]-ay.a);   y*=y;
// throw away this ellipse if not
if (x+y>ab.a*ab.a) { e=2.0*area; break; }
}
}
}

ellipse_x =ax.aa;
ellipse_y =ay.aa;
ellipse_rx=aa.aa;
ellipse_ry=ab.aa;
}
//---------------------------------------------------------------------------

Даже этот простой пример с 15 точками занял около 2 секунд. Вы можете улучшить производительность, добавив эвристику, например, тестировать только области ниже, чем circle_r^2 и т.д., или лучше выберите область решения с некоторым математическим правилом. Если вы используете метод грубой силы вместо поиска с приближением, который ожидает, что время вычислений может составлять даже минуты или более, следовательно, O(scary)…

Осторожно, этот пример не будет работать для любого соотношения сторон точек, так как я жестко закодировал верхнюю границу для rx,ry в 2.0*circle_r что может быть недостаточно. Вместо этого вы можете вычислить верхнюю границу из соотношения сторон точек или условия rx*ry<=circle_r^2…

Существуют также другие («более быстрые») методы, например, вариация CCD (циклическая координата по убыванию). Но такие методы обычно не могут гарантировать, что будет найдено оптимальное решение или вообще нет …

Вот краткий обзор примера вывода:

Точки являются отдельными точками от pnt[n]Серые пунктирные штемпели ограничивают коробку и используют разрисованный круг. Зеленый эллипс найден решением.