# Computing the intersection area of two disks in 3D

Closely related to this question about highlighting intersection of two disks, I am trying to figure out if one can do so similarly for disks embedded in $$3D$$ (e.g. in a bounding box). The difference is that, in $$3D$$ the orientation of the disks matters in how much of overlap/orthogonal-projection there is between them. The orientation of a disk is simply the vector normal to its surface and centered at its center. Therefore, each disk has a center vector (for its position) $$\mathbf v$$ and a normal vector $$\mathbf n$$ for its orientation. As an example, 2 disks $$i,j$$ have maximal overlap if $$\mathbf n_i \parallel \mathbf n_j$$ and the difference vector of their center positions $$\mathbf v_j-\mathbf v_i$$ also being parallel to their normal, then the overlap area is exactly $$\pi r^2,$$ $$r$$ being the radius of the disks.

Intuitively, computing such projection is as if we computed the shadow two drawn particles (here disks) create onto one another when visualizing them.

• But is there a way we can quantify the overlap area between two $$3D$$-embedded disks in Mathematica? Can RegionIntersection be made use of for such application?

To clarify how the overlap between the disks is defined or at least what I mean by it, the idea is to compute the orthogonal projection of their respective surfaces onto one another. For instance given $$2$$ disks $$i,j$$ with their position and normal vectors $$\mathbf v_i,\mathbf n_i$$ and $$\mathbf v_j,\mathbf n_j$$, we can take the average of the orthogonal-surface-projection of disk $$i$$ onto plane of disk $$j$$ with that of disk $$j$$ onto plane of disk $$i$$ which yields a symmetrized definition of overlap or intersection between the disks, taking into account not only their orientations but also relative positions.

Stealing from J. M.'s answer here (its first part), here's an image of one such disk within its plane and its orientation vector visualised (the normal to the plane centered at center of disk):

An attempt to visualize DaveH's suggestion which was very briefly put in their answer:

Say we have one disk centered at v1 and with normal vector n1 and another with v2,n2 as given by (both with diameter d):

v1 = {0.5, 0.5, 0.5}
n1 = {1, 1, 1}
v2 = {1, 1.5, 0}
n2 = {1, 1, 0}
d = 4


then we create cylinders out of the disks, with end-points of each ceylinder given by $$\pm 5 \mathbf n_i$$ to respective center position of disk $$i$$:

cyl1 = Cylinder[{v1 - 5*n1, v1 + 5*n1}, d/2]
cyl2 = Cylinder[{v2 - 5*n2, v2 + 5*n2}, d/2]


and visualizing Graphics3D[{Opacity[.5], cyl1, cyl2}]:

But I don't know how much this approach helps in computing the overlap area of interest (and if computationally feasible).

• This requires a precise definition of "intersection area" that is more general than the given example. Commented Nov 25, 2019 at 12:59
• @HenrikSchumacher Good point! Sorry I was a bit unclear about this in my OP. What I meant is in essence: the orthogonal projection of a surface onto a plane. In this case, the surface is that of one of the two disks, and the plane is the plane of the other disk (the plane in which its surface lies). For related example, here in math SE the case of projecting the surface of a doughnut onto a plane is discussed, in our case it is always 2D disks in $3D.$ I hope this clarifies the aim of the question.
– user52181
Commented Nov 25, 2019 at 14:50
• Hmm. So this notion of intersection are is not symmetric in the two disks? Commented Nov 25, 2019 at 14:53
• @HenrikSchumacher It should be ultimately (to get the desired overlap quantification described in the OP). So I figure, by taking the product of the projection of disk $i$ onto plane of disk $j$ with that of disk $j$ onto plane of disk $i$ will yield a symmetrized intersection area definition.
– user52181
Commented Nov 25, 2019 at 14:54
• I understood the problem as an occlusion between objects or of a projected shadow of a disk over the other. In this case a point of view or a light source and a perspective are required. Commented Nov 25, 2019 at 16:56

Here's my take on solving it algebraically:

DiskRadius[Disk3D[_, _, radius_]] := radius;
RotateZToNormal[Disk3D[_, n_, _]] := RotationTransform[{{0, 0, 1}, n}];
MoveToDiskCenter[Disk3D[p_, _, _]] := TranslationTransform[p];
TransformUnitDiskTo[d_Disk3D] := RightComposition[RotateZToNormal[d], MoveToDiskCenter[d]]
Project2D = Most;(*leave out z component to project into 2D*)
CartesianFromPolar = (# /. {r -> Sqrt[x^2 + y^2], \[Phi] -> ArcTan[x, y]} &);
UnitDiskToProjectedEllipseTransform[to_Disk3D] := Function[from,
Composition[
AffineTransform[Reverse[##]] &, (*construct 2d affine transform unitdisk -> projected disk/ellipse*)
CoefficientArrays[#, {x, y}] &, (*extract 2d ellipse linear transformation coefficients*)
Simplify, CartesianFromPolar, Project2D,
InverseFunction[TransformUnitDiskTo[to]],
TransformUnitDiskTo[from]
][r {Cos[\[Phi]], Sin[\[Phi]], 0}]
]
ProjectDiskOnto[to_Disk3D] := Function[from,
Composition[
InverseFunction[UnitDiskToProjectedEllipseTransform[to][from]]
][{x, y}]
]
ProjectedDiskRegion[to_Disk3D] := Function[from,
RegionIntersection[
(ImplicitRegion[#, {x, y}] & @* ProjectDiskOnto[to]) /@ {to, from}
]
]
DiskIntersectionArea[disk1_Disk3D, disk2_Disk3D] :=
Mean[Area /@ {ProjectedDiskRegion[disk1][disk2],
ProjectedDiskRegion[disk2][disk1]}]


Let's have a look at an example:

d1 = Disk3D[{-2, 3, -3}, {2, -3, 6}/7, 1]
d2 = Disk3D[{-1, 1, 1}, {-1, 2, -2}/3, 4/5]


Here we encoded the disks by their center point, their normal and their radius with a custom Disk3D head. We can plot these to get an idea

PlotDisk3D[d_Disk3D] := ParametricPlot3D[
TransformUnitDiskTo[d][r {Cos[\[Phi]], Sin[\[Phi]], 0}],
{r, 0, DiskRadius[d]}, {\[Phi], 0, 2 \[Pi]}, Mesh -> None
]
Show[PlotDisk3D /@ {d1, d2}]


The idea of the solution is to first get an implicit 2d equation of each disk transformed into the reference frame of the other disk and then project it into the xy plane. We do that by creating the function TransformUnitDiskTo which produces an AffineTransform that would transform a unit disk sitting in the xy-plane into any given to disk. Next we start with a parametric polar representation of a unit disk, which we first transform into our (from) disk which we want to project, and then follow it by an inverse affine transform to get it into the reference frame of our to disk. In this reference frame we can project it into 2D and after that convert back to cartesian coordinates and into an implicit representation instead of parametric. Our two example disks in the other reference frame now look like this:

ProjectDiskOnto[d1][d2]


$$\left(\frac{43 x}{52}-\frac{8 y}{13}+\frac{1}{4}\right)^2+\left(\frac{37 x}{65}+\frac{11 y}{13}+\frac{1}{5}\right)^2\leq \frac{16}{25}$$

ProjectDiskOnto[d2][d1]


$$\left(\frac{11 x}{13}-\frac{8 y}{13}+\frac{18}{13}\right)^2+\left(\frac{37 x}{65}+\frac{43 y}{52}-\frac{47}{65}\right)^2\leq 1$$

Projecting a disk onto itself naturally always gives back the disk unaltered:

ProjectDiskOnto[d1][d1]


$$x^2+y^2\leq 1$$

ProjectDiskOnto[d2][d2]


$$x^2+y^2\leq \frac{16}{25}$$

Now we can perform the region intersection inside ImplicitRegions

and finally take the average of the Region Areas, which Mathematica happily performs for us symbolically and we end up with an exact expression, which we can either simplify a bit via RootReduce on the algebraic parts or just get a numeric approximation with desired accuracy:

DiskIntersectionArea[d1, d2] // N
(* 0.9875 *)


# Update

• Fixed a bug in the construction of the implicit ellipse representation.
• I like this solution a lot. Commented Nov 28, 2019 at 23:20
• @user929304 i) this is true when our disk normal is {0,0,1}. Generally the projection always happens along the disk normal, which means you have to translate a disk along the other disks normal to keep the projected area constant. ii) No, it should be easily possible to extend it to disks with radii different from one. I'll try to expand the code when i have a bit more free time. Commented Nov 29, 2019 at 17:10
• @user929304 iii) I think mainly you need to be consistent with your units. That means you either could convert the radius into the same unnormalized units as you use for positions and multiply the resulting area with your normalization factor^2, or you could scale your unnormalized positions to 'real' position units before you feed it to the area calculation. Commented Nov 29, 2019 at 17:12
• Thanks a lot for the update Thies! This indeed seems to work really well and results are in agreement with e.g. Henrik's approach. There are sometimes cases where the algorithm pauses longer than usual and returns a non-numeric result, e.g. with this example: d1test = Disk3D[{-148.523, -19.2954, 142.73}, {0.696493, -0.319583, -0.642467}, 12.5]; d2test = Disk3D[{148.239, -43.1213, -149.823}, {-0.415674, 0.313735, 0.85369}, 12.5]; DiskIntersectionArea[d1test, d2test] // N Do you know what might be causing this?
– user52181
Commented Dec 5, 2019 at 17:37
• @user929304 Oh, thanks for pointing that out, the code wasn't robust enough to handle the empty RegionIntersection case. I updated ProjectedDiskRegion to support those cases properly, so now it should return EmptyRegion[2] instead and DiskIntersectionArea returns 0 in those cases. Commented Dec 6, 2019 at 16:20

A circle with center at $$p$$ radius $$r$$ and orientation $$\vec n$$ normalized, can be represented as

$$c_i \to \{p_i, r_i \vec n_i\}$$ A circle

$$c_0 \to \{p_0, r_0 \vec n_0\}$$

can be drawn with the parametric

$$p = p_0 + r_0 \vec e_1\cos\mu +r_0\vec e_2\sin\mu,\ \ \mu\in (0,2\pi]$$

where $$\vec n_0, \vec e_1, \vec e_2,$$ form an orthonormal basis.

An ellipse can be drawn with the parametric

$$p = p_0 + a_0 \vec e_1\cos\mu +b_0\vec e_2\sin\mu,\ \ \mu\in (0,2\pi]$$

where $$a_0, b_0$$ are the main axes.

With those facts follow two modules: one which plot a circle drawcircle and other that giving two circles, draw the projection for a circle $$\{p_1,r_1\vec n_1\}$$ onto the plane defined by $$\{p_2,\vec n_2\}$$ projection

With those modules we can verify the projection overlapping. The overlapping area computation is not considered here.

NOTE

The projections are represented by dashed lines.

drawcircle[p_, n_, color_] := Module[{t1, t2, t3, tau, tau0, n0, v, gr1, gr2, equ1, equ2, equ3, sols, r},
r = Norm[n];
n0 = n/r;
tau = {t1, t2, t3};
equ1 = n0.tau == 0;
equ2 = tau.tau == 1;
equ3 = t1 + t2 + t3 == 0;
sols = Quiet@Solve[{equ1, equ2, equ3}, tau][[1]];
tau0 = tau /. sols;
v = Normalize[Cross[tau0, n]];
gr1 = ParametricPlot3D[p + r tau0 Cos[mu] + r v Sin[mu], {mu, 0, 2 Pi}, PlotStyle -> color];
gr2 = Graphics3D[{Black, Arrow[{p, n + p}]}, Axes -> True,
Boxed -> False];
Return[Show[gr1, gr2]]
]

projection[p1_, n1_, p2_, n2_, color_] := Module[{n10, n20, v1, v2, e1, e2, ex, ey, ez, equ1, equ2, equ3, sols,e20, p1p2, grep, axis, mu, r1, r2},
n10 = Normalize[n1];
r1 = Norm[n1];
r2 = Norm[n2];
n10 = n1/r1;
n20 = n2/r2;
v1 = (n10.n20) n20;
v2 = n10 - v1;
e1 = Normalize[v2];
e2 = {ex, ey, ez};
equ1 = e2.n2 == 0;
equ2 = e1.e2 == 0;
equ3 = e2.e2 == 1;
sols = Quiet@Solve[{equ1, equ2, equ3}, e2][[1]];
e20 = e2 /. sols;
p1p2 = p1 - ((p1 - p2).n2) n2;
grep = ParametricPlot3D[p1p2 + r1 n1.n2 e1 Cos[t] + r1 e20 Sin[t], {t, 0, 2 Pi}, PlotStyle -> {Dashed, color}];
axis = ParametricPlot3D[(1 - mu) p1p2 + mu p1, {mu, 0, 1}, PlotStyle -> {Dotted, Gray}];
Return[Show[grep, axis]]
]

p1 = {1/2, -1/2, 1/2};
p2 = {0, 0, 0};
n1 = Normalize[{1, 0, 1}];
n2 = Normalize[{0, 1, -1}];

c1 = drawcircle[p1, n1, Blue];
c2 = drawcircle[p2, n2, Red];
Show[c1, c2, PlotRange -> All, AspectRatio -> 1]

pc1c2 = projection[p1, n1, p2, n2, Blue];
Show[c1, c2, pc1c2, PlotRange -> All, AspectRatio -> 1]

pc2c1 = projection[p2, n2, p1, n1, Red];
Show[c1, c2, pc2c1, PlotRange -> All, AspectRatio -> 1]


Maybe this function does the trick:

F[{p1_, n1_, r1_}, {p2_, n2_, r2_}] :=
Block[{A1, A2, v, w1, w2, area1, area2, angle},
A1 = Orthogonalize[Join[{n1}, IdentityMatrix[3]]][[2 ;; 3]];
A2 = Orthogonalize[Join[{n2}, IdentityMatrix[3]]][[2 ;; 3]];
angle = Min[VectorAngle[n1, n2], VectorAngle[n1, -n2]];
If[angle < 1. 10^-12,
{w1, w2} = r1 IdentityMatrix[2]
,
v = r1 Normalize@Cross[n1, n2];
w1 = A2.v;
w2 = A2.Cross[n1, v];
];
area1 = Area@RegionIntersection[
Ellipsoid[A2.(p1 - p2), KroneckerProduct[w1, w1] + KroneckerProduct[w2, w2]],
Disk[{0, 0}, r2]
];
If[angle < 1. 10^-12,
{w1, w2} = r2 IdentityMatrix[2]
,
v = r2 Normalize@Cross[n2, n1];
w1 = A1.v;
w2 = A1.Cross[n2, v];
];
area2 = Area@RegionIntersection[
Ellipsoid[A1.(p2 - p1), KroneckerProduct[w1, w1] + KroneckerProduct[w2, w2]],
Disk[{0, 0}, r1]
];
(area1 + area2)/2
]


The procedure works like this:

First, we compute two linear projections A1 and A2 onto the plane, one with nullspace spanned by the normal n1, the other spanned by the normal n2 (and each projections is isometric in the orthogonal complement of its nullspace).

For defining an Ellipsoid, it suffices to know its center and its covariance matrix. The latter can then be obtained from the sum of the Kronecker products of the half-axis vectors. And half-axis vectors (w1 and w2) can be computed as in the code above.

After having projected everything into the 2D-plane, we can simply apply RegionIntersection and Area. Here, the intersection and its area are computed numerically, which might make this a rather slow approach (if you have to run this computation a zillion times). Maybe the intersection area for a Disk and a general Ellipsoid can be computed symbolically and offline? Then one can make this faster. But I don't really know how to so that. Maybe this might help.

# Edit

I have updated the code above to treat the cases when n1 and n2 are almost parallel in a different way. I have not tested it for correctness and for robustness, yet. So use it with caution.

• Many thanks Henrik, this seems to work out quite well and I really like your approach! I'm gonna try it on a set of different examples and get back to you! Quick question, so r1, r2 are here our diameters right? in case (though happens rarely) the normals are parallel or exactly the same the algorithm has an exception as v becomes zero. But slightly nudging one component fixes it (so e.g. n1=n2={1, 1, 0}; crashes but n1={1, 1, 0}; and n2={0.99, 1, 0}; is sufficient to workaround it.
– user52181
Commented Nov 30, 2019 at 15:10
• r1 and r1 are of course the radii. ;) Yeah, one certainly has to work out some degenerate cases. However, that's quite simple: If n1 and n2 are almost parallel or antiparallel, you may choose {w1, w2} = r1 IdentityMatrix[2] and {w1, w2} = r2 IdentityMatrix[2]. Notice also that I assumed the normals n1 and n2 to be normalized. If your input vectors are not normalized, just normalize them with Normalize. Commented Nov 30, 2019 at 15:20
• The w-vectors are the half-axis vectors of the ellipsoid that is created by projection one circle onto the plain of the orther (and than projected to $\mathbb{R}^2$ via the appropriate A-matrix). If the planes of the circles are parallel, then this ellipsoid is a round circle and actually every two vectors that are orthogonal to each other and have the approriate length may serve as w-vectors. Commented Dec 1, 2019 at 15:06
• Thank you for getting back to me so promptly, and for the edit (I've tested it, works well ;) )! Just to sanity check my understanding, here are my grouped thoughts 1) to find the projected circle (ellipsoid), the diameters of the circle are projected onto the diameters of the ellipsoid. 2) To find the large axis vector of the corresponding ellipsoid, we have to find the vector along the principal line of the plane of n2 at which the two planes intersect, therefore, the large diameter of the ellipsoid has the same length as the diameter of the circle. The intersection line of the(..)
– user52181
Commented Dec 2, 2019 at 11:45
• is by definition orthogonal to both normal vectors (because it belongs to both planes), thus $\mathbf w_1 \propto \mathbf n_1 \times \mathbf n_2.$ 3) then comes the minor axis vector $\mathbf w_2$, which is $\perp \mathbf w_1,$ and it must be perpendicular to $\mathbf n_2$ because it lies in the plane to which $\mathbf n_2$ is normal. So I would have thought, $\mathbf w_2 \propto \mathbf n_2\times \mathbf v,$ instead we used (?) $\propto \mathbf n_1 \times \mathbf v$ which I don't get yet :/ 4) Once w1, w2 computed, we collapse them to 2D and (...)
– user52181
Commented Dec 2, 2019 at 11:46

This seems like it ought to work, but throws an exception right at the end. Building on the the example in the question, but with the length of the cylinders defined by a variable len that I can play with:

v1  = {0.5, 0.5, 0.5};
n1  = {1, 1, 1};
v2  = {1, 1.5, 0};
n2  = {1, 1, 0};
d   = 4;
len = 10;


Draw the cylinders and also a hperplane located at the origin of one of the cylinders.

cyl1 = Cylinder[{v1 - len*n1, v1 + len*n1}, d/2];
cyl2 = Cylinder[{v2 - len*n2, v2 + len*n2}, d/2];
hp   = Hyperplane[n1, v1];
Graphics3D[{Opacity[.5], cyl1, cyl2, hp}]


Now intersect them.

proj = RegionIntersection[Region@cyl1, Region@cyl2, Region@hp]


This is a proper region according to the test

RegionQ@proj
(* True *)


But it burps when I try to compute over it, using RegionMeasure[ ] or other methods.

RegionMeasure@proj


• This is a bug and it were good if you could report it to support. Commented Nov 29, 2019 at 15:49
• Small typo in the cyl definitions, they should be: cyl1 = Cylinder[{v1 - len*n1, v1 + len*n1}, d/2]; cyl2 = Cylinder[{v2 - len*n2, v2 + len*n2}, d/2]; (now corrected in the OP too)
– user52181
Commented Nov 29, 2019 at 15:51
• Updated to match... Commented Nov 29, 2019 at 16:04
• @user21, bug reported. Commented Nov 29, 2019 at 16:29
• FYI, Wolfram acknowledged the issue as a bug. Commented Dec 3, 2019 at 14:51

To see the orthogonal projection, just make the disks into cylinders in turn. Most of the time the projection is going to be the empty set, though.

• Apologies for being terse, I'm on mobile. Make disk 1 into a cylinder but leave disk 2 alone before using region intersect. Then change roles of disk 1 and disk 2. Commented Nov 26, 2019 at 16:32
• Dear Dave, could you please kindly expand upon your answer when time allows. As it stands it is somewhat unclear and I'm not sure I've interpreted it correctly.
– user52181
Commented Nov 29, 2019 at 12:09
• Nothing I tried panned out: geometric scene not set up for 3d, ellipses, or exact calculations. The projection onto a disc as you've described is going to be an ellipse that will almost certainly have an x y term. Wlig, let the disk being projected onto is the unit circle on the x y plane. Solving simultaneous non linear equations with the unit circle and the general ellipse surprisingly does give exact values for x and y. I didnt see any symmetry and see no obvious way to set up integrals to get an exact solution without case by case situations. Commented Dec 1, 2019 at 20:30