What is a clever bounding box for a 3D monte carlo integration?

Question

For a project, I need to calculate the overlap volume of two overlapping ellipsoids in 3D. The method itself isn't a problem: I'm basically picking random points within a bounding box and checking whether they're in both ellipsoids simultaneously.

In my never-ending quest to optimize the program in terms of runtime, a smaller bounding box would obviously be advantageous. Right now, the "box" is a sphere which is centered around the midpoint between the ellipsoids' centers of mass and has a diameter corresponding to the longest ellipsoid axis. This is completely arbitrary and I'm fairly certain the overlapping volume will always be contained in this sphere, but I'd really like to find some way to optimize the entire process.

Is there some general method to optimize the bounding volume?

Answer 1

Isn't this a math problem, rather than a CS problem? What you seem to be asking for is a formula defining the possible overlap of two ellipsoids. You intend to use that formula to make your code more efficient, sure, but that is not related to the core question as far as I can tell. If you had the formula and were trying to figure out how to write it as code, that would be another story. Maybe you should consider posting this to https://math.stackexchange.com/

It seems to me you can redefine your question as "What is formula defining the volume where overlap is possible between 2 ellipsoids in a 3D space?", and that has no reference to computing whatsoever.

Answer 2

This is maybe more of a math than programming question but I'm thinking the intersection of 2 ellipsoids is a convex shape so maybe you could try another approach than Monte Carlo to deal with such a relatively nicely behaved object.

First you need to (somehow) find a point inside both ellipsoids.

Now create a cube centered at that point, large enough to completely contain both ellipsoids. Subdivide the cube along all 3 axes into 8 equally sized sub-cubes (first step of creating an octree).

After this, apply recursion to the subcubes. If some cube corners are inside both ellipsoids and some are not, subdivide the cube again. Also if some areas of cube faces are inside both ellipsoids and others not, subdivide again.

Repeat to desired recursion depth. You can track the calculation error by summing the volume of cubes that were subdivided.

Now the challenging operation here is resolving "if some areas of a cube face are inside both ellipsoids". An intersection between an ellipsoid and a plane is an ellipse. To inspect each of the 6 cube faces, you need to intersect them with the ellipsoids and then intersect the remaining ellipses with each other if they both overlap the cube face. Luckily for 2D ellipse-ellipse overlap test you can find references more easily...

Monte Carlo will give you n bits of accuracy with O(2^n) effort so for example 24 bits takes 16.7 million samples. If you focus on the shape surface with such a divide and conquer method, you should only need O(2^(n*2/3)) samples which is about 65 thousand for 24 bits.

Also your concern about initial bounding box size becomes irrelevant since excess is quickly removed in initial recursion steps.

Answer 3

Here's another idea if you want to stick to Monte Carlo:

Project your ellipsoids on the XZ plane to get 2 ellipses. Calculate the intersection of their bounding rectangles to use as a bounding rectangle for your calculations.

Then shoot rays through the XZ plane within your bounding rectangle, along the Y axis. For each ray, calculate intersections with the ellipsoids and calculate the length of the ray that is inside both. Add this to a variable "hits_length".

The volume of the ellipsoid intersection should be approximated by:

hits_length / rays_shot * bounding_rectangle_area

Again this way you're calculating the volume by sampling a 2D space which should make it converge faster.