Computing the intersection area of two disks in 3D

Question

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?

---

Additional clarifications after comments:

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).

Answer 1

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[
    #.# <= DiskRadius[from]^2 &, 
    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

• Added support for arbitrary disk radii.
• Fixed a bug in the construction of the implicit ellipse representation.

Answer 2

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]

Answer 3

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.

Answer 4

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

Answer 5

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.