Как вычислить ограничивающую рамку эллипса, выровненную по оси?
https://stackoverflow.com/questions/87734
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianВопрос
Если большая ось эллипса вертикальна или горизонтальна, вычислить ограничивающую рамку несложно, но что будет, если эллипс повернуть?
Единственный способ, который я могу придумать на данный момент, - это вычислить все точки по периметру и найти максимальные/минимальные значения x и y.Кажется, должен быть более простой способ.
Если существует функция (в математическом смысле), описывающая эллипс под произвольным углом, то я мог бы использовать ее производную, чтобы найти точки, где наклон равен нулю или не определен, но я не могу ее найти.
Редактировать:чтобы уточнить, мне нужна ограничивающая рамка, выровненная по оси, т.е.он не должен вращаться вместе с эллипсом, а должен оставаться выровненным по оси X, поэтому преобразование ограничивающей рамки не будет работать.
Решение
Вы можете попробовать использовать параметризованные уравнения для эллипса, повернутого на произвольный угол:
x = h + a*cos(t)*cos(phi) - b*sin(t)*sin(phi) [1]
y = k + b*sin(t)*cos(phi) + a*cos(t)*sin(phi) [2]
...где эллипс имеет центр (h,k) большую полуось a и малую полуось b и повернут на угол фи.
Затем вы можете дифференцировать и найти градиент = 0:
0 = dx/dt = -a*sin(t)*cos(phi) - b*cos(t)*sin(phi)
=>
tan(t) = -b*tan(phi)/a [3]
Это должно дать вам множество решений для t (два из которых вас интересуют), подключите их обратно к [1], чтобы получить максимальное и минимальное значения x.
Повторите для [2]:
0 = dy/dt = b*cos(t)*cos(phi) - a*sin(t)*sin(phi)
=>
tan(t) = b*cot(phi)/a [4]
Давайте попробуем пример:
Рассмотрим эллипс в точке (0,0) с a=2, b=1, повернутый на PI/4:
[1] =>
x = 2*cos(t)*cos(PI/4) - sin(t)*sin(PI/4)
[3] =>
tan(t) = -tan(PI/4)/2 = -1/2
=>
t = -0.4636 + n*PI
Нас интересуют t = -0,4636 и t = -3,6052.
Итак, мы получаем:
x = 2*cos(-0.4636)*cos(PI/4) - sin(-0.4636)*sin(PI/4) = 1.5811
и
x = 2*cos(-3.6052)*cos(PI/4) - sin(-3.6052)*sin(PI/4) = -1.5811
Другие советы
Я нашел простую формулу на http://www.iquilezles.org/www/articles/ellipses/ellipses.htm (и игнорировал ось z).
Я реализовал это примерно так:
num ux = ellipse.r1 * cos(ellipse.phi);
num uy = ellipse.r1 * sin(ellipse.phi);
num vx = ellipse.r2 * cos(ellipse.phi+PI/2);
num vy = ellipse.r2 * sin(ellipse.phi+PI/2);
num bbox_halfwidth = sqrt(ux*ux + vx*vx);
num bbox_halfheight = sqrt(uy*uy + vy*vy);
Point bbox_ul_corner = new Point(ellipse.center.x - bbox_halfwidth,
ellipse.center.y - bbox_halfheight);
Point bbox_br_corner = new Point(ellipse.center.x + bbox_halfwidth,
ellipse.center.y + bbox_halfheight);
Это относительно просто, но немного сложно объяснить, поскольку вы не рассказали нам, как вы представляете свой эллипс.Есть так много способов сделать это..
В любом случае, общий принцип звучит так:Вы не можете вычислить граничную рамку, выровненную по оси, напрямую.Однако вы можете вычислить экстремумы эллипса по x и y как точки в 2D-пространстве.
Для этого достаточно взять уравнения x(t) = ellipse_equation(t) и y(t) = ellipse_equation(t).Получите его производную первого порядка и найдите ее корень.Поскольку мы имеем дело с эллипсами, основанными на тригонометрии, это просто.В итоге у вас должно получиться уравнение, корни которого либо получаются через atan, acos или asin.
Намекать:Чтобы проверить свой код, попробуйте использовать неповернутый эллипс:Вы должны получить корни в точках 0, Пи/2, Пи и 3*Пи/2.
Сделайте это для каждой оси (x и y).Вы получите не более четырех корней (меньше, если ваш эллипс вырожден, напримеродин из радиусов равен нулю).Оцените положения корней и получите все крайние точки эллипса.
Теперь вы почти у цели.Получить граничную рамку эллипса так же просто, как сканировать эти четыре точки на предмет xmin, xmax, ymin и ymax.
Кстати, если у вас возникли проблемы с поиском уравнения вашего эллипса:попытайтесь свести это к случаю, когда у вас есть эллипс, выровненный по оси, с центром, двумя радиусами и углом поворота вокруг центра.
Если вы это сделаете, уравнения примут вид:
// the ellipse unrotated:
temp_x (t) = radius.x * cos(t);
temp_y (t) = radius.y = sin(t);
// the ellipse with rotation applied:
x(t) = temp_x(t) * cos(angle) - temp_y(t) * sin(angle) + center.x;
y(t) = temp_x(t) * sin(angle) + temp_y(t) * cos(angle) + center.y;
Брилиан Йохан Нильссон.Я расшифровал ваш код на C# - ellipseAngle теперь в градусах:
private static RectangleF EllipseBoundingBox(int ellipseCenterX, int ellipseCenterY, int ellipseRadiusX, int ellipseRadiusY, double ellipseAngle)
{
double angle = ellipseAngle * Math.PI / 180;
double a = ellipseRadiusX * Math.Cos(angle);
double b = ellipseRadiusY * Math.Sin(angle);
double c = ellipseRadiusX * Math.Sin(angle);
double d = ellipseRadiusY * Math.Cos(angle);
double width = Math.Sqrt(Math.Pow(a, 2) + Math.Pow(b, 2)) * 2;
double height = Math.Sqrt(Math.Pow(c, 2) + Math.Pow(d, 2)) * 2;
var x= ellipseCenterX - width * 0.5;
var y= ellipseCenterY + height * 0.5;
return new Rectangle((int)x, (int)y, (int)width, (int)height);
}
Я думаю, что самая полезная формула — вот эта.Эллипсис, повернутый на угол фи от начала координат, имеет уравнение:
где (h,k) — центр, a и b — размер большой и малой оси, а t варьируется от -pi до pi.
Исходя из этого, вы сможете определить, для какого значения t dx/dt или dy/dt становится равным 0.
Вот формула для случая, если эллипс задан своим очаги и эксцентриситет (для случая, когда оно задано длинами осей, центром и углом, см.г.ответ пользователя 1789690).
А именно, если фокусы равны (x0, y0) и (x1, y1), а эксцентриситет равен e, то
bbox_halfwidth = sqrt(k2*dx2 + (k2-1)*dy2)/2
bbox_halfheight = sqrt((k2-1)*dx2 + k2*dy2)/2
где
dx = x1-x0
dy = y1-y0
dx2 = dx*dx
dy2 = dy*dy
k2 = 1.0/(e*e)
Я вывел формулы из ответа пользователя 1789690 и Йохана Нильссона.
Если вы работаете с OpenCV/C++ и используете cv::fitEllipse(..) вам может понадобиться ограничивающий прямоугольник эллипса.Здесь я принял решение, используя ответ Майка:
// tau = 2 * pi, see tau manifest
const double TAU = 2 * std::acos(-1);
cv::Rect calcEllipseBoundingBox(const cv::RotatedRect &anEllipse)
{
if (std::fmod(std::abs(anEllipse.angle), 90.0) <= 0.01) {
return anEllipse.boundingRect();
}
double phi = anEllipse.angle * TAU / 360;
double major = anEllipse.size.width / 2.0;
double minor = anEllipse.size.height / 2.0;
if (minor > major) {
std::swap(minor, major);
phi += TAU / 4;
}
double cosPhi = std::cos(phi), sinPhi = std::sin(phi);
double tanPhi = sinPhi / cosPhi;
double tx = std::atan(-minor * tanPhi / major);
cv::Vec2d eqx{ major * cosPhi, - minor * sinPhi };
double x1 = eqx.dot({ std::cos(tx), std::sin(tx) });
double x2 = eqx.dot({ std::cos(tx + TAU / 2), std::sin(tx + TAU / 2) });
double ty = std::atan(minor / (major * tanPhi));
cv::Vec2d eqy{ major * sinPhi, minor * cosPhi };
double y1 = eqy.dot({ std::cos(ty), std::sin(ty) });
double y2 = eqy.dot({ std::cos(ty + TAU / 2), std::sin(ty + TAU / 2) });
cv::Rect_<float> bb{
cv::Point2f(std::min(x1, x2), std::min(y1, y2)),
cv::Point2f(std::max(x1, x2), std::max(y1, y2))
};
return bb + anEllipse.center;
}
Этот код основан на коде, представленном пользователем user1789690 выше, но реализован в Delphi.Я проверил это и, насколько я могу судить, работает отлично.Я потратил целый день на поиск алгоритма или какого-то кода, протестировал некоторые из них, которые не работали, и был очень рад, наконец, найти приведенный выше код.Я надеюсь, что кто-то найдет это полезным.Этот код рассчитает ограничивающую рамку повернутого эллипса.Ограничивающая рамка выровнена по оси и НЕ вращается вместе с эллипсом.Радиусы указаны для эллипса до его вращения.
type
TSingleRect = record
X: Single;
Y: Single;
Width: Single;
Height: Single;
end;
function GetBoundingBoxForRotatedEllipse(EllipseCenterX, EllipseCenterY, EllipseRadiusX, EllipseRadiusY, EllipseAngle: Single): TSingleRect;
var
a: Single;
b: Single;
c: Single;
d: Single;
begin
a := EllipseRadiusX * Cos(EllipseAngle);
b := EllipseRadiusY * Sin(EllipseAngle);
c := EllipseRadiusX * Sin(EllipseAngle);
d := EllipseRadiusY * Cos(EllipseAngle);
Result.Width := Hypot(a, b) * 2;
Result.Height := Hypot(c, d) * 2;
Result.X := EllipseCenterX - Result.Width * 0.5;
Result.Y := EllipseCenterY - Result.Height * 0.5;
end;
Это моя функция для поиска плотно прилегающего прямоугольника к эллипсу с произвольной ориентацией.
У меня есть opencv rect и точка для реализации:
cg - центр эллипса
размер - большая, малая оси эллипса
угол - ориентация эллипса
cv::Rect ellipse_bounding_box(const cv::Point2f &cg, const cv::Size2f &size, const float angle) {
float a = size.width / 2;
float b = size.height / 2;
cv::Point pts[4];
float phi = angle * (CV_PI / 180);
float tan_angle = tan(phi);
float t = atan((-b*tan_angle) / a);
float x = cg.x + a*cos(t)*cos(phi) - b*sin(t)*sin(phi);
float y = cg.y + b*sin(t)*cos(phi) + a*cos(t)*sin(phi);
pts[0] = cv::Point(cvRound(x), cvRound(y));
t = atan((b*(1 / tan(phi))) / a);
x = cg.x + a*cos(t)*cos(phi) - b*sin(t)*sin(phi);
y = cg.y + b*sin(t)*cos(phi) + a*cos(t)*sin(phi);
pts[1] = cv::Point(cvRound(x), cvRound(y));
phi += CV_PI;
tan_angle = tan(phi);
t = atan((-b*tan_angle) / a);
x = cg.x + a*cos(t)*cos(phi) - b*sin(t)*sin(phi);
y = cg.y + b*sin(t)*cos(phi) + a*cos(t)*sin(phi);
pts[2] = cv::Point(cvRound(x), cvRound(y));
t = atan((b*(1 / tan(phi))) / a);
x = cg.x + a*cos(t)*cos(phi) - b*sin(t)*sin(phi);
y = cg.y + b*sin(t)*cos(phi) + a*cos(t)*sin(phi);
pts[3] = cv::Point(cvRound(x), cvRound(y));
long left = 0xfffffff, top = 0xfffffff, right = 0, bottom = 0;
for (int i = 0; i < 4; i++) {
left = left < pts[i].x ? left : pts[i].x;
top = top < pts[i].y ? top : pts[i].y;
right = right > pts[i].x ? right : pts[i].x;
bottom = bottom > pts[i].y ? bottom : pts[i].y;
}
cv::Rect fit_rect(left, top, (right - left) + 1, (bottom - top) + 1);
return fit_rect;
}
