We all love ODF, the best and the most vendor-neutral file-format in the Universe and its surroundings. But for sure, we have some spots where we would prefer it to be somehow less cumbersome. My favourite spot is the need to compute the bounding box of the path element when one writes the draw:path into an OpenDocument Graphics file. Without proper svg:x, svg:y, svg:height, svg:width and svg:viewBox values the path will not be correctly placed.
Computing bounding boxes is not so complicated work when everything is a polygon (which is the case in the internal model of basegfx module), but it becomes a bit more complicated if an application wants to generate paths including elliptical arcs.
I hit this problem about a year ago when I was working during my hackweek on an improvement of libwpg. I tried first to implement the bounding box of an elliptical arc the lazy hacker way, by googling for what other people did. And to my surprise, there is a huge vacuum in what concerns computation of a bounding box of the "A" element of an SVG path. So, I settled for the lazy hacker's plan B: I abandoned the idea by saying I will handle it later, in the first moment, and by implementing a suboptimal solution in the second time. But, since Eilidh did some spectacular progress last week in handling elliptic arcs, my lazyness became the bottle-neck of the progress. So, it was time to remember those nice times when I was warming the benches of the University, dust off my knowledge of mathematical analysis and analytical algebra (or the lack thereof), and try to compute the bounding box of an elliptical arc myself.
And for the purpose of people that might be as lazy as me, I decided to fight my lazyness the second time to give Uncle Google the opportunity to spit out something meaningful, when someone asks it about "Bounding box of an elliptical arc". Here are the notes:
Compute extremes using parametric description of ellipse
Let us start from this parametric description of an ellipse:
x(theta) = cx + rx*cos(theta)*cos(phi) - ry*sin(theta)*sin(phi)y(theta) = cy + rx*cos(theta)*sin(phi) + ry*sin(theta)*cos(phi)
where cx and cy are the coordinates of the centre of the ellipse, rx and ry are the radii and phi is the x-axis-rotation.
To compute the bounding box of the whole ellipse we need to find for which value of theta the above mentioned functions reach the local extremes. It means where the first derivatives of x and y according to theta are zero. We will get this two equations:
0 = -rx*sin(theta)*cos(phi) - ry*cos(theta)*sin(phi)0 = -rx*sin(theta)*sin(phi) - ry*cos(theta)*cos(phi)
which give us two solutions for x:
theta = -atan(ry*tan(phi)/rx) and theta = M_PI -atan(ry*tan(phi)/rx)
and two solutions for y:
theta = atan(ry/(tan(phi)*rx)) and theta = M_PI + atan(ry/(tan(phi)*rx))
Compute the center of the ellipse
Since we know now the values of theta describing the extremes of our ellipse, we can compute the x and y values of the bounding box of the whole ellipse. Just to do that, we still need to know the coordinates of the center of the ellipse, cx and cy.
The computation of the center of the ellipse is pretty well described in the Appendix F.6.5 of the SVG standard and does not need to be reproduced here. Just note that for this we need the coordinates of the starting point of the arc that correspond to the end point of the previous path segment.
Determine the bounding box of the whole ellipseCompute the xmin, xmax, ymin and ymax using the values of theta for the local extremes and the newly computed cx and cy coordinates. Like this not only we will know the bounding box of the whole ellipse, but we will also know which value of theta corresponds to maximum and which one to minimum. This knowledge will be later valuable for determining the tightest possible bounding box of a given elliptical arc.
Tightest possible bounding box
By calculation of the bounding box of the whole ellipse, we now know the rectangle that will contain the ellipse and thus our elliptical arc too. Nonetheless, this rectangle is too big for our arc. So, the next thing is to trim it down so that it becomes the tightest possible rectangle that will still contain the whole arc.
For this task, we will use the polar coordinates rather then the cartesian ones. The principle is that if any of the points corresponding to xmin, xmax, ymin or ymax of the whole ellipse, lie on the arc they will be be the extremes of the arc too. Nevertheless, if for instance the point ymin does not lie on the arc, the new ymin will be the minimum of the y coordinates of the starting and ending points. In the same way, if the point xmax does not lie on the arc, the new xmax will be the maximum of the x coordinates of the starting and ending points. Whether an extreme does or does not lie on our arc is something trivial to see once the arc is drawn, to determine it programatically will require some efforts.
First, we will compute the coordinates of the points where the whole ellipse touches the bounding box using the parametric description of the ellipse and the values of the theta that we found out in the previous steps. And for determination whether they lie or not on our arc we will use their position in polar coordinates. We will thus need to compute the angles with the x-axis of the lines going through the center of the ellipse and our extreme points. In other terms, we will compute the angle between vector (1,0) and vector (xextreme-cx, yextreme-cy).
The formula for computing the angle between two vectors is known and mentioned inter alia as formula F.6.5.4 of the SVG standard. Generally, the expression to calculate the angle between a vector (ax,ay) and a vector (bx,by) is:
(ax * by > ay * bx ? 1.0 : -1.0) * acos( (ax * bx) + (ay * by) / ( sqrt(ax * ax + ay * ay) * sqrt(bx * bx + by * by) ) ).
But since we already know that the first vector is (1,0), we can simplify it:
(by > 0.0 ? 1.0 : -1.0) * acos( bx / sqrt(bx * bx + by * by) ), which could be eventually simplified to atan(by / bx), but this expression has a potential division by zero and the code would have to handle those border situations in a special way.
Once we know the angles of the extremes, we still need to calculate the angles of the starting and the end points or our arc using exactly the same formula. So we get angle1 corresponding to our starting point and angle2 corresponding to the endpoint. It is necessary to normalize all angles so that they lie in the interval of [0.0, 2.0*M_PI).
In case the sweep flag is 0, the angles are decreasing when the ellipse is drawn. But, for the computation of bounding box the direction of rotation is irrelevant and only complicates the situation. So we swap the angles if the sweep flag is not set. In this way, we can just check for the absence of the extreme points on our elliptical arc, rotating from angle1 to angle2. Nevertheless, we have another difficulty with the fact that the angle of 0 radians is the same as the one of 2*M_PI radians. This passage through the 2*M_PI / 0 border is not very easy to handle directly. That is why we swap the points in case where angle1 > angle2 and will not look in this case for absence of the extreme points on the arc, but for their presence on the complement arc that would close the ellipse.
And as my teachers used to say: "Grey is the theory, but green is the tree of life," here is what it looks like in a plain C++:
#include <algorithm>
#include <cmath>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
inline double getAngle(double bx, double by)
{
return fmod(2*M_PI + (by > 0.0 ? 1.0 : -1.0) * acos( bx / sqrt(bx * bx + by * by) ), 2*M_PI);
}
void EllpArcBBox(double x1, double y1,
double rx, double ry, double phi, bool largeArc, bool sweep, double x2, double y2,
double &xmin, double &ymin, double &xmax, double &ymax)
{
if (rx < 0.0)
rx *= -1.0;
if (ry < 0.0)
ry *= -1.0;
if (rx == 0.0 || ry == 0.0) {
xmin = (x1 < x2 ? x1 : x2);
xmax = (x1 > x2 ? x1 : x2);
ymin = (y1 < y2 ? y1 : y2);
ymax = (y1 > y2 ? y1 : y2);
return;
}
const double x1prime = cos(phi)*(x1 - x2)/2 + sin(phi)*(y1 - y2)/2;
const double y1prime = -sin(phi)*(x1 - x2)/2 + cos(phi)*(y1 - y2)/2;
double radicant = (rx*rx*ry*ry - rx*rx*y1prime*y1prime - ry*ry*x1prime*x1prime);
radicant /= (rx*rx*y1prime*y1prime + ry*ry*x1prime*x1prime);
double cxprime = 0.0;
double cyprime = 0.0;
if (radicant < 0.0) {
double ratio = rx/ry;
double radicant = y1prime*y1prime + x1prime*x1prime/(ratio*ratio);
if (radicant < 0.0) {
xmin = (x1 < x2 ? x1 : x2);
xmax = (x1 > x2 ? x1 : x2);
ymin = (y1 < y2 ? y1 : y2);
ymax = (y1 > y2 ? y1 : y2);
return;
}
ry=sqrt(radicant);
rx=ratio*ry;
} else {
double factor = (largeArc==sweep ? -1.0 : 1.0)*sqrt(radicant);
cxprime = factor*rx*y1prime/ry;
cyprime = -factor*ry*x1prime/rx;
}
double cx = cxprime*cos(phi) - cyprime*sin(phi) + (x1 + x2)/2;
double cy = cxprime*sin(phi) + cyprime*cos(phi) + (y1 + y2)/2;
double txmin, txmax, tymin, tymax;
if (phi == 0 || phi == M_PI) {
xmin = cx - rx;
txmin = getAngle(-rx, 0);
xmax = cx + rx;
txmax = getAngle(rx, 0);
ymin = cy - ry;
tymin = getAngle(0, -ry);
ymax = cy + ry;
tymax = getAngle(0, ry);
} else if (phi == M_PI / 2.0 || phi == 3.0*M_PI/2.0) {
xmin = cx - ry;
txmin = getAngle(-ry, 0);
xmax = cx + ry;
txmax = getAngle(ry, 0);
ymin = cy - rx;
tymin = getAngle(0, -rx);
ymax = cy + rx;
tymax = getAngle(0, rx);
} else {
txmin = -atan(ry*tan(phi)/rx);
txmax = M_PI - atan (ry*tan(phi)/rx);
xmin = cx + rx*cos(txmin)*cos(phi) - ry*sin(txmin)*sin(phi);
xmax = cx + rx*cos(txmax)*cos(phi) - ry*sin(txmax)*sin(phi);
if (xmin > xmax) {
std::swap(xmin,xmax);
std::swap(txmin,txmax);
}
double tmpY = cy + rx*cos(txmin)*sin(phi) + ry*sin(txmin)*cos(phi);
txmin = getAngle(xmin - cx, tmpY - cy);
tmpY = cy + rx*cos(txmax)*sin(phi) + ry*sin(txmax)*cos(phi);
txmax = getAngle(xmax - cx, tmpY - cy);
tymin = atan(ry/(tan(phi)*rx));
tymax = atan(ry/(tan(phi)*rx))+M_PI;
ymin = cy + rx*cos(tymin)*sin(phi) + ry*sin(tymin)*cos(phi);
ymax = cy + rx*cos(tymax)*sin(phi) + ry*sin(tymax)*cos(phi);
if (ymin > ymax) {
std::swap(ymin,ymax);
std::swap(tymin,tymax);
}
double tmpX = cx + rx*cos(tymin)*cos(phi) - ry*sin(tymin)*sin(phi);
tymin = getAngle(tmpX - cx, ymin - cy);
tmpX = cx + rx*cos(tymax)*cos(phi) - ry*sin(tymax)*sin(phi);
tymax = getAngle(tmpX - cx, ymax - cy);
}
double angle1 = getAngle(x1 - cx, y1 - cy);
double angle2 = getAngle(x2 - cx, y2 - cy);
if (!sweep)
std::swap(angle1, angle2);
bool otherArc = false;
if (angle1 > angle2) {
std::swap(angle1, angle2);
otherArc = true;
}
if ((!otherArc && (angle1 > txmin || angle2 < txmin)) || (otherArc && !(angle1 > txmin || angle2 < txmin)))
xmin = x1 < x2 ? x1 : x2;
if ((!otherArc && (angle1 > txmax || angle2 < txmax)) || (otherArc && !(angle1 > txmax || angle2 < txmax)))
xmax = x1 > x2 ? x1 : x2;
if ((!otherArc && (angle1 > tymin || angle2 < tymin)) || (otherArc && !(angle1 > tymin || angle2 < tymin)))
ymin = y1 < y2 ? y1 : y2;
if ((!otherArc && (angle1 > tymax || angle2 < tymax)) || (otherArc && !(angle1 > tymax || angle2 < tymax)))
ymax = y1 > y2 ? y1 : y2;
}
