How to find or estimate the area of intersection of an ellipse with each square of a grid

Question

I am trying to estimate the overflow of water from one container into its adjacent containers for a computer simulation. To do this I have modeled a situation where the location the water will fall in is modeled by an elliptical cylinder and each cell is a 1x1 unit area in a coordinate system.

As such I have the equation ((x-a)/w)^2+((y-b)/h)^2=1 representing the water where (a,b) is the location of the center of the ellipse and w is the width and h is the height of the elliptical cylinder. There are two constraints on these dimensions. a and b will be between 0 and 1 and w and h will be between 1 and 2.

Here is a sort of diagram of what is going on:

I need to find an estimate the percentage of the water that falls in each 1x1 square in the coordinate grid in order to determine how much water has fallen in that particular container in the simulation.

Of course I could simulate this situation and estimate the values for a specific ellipsoid but it would be great to generalize it to an equation.

Answer 1

The best way to do this, I believe, is to use line integrals and compute the areas with Green's Theorem. For each square, we evaluate the integral $$\int_C-ydx,$$ where $C$ is the boundary the intersection of the square with the ellipse, oriented counterclockwise. $C$ is composed of line segments and elliptical arcs. We need to figure out where they start and end, and we need to parameterize the various portions of the path.

Let's say the vertices of the square are $(x,y), (x+1,y), (x+1,y+1), (x,y+1),$ where you'll note that I've listed them in counterclockwise order. The integral is very easy to evaluate over a vertical or horizontal line segment.

Over a horizontal segment from $(x_1,y)$ to $(x_2,y)$ the integral becomes $$ -\int_{x_1}^{x_2}{ydx}= y(x_1-x_2).$$ Notice that when we go from left to right, we get the opposite sign from when we go from right to left.

The integral is even easier to evaluate over a vertical segment. Since $x$ is constant, the integral is $0$. You should check that if the square is entirely inside the ellipse, so that $C$ is the boundary of the square, then these rules reassuringly give $1$ as the area of the square.

Of course, the interesting part is when the boundary of the ellipse intersects the square. (Actually the ellipse is the boundary, but I'm calling the region bounded by the curve "the ellipse," because I don't know a good name for it. The term "ellipsoid" that you use in your question should actually be "elliptical cylinder," as noted by Jean Marie.)

Back to business. We can parameterize the ellipse $$ \left(\frac{x-a}{w}\right)^2+\left(\frac{y-b}{h}\right)^2=1$$ by $$ \begin{align} x &=a+w\cos \theta,\\ y &=b+h\sin \theta,\\ 0&\le\theta\le 2\pi \end{align} $$

When $C$ is the arc of the ellipse from $\theta_1$ to $\theta_2,$ the integral becomes $$ \int_{\theta_1}^{\theta_2}{(b+h\sin\theta)(w\sin\theta)d\theta} =F(\theta_2)-F(\theta_1),\\ \text{where } F(\theta)= -bw\cos\theta + \frac{hw\theta}{2}-\frac{hw\sin 2\theta}{4} $$ You should check that when $\theta_1=0, \theta_2=2\pi,$ so that we are integrating over the full boundary of the ellipse, $F(\theta_2)-F(\theta_1)=hw\pi,$ the area of the ellipse.

Now it's mostly a programming problem. The details are tedious, and I can't detail the full algorithm here, but I'll outline the salient points. Before anything else, let me mention that I would certainly add up all the areas computed and checked that they sum to $hw\pi.$

First thing is to figure out what grid squares might possibly intersect the (interior of the) ellipse. I will designate a square by the coordinates of its lower-left hand (southwest)corner. We want all squares $(x,y)$ where $$ \lfloor a-w \rfloor\le x \le \lfloor a+w \rfloor\,\text{ and } \lfloor b-h \rfloor\le y \le \lfloor b+h \rfloor\ $$ As an aside, check that whatever language you are using computes the floor correctly for negative numbers. Some languages truncate towards $0$ so that $\lfloor -2.4\rfloor$ is computed as $-2$ instead of $-3$. If this is the case in your language of choice, you need to subtract $1.$

Now, for each such square, we call a function to compute the area. I shall describe this function next. We need to know which corners of the square a are inside the ellipse. This is easy to tell. For each vertex $(x,y)$, we compute $\left(\frac{x-a}{w}\right)^2+\left(\frac{y-b}{h}\right)^2.$ The vertex is inside if and only if the value computed is $\le 1.$

The intricate part comes when some of the corners are inside, and some not. I've been thinking about how to do this, and I think it's best to find a corner inside the ellipse, and start from there, working around the boundary counterclockwise. I would call a different routine depending on which was the first corner inside the ellipse that I encountered. So, if $(x,y),$ the lower left-hand corner were inside the ellipse I would call a function southwest. If $(x,y)$ were outside, but $(x+1,y)$ were inside I would call a function southeast and so on. If none of the corners were inside, the area function would return $0.$

area(x,y,a,b,w,h):
    if inside(x,y,a,b,w,h) return southwest(x,y,a,b,w,h)
    if inside(x+1,y,a,b,w,h) return southeast(x,y,a,b,w,h)
    if inside(x+1,y+1,a,b,w,h) return northeast(x,y,a,b,w,h)
    if inside(x,y+1,a,b,w,h) return northwest(x,y,a,b,w,h)
    return 0   

Now, I'll try to describe the southwest function, and leave it up to you to work out the others. First, we deal with the lower edge. If $(x+1,y)$ is inside, the integral is just $-y.$ If it is outside, then we need to find the value $\xi$ so that the lower edge intersects the ellipse in $(\xi,y).$ This is easy. We substitute $y$ in the equation for the ellipse and solve for $x$. We get $2$ values, and we choose the one between $x$ and $x+1.$ Now the integral on the lower edge goes from $x$ to $\xi$ so we get $y(\xi-x).$

At this point, we will have to integrate over an arc of the ellipse. We need to know what is the value of $\theta_1$ that corresponds to $(\xi,y).$ This is straightforward. We have to solve $$ \begin{align} \xi &=a+w\cos \theta\\ y &=b+h\sin \theta\\ \end{align} $$ for $\theta.$ Just remember that the range of the inverse trigonometric functions is limited, so you may need to adjust the result to get $\theta_1$ in the proper quadrant. You can tell what quadrant it should be in by reference to $(a,b).$

Now that we've found $\theta_1,$ where does this elliptical arc end? That depends on which is the next vertex inside the ellipse. If $(x+1,y+1)$ is inside the ellipse, then the ellipse must intersect the line segment from $(x+1,y)$ to $(x+1,y+1)$ in some point $(x+1,\eta).$ We find $\eta$ in a manner exactly analogous to the way we found $\xi,$ and we find $\theta_2$ corresponding to $(x+1,\eta)$ in the same way that we found $\theta_1.$ Now we can compute $F(\theta_2)-F(\theta_1)$ and add it to the integral we computed previously. The line segment from $(x+1,\eta)$ to $(x+1,y+1)$ is vertical, so we can skip it. Now we proceed in the same manner from $(x+1,y+1)$ until we get back to $(x,y).$

The foregoing was the case where $(x+1,y+1)$ is inside. If it is outside, then we need to check $(x,y+1).$ If $(x,y+1)$ is inside, then the elliptical arc intersects the line segment from $(x+1,y+1)$ to $(x,y+1)$ in some point $(\xi',y+1).$ Solve for $\xi'$ and $\theta_2$ as before as compute the integral with $F$. Now the line segment from $(\xi',y+1)$ to $(x,y+1)$ is horizontal, so we get $y(\xi'-x)$ for the integral. In this scenario, both $(x,y+1)$ and $(x,y)$ are inside, so the last part of the path is a vertical segment, which we can ignore.

There remains the case where $(x,y)$ is the only vertex inside. In that case, the ellipse intersects the line segment from $(x,y+1)$ to $(x,y)$ in some point $(x,\eta)$ and we know how to handle this. In this case, the last part of the path is the vertical segment from $(x,\eta)$ to $(x,y),$ which we can ignore.

If $(x,y)$ is outside the circle but $(x+1,y)$ is inside the circle, then we would call the southeast function. This would be analogous to the southwest function but would start at $(x+1,y)$. A small efficiency gain is possible in that we don't have to check whether $(x,y)$ is inside; we know it isn't.

I hope this outline is enough to allow you to write the code. I'm hampered by my inability to draw pictures.

EDIT Rats. This isn't quite right. I just realized that it's possible for a square to intersect the ellipse even though none of its corners are inside. Imagine a square that just clips off a bit of the ellipse at the top. The general approach is still correct. We just need a way to detect the situation and another set of routines top, bottom, right, left to deal with the possibilities.

I don't really think this proliferation of routines is a problem, since there will be much code reuse possible. I think it will probably be easier to write the program than to describe it in words.

EDIT I think it will be best to assume we call four routines top, bottom, left, and right to compute the missing areas, if any. There will be a minuscule loss of efficiency, but a great gain in clarity, especially in the context of an answer on MSE. I'll describe top in some detail, and just sketch the others.

We first take the largest value $\eta$ that could be the $y-$coordinate of the lower edge of a square. Then we find the points where the line $y=\eta$ intersects the ellipse. Say the intersection points are $(x_1,\eta), (x_2, \eta).$ For each $x$ that could be the $x-$coordinate of the left-hand side of a square, we check whether $x\le x_1$ and $x_2\le x+1.$ If so, we have a square that intersects the ellipse that has no corners inside, and we compute the area. There can't be more than one such square, so we don't need to check the others. The area we want to compute is bounded by the line segment from $(x_1,\eta) \text{ to }(x_2, \eta)$ and by the elliptical arc from $(x_2,\eta) \text{ to }(x_1, \eta).$ Note the difference in orientation, in order to orient the path counterclockwise.

top(a,b,w,h):
    eta = floor(b+h)
    delta = w*sqrt(1 - ((eta-b)/h)**2)
    x1 = a-delta
    x2 = a+delta
    for x = floor(a-w) to floor(a+w)
        if x <= x1 and x2 <= x+1:
            return topArea(x1, x2, eta)
    return 0

The function bottom differs in two respects. First, we want $\eta$ to be the smallest value that could be the top edge of a square, so we have $\eta = 1+\lfloor b-h\rfloor.$ Second, to orient the path counterclockwise, we need then line segment from $(x_2,\eta) \text{ to }(x_1, \eta)$ and the elliptical arc from $(x_1,\eta) \text{ to }(x_2, \eta),$ exactly the reverse of the top. Alternatively, you can have just one function, topBottomArea and have it return the absolute value. When you use the clockwise orientation of the path, you get the negative of the area.

The functions left and right are similar, but you fix $\xi$ and solve for $y_1, y^2.$ On the right, we want $\xi = \lfloor a+w\rfloor$ and on the left, we want $\xi = 1+\lfloor a-w\rfloor$. We still have to worry about the orientation of the elliptical path (or take absolute values) by we can ignore the vertical line segment, as always.

I think this is a complete algorithm now. I have not considered roundoff error. It's not likely to matter much because there re very few squares involved. With the conditions you gave, there are a maximum of $25$ squares, and some will be wholly inside or wholly outside the ellipse. There's always the possibility of wrongly deciding whether a point is inside or outside the ellipse, because of the imprecision of floating point calculation, but again that isn't likely to matter much. The point will have to be very close to the boundary for that that to happen, so that the area we compute will be very close to the true area.

I can see one difficulty. If the point is actually outside the ellipse, and we wrongly decide it is inside, we may look for the intersection of the ellipse with the next side, when no such intersection actually exists. Probably the inside function should return True if the point is very close to the ellipse. Then we'll use a line segment, and there's no chance of trying to compute the square root of a negative number.