如何计算椭圆的轴对齐边界框?

38

如果椭圆的主轴是垂直或水平的,那么计算其边界框就很容易,但如果椭圆被旋转了呢?

目前我能想到的唯一方法就是计算周长上的所有点,并找到最大/最小 x 和 y 值。看起来应该有更简单的方法。

如果存在一个描述任意角度下椭圆的函数(数学意义上的),那么我可以使用它的导数来找到斜率为零或未定义的点,但我似乎找不到这样的函数。

编辑:为了澄清,我需要轴对齐的边界框,即它不会随着椭圆旋转而旋转,而是与 x 轴保持对齐,因此转换边界框是行不通的。

13个回答

45

您可以尝试使用参数方程式来表示任意角度旋转的椭圆:

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,并绕角度phi旋转。

然后可以求导并解出梯度=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) 处、长轴为 2,短轴为 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
and
x = 2*cos(-3.6052)*cos(PI/4) - sin(-3.6052)*sin(PI/4) = -1.5811

@sh1ftst0rm 你是对的。这是个奇怪的错误,因为我记录了正确的结果。按照所写的计算会得出-0.9487,那是相当错误的。我会纠正它。 - Mike Tunnicliffe
1
这个答案有一件事情我不太明白:为什么我们只关心 t = -0.4636 和 t = -3.6052 ?我知道这些结果是通过将变量 'n' 分别替换为 0 和 -1,带入方程式 t = -0.4636 + n*PI 中得出的,但不知道这些特定值的重要性。是因为插入原始方程中会产生重复的笛卡尔坐标吗? - Caleb
1
@Caleb: 完全正确。在单个完整的2PI旋转中,只有两个唯一的解决方案适用(因为我们有...+nPI条款)。通过选择任何2个连续的t值,我们确保从将这些值插入方程式中得到两个唯一的答案。您可以尝试其他t值来确认这一点。 - Mike Tunnicliffe
1
显然,你的解决方案对于phi=90°和phi=270°无效。 - Maksym Ganenko
5
这是Python的实现,可以在此处找到:https://gist.github.com/smidm/b398312a13f60c24449a2c7533877dc0。 - Matěj Šmíd
显示剩余5条评论

12

我在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);

JavaScript 实现 https://gist.github.com/kolosovsky/b56d90ad3a507876ae2da96e5bb8ff71 - Eduard Kolosovskyi

7
这很简单,但有点难以解释,因为您没有告诉我们如何表示您的椭圆。 有很多方法可以做到这一点.. 无论如何,总的原则是:您无法直接计算轴对齐边界框。 但是,您可以将椭圆在x和y中的极值作为2D空间中的点来计算。 对此,只需使用方程x(t)= ellipse_equation(t)和y(t)= ellipse_equation(t)。 获取它的一阶导数并将其解决为其根。 由于我们处理的是基于三角函数的椭圆,因此这是直截了当的。 您应该得到一个方程,可以通过atan,acos或asin获得根。

提示:要检查您的代码,请尝试使用未旋转的椭圆:您应该在0,Pi / 2,Pi和3 * Pi / 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;

1
你的代码第三行有一个错别字。第二个“=”应该是“*”。我认为这对每个人来说都很明显。我无法编辑它,因为编辑必须至少有6个字符,而这只有一个字符需要编辑。 - miile7

4

Brilian Johan Nilsson。 我已经将您的代码转换为C# - 椭圆角度现在以度为单位:

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);
}

我相信应该是 ellipseCenterY - height * 0.5(减号,而不是加号)。否则完美无缺。 - JeremyDouglass

2
我认为最有用的公式是这个。从原点旋转角度phi的省略号具有以下方程式:
(x-h)^2/a^2 + (y-k)^2/b^2 = 1
其中(h,k)是中心,a和b是长轴和短轴的大小,t的变化范围为-pi到pi。
通过这个公式,你应该能够推导出dx/dt或dy/dt何时为0。

我现在感觉很慢,花了很长时间才写完我的答案 T.T - Mike Tunnicliffe

1
如果您使用OpenCV/C++并使用cv::fitEllipse(..)函数,则可能需要椭圆的边界矩形。这里我使用Mike的答案提供了一个解决方案:
// 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;
}

1
这是一个基于上面答案的TypeScript函数。
export function getRotatedEllipseBounds(
  x: number,
  y: number,
  rx: number,
  ry: number,
  rotation: number
) {
  const c = Math.cos(rotation)
  const s = Math.sin(rotation)
  const w = Math.hypot(rx * c, ry * s)
  const h = Math.hypot(rx * s, ry * c)

  return {
    minX: x + rx - w,
    minY: y + ry - h,
    maxX: x + rx + w,
    maxY: y + ry + h,
    width: w * 2,
    height: h * 2,
  }
}

1
以下是椭圆的公式,如果通过其焦点和离心率给出(如果通过轴长度、中心和角度给出,请参见用户1789690的答案)。
即,如果焦点为(x0,y0)和(x1,y1),离心率为e,则
bbox_halfwidth  = sqrt(k2*dx2 + (k2-1)*dy2)/2
bbox_halfheight = sqrt((k2-1)*dx2 + k2*dy2)/2

where

dx = x1-x0
dy = y1-y0
dx2 = dx*dx
dy2 = dy*dy
k2 = 1.0/(e*e)

我从用户1789690和Johan Nilsson的答案中推导出了这些公式。

0
这是我用于找到任意方向椭圆的紧密拟合矩形的函数。
我使用OpenCV中的矩形和点进行实现:
cg - 椭圆的中心
size - 椭圆的长轴和短轴
angle - 椭圆的方向
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;
}

0

这是Pranay Soni代码的另一个版本,使用js实现 codepen 我希望有人会发现它有用

/**
* @param {Number} rotation
* @param {Number} majorAxis
* @param {Nmber} minorAxis
* @pivot {Point} pivot {x: number, y: number}
* @returns {Object}
*/
export function getElipseBoundingLines(ratation, majorAxis, minorAxis, pivot) {
   const {cos, sin, tan, atan, round, min, max, PI} = Math;

   let phi = rotation / 180 * PI;

   if(phi === 0) phi = 0.00001;
   // major axis
   let a = majorAxis;
   //minor axis
   let b = minorAxis;

   const getX = (pivot, phi, t)  => {
      return round(pivot.x + a * cos(t) * cos(phi) - b * sin(t) * sin(phi))
   }
   const getY = (pivot, phi, t)  => {
      return round(pivot.y + b * sin(t) * cos(phi) + a * cos(t) * sin(phi))
   }

   const X = [], Y = [];

   let t = atan(-b * tan(phi) / a);
   X.push(getX(pivot, phi, t));
   Y.push(getY(pivot, phi, t));

   t = atan(b * (1 / tan(phi) / a));
   X.push(getX(pivot, phi, t));
   Y.push(getY(pivot, phi, t));

   phi += PI;

   t = atan(-b * tan(phi) / a);
   X.push(getX(pivot, phi, t));
   Y.push(getY(pivot, phi, t));

   t = atan(b * (1 / tan(phi)) / a);
   X.push(getX(pivot, phi, t));
   Y.push(getY(pivot, phi, t));

   const left    = min(...X);
   const right   = max(...X);
   const top     = min(...Y);
   const bottom  = max(...Y);

   return {left, top, right, bottom};
}

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接