Surface area of a sphere sliced by two orthogonal planes

I'm a programmer doing a graphics effect, and I've arrived at a math problem that my old dumb brain can't quite figure out the best way to solve.

Here's the problem: take a sphere and slice off a bit of it with a plane. You get a sphere cap and sliced sphere. Now, take another plane that is exactly orthogonal to the slicing plane and slice off another sphere cap. What is the surface area of the sliced sphere left over?

If the two sphere caps don't intersect, then the problem is simple: just subtract the surface area of the two sphere caps from the surface area of the sphere, and you're done. However, if the two planes are close enough to the center that the two sphere caps intersect, then there's a sliver of the sphere that gets sliced off twice, and you have to add it back. How do you do that?

In case anyone is curious about the actual thing I'm trying to do: I want to bake ambient occlusion into a texture that I'm going to use on a voxel mesh. The idea with AO is that you take a point on a surface, and a hemisphere centered around the surface normal, and you check how much of that surface area is occluded by objects. In a voxel world, that essentially reduces to the problem above, since the only thing that could occlude is neighboring voxels. I've googled this and the stuff I find are some variation of this method which uses vertex coloring (essentially) and looks quite bad. I want to do it properly.

• Incidentally: this is the right place to post this, right? Not MathOverflow? It seems elementary enough not to belong there, but I'm not totally clear on the distinction. – Oskar Feb 5 at 11:40
• elementary enough? You must be a genius! – pendermath Feb 5 at 12:09
• I guess I meant "too elementary for MathOverflow" :) – Oskar Feb 5 at 12:12
• This is the right place to post this. It can be set up as an advanced calculus problem. See for example en.wikipedia.org/wiki/Surface_area . If I have more time later on I might write up the functions and how to calculate it. – quarague Feb 5 at 12:30
• @quarague i figured as much, and once upon a time i feel like i could do calculus problems like this without much trouble, but it's been a little while :) I would appreciate any help – Oskar Feb 5 at 12:43

Just to present another way with respect to the good answer above.

Let's refer to the attached 2-D orthogonal projection We take a sphere of unitary radius, where the two cuts are individuated by the planes $$x=a$$ and $$y=b$$, or equivalently by the latitudinal angles $$\alpha, \, \beta$$.

The surface of the required "partial spherical wedge" will be described by the infinitesimal arc $$d \phi$$, rotated around the circle of radius $$r$$, at constant latitude $$\phi$$, by the logitudinal arc $$2 r \theta$$.

We have the following relations among the parameters. $$\left\{ \matrix{ a = \cos \alpha \hfill \cr b = \sin \beta \hfill \cr r = \cos \phi \hfill \cr r\cos \theta = a \hfill \cr} \right.$$

So, clearly, we will have $$S(\alpha ,\beta ) = 2\int_{\,\phi \, = \,\beta }^{\;\alpha } {r\,\theta \,d\phi } = 2\int_{\,\phi \, = \,\beta }^{\;\alpha } {\cos \phi \arccos \left( {\cos \alpha /\cos \phi } \right)d\phi }$$ where, if we want to preserve the positive sign for the integral we shall impose the inequality $$\beta \le \alpha$$.

I do not know (and doubt) whether the integral might be put into a closed form. Most probably you shall compute it numerically.

• Can you do this integral in WolframAlpha for x>0.4, y>0.75, so we can compare? (The values are chosen as 1/20th of an 8-15-17 triangle.) – user210229 Feb 7 at 23:26
• I have numerically tested this formula and I can confirm that it gives the same results as Matt's answer for all $0 < a < 1$ and $0 < b < 1$. – Paul Feb 8 at 14:04
• Thank you so much for your assistance! I think i follow the reasoning here (and even if I don't, the final answer is easy enough to implement). You have been extremely helpful! – Oskar Feb 8 at 15:44
• @Oskar: glad to help, and thanks for the bounty. – G Cab Feb 8 at 19:14

Update: An exact formula for the area on a unit sphere with $$x>a$$, $$y>b$$ is $$A=\big(\frac 12-a-b\big)\pi + 2b\arctan(\frac ac)+2a\arctan(\frac bc)+\arctan(\frac{c^2-a^2 b^2}{2abc})$$ for positive $$a$$ and $$b$$ with $$a^2+b^2<1$$ and $$c=\sqrt{1-a^2-b^2}$$.

Finding this formula was difficult, requiring both some Mathematica and some guesswork for me. Fortunately we can verify it just by checking its values at endpoints and its changes in the middle.

In the limit as $$a\rightarrow 0^+$$, the last term goes to $$\pi/2$$, so this gives $$A=(1-b)\pi$$ and agrees with Archimedes's hatbox theorem. Similarly this is correct in the limit as $$b\rightarrow 0^+$$.

In the middle, $$\frac{d^2A}{da\ db}=\frac{2}{\sqrt{1-a^2-b^2}}=2\sqrt{1+f_x^2+f_y^2},$$ which is correctly twice the element of surface area for $$f(x,y)=\sqrt{1-x^2-y^2}$$, corresponding to the upper and lower hemispheres.

As an example, if $$a=2/5$$ and $$b=3/4$$, then $$c=\sqrt{111}/20$$ and $$A = \frac{-3\pi}{20} + \frac 45\arctan(\frac{15}{\sqrt{111}}) + \frac{3}{2}\arctan(\frac{8}{\sqrt{111}})-\arctan(\frac{4\sqrt{111}}{25}) \simeq .234333$$

• interesting formulation(+1) – G Cab Feb 7 at 22:39
• Thank you so much for your answer. I wish I could split the bounty between you and G Cab, but unfortunately that's not the case. I'm extremely grateful for you time and good work. I will study your answer detail, if for no other reason than that I just might learn something :) – Oskar Feb 8 at 15:47
• I tried to make a nice closed form. Got as far as $\pi(1-a-b)-2 \sin ^{-1}\left(b \sqrt{\frac{1}{\left(1-a^2\right) \left(a^2+b^2\right)}}\right)+2 b \tan ^{-1}\left(\frac{a}{\sqrt{1-a^2-b^2}}\right)+2 a \tan ^{-1}\left(\frac{b}{\sqrt{1-a^2-b^2}}\right)+ i \left( \log(1-b^2) + \log(a^2+b^2)-2 \log\left(a+b\sqrt{a^2+b^2-1}\right)\right)$ (under condition $a\ge0,b\ge0,a^2+b^2<1$) – Paul Feb 8 at 16:24
• @Paul, I would only trust a formula where all the terms are real -- here $\sqrt{a^2+b^2-1}$ is imaginary, and its not obvious that the whole expression inside the big parenthesis is purely imaginary to cancel the $i$ outside. Ideally the formula should show symmetry in $a$ and $b$ too. – user210229 Feb 8 at 17:12
• One reason it’s so hard to find the formula is that it’s singular on all three parts of the boundary region ($a=0, b=0, c=0$), and replacing $\arctan(x)$ with $\pi/2-\arctan(1/x)$ leads to a singularity on the interior curve with $ab=c$ instead. – user210229 Feb 9 at 17:32

Area of piece $$P$$ in $$\mathbb{S}^2$$ :

i) Recall Gauss-Bonnet Theorem : $$\int_{\Omega}\ K +\int_{\partial \Omega}\ k_g = 2\pi \chi(\Omega)$$ where $$k_g$$ is a geodesic curvature.

ii) Geodesic curvature on the boundary cap : $$C_z$$ is a cap wrt $$z=\delta$$ so that $$( {\rm area}\ C_z)_{=2\pi \delta} +\int_{\partial C_z}\ k_g =2\pi$$

Remark - Cone : And there is a cone $$T_z$$ which is tangent to $$\mathbb{S}^2$$ along $$\partial C_z$$. Note that a central angle of $$T_z$$ is $$2\pi(1-\delta )$$.

iii) Define two orthogonal planes $$z=\delta,\ y=-1+\varepsilon$$ in $$\mathbb{E}^3$$. Assume that two planes cut $$P$$ from $$\mathbb{S}^2$$

For $$\partial C_z$$, we have a parametrization $$c_z(t)= \sqrt{2\delta-\delta^2 }(\cos\ t,\sin\ t,0) +(0,0,1-\delta)$$

Similarly, for $$\partial C_y$$ we have $$c_y(s)=\sqrt{2\varepsilon -\varepsilon^2 }(\cos\ s,0,\sin\ s) +(0,-1+\varepsilon ,0)$$

Find intersection between $$c_z$$ and $$c_y$$ : $$\sin\ t_i= \frac{-1+\varepsilon }{\sqrt{2\delta -\delta^2 }},\ \pi

$$\sin\ s_i = \frac{1-\delta}{\sqrt{2\varepsilon -\varepsilon^2}},\ 0

Along smooth arc in $$\partial P$$, the integral of geodesic curvature is $$|t_1-t_2|(1-\delta) + |s_1-s_2|(1-\varepsilon ) \ (:=T)$$

iv) Consider angle between two tangents at intersection point :

$$c_z'(t_2)=\sqrt{2\delta-\delta^2 }(-\sin\ t_2,\cos\ t_2,0) = ( \varepsilon -1 ,\sqrt{ 2\delta-\delta^2 -(1-\varepsilon)^2 } ,0)$$

$$c_y'(s_2)=\sqrt{2\varepsilon -\varepsilon^2 }(-\sin\ s_2,0,\cos\ s_2) = (\delta-1,0,\sqrt{2\varepsilon- \varepsilon^2 -(1-\delta)^2})$$

Hence angle $$\theta$$ between two vectors is $$\cos\ \theta=\frac{(1-\varepsilon )(1-\delta )}{\sqrt{2\varepsilon - \varepsilon^2}\sqrt{2 \delta -\delta^2}}$$

Hence Gauss-Bonnet theorem implies that $${\rm area}\ P= 2\pi- 2\theta-T$$