As you've pointed out, the only parameter that matters here is the angle $\theta$ between $x$ and $y$. To see how, consider instead the Gaussian integral:
$$
I(x,y)=\frac{1}{(2\pi)^{(d+1)/2}}\int_{u\in\mathbb{R}^{d+1}}\max(0,x^Tu)\cdot\max(0,y^Tu)\exp(-\frac{1}{2}u\cdot u)du
$$
The integral you are interested in is obtained by changing to $(d+1)$ dimensional spherical coordinates, and easily integrating over the radial coordinate. As such, it will be enough to just focus on the Gaussian integral.
I like computing these types of integrals by forming an orthonormal basis of $\mathbb{R}^{d+1}$ like so:
$$
b_1=(x+y)/\|x+y\|\\
b_2=(x-y)/\|x-y\|
$$
...and where the remaining basis elements are chosen to be orthogonal to $b_1$ and $b_2$. With respect to this basis we have:
$$
x=(a,b,0,...,0)\\
y=(a,-b,0,...,0)
$$
where $a=|\cos(\theta/2)|$ and $b=|\sin(\theta/2)|$. Of course the components of $u$ change as does the Gaussian measure, but the Gaussian measure is invariant under rotations so I'll suppress relabeling the components of $u$ with respect to the new basis. After easily integrating over the last $(d-1)$ coordinates of $u$ we write the integral $I(x,y)$ as:
$$
\frac{1}{2\pi}\int_{(u_1,u_2)\in R}(a\cdot u_1+b\cdot u_2)\cdot (a\cdot u_1-b\cdot u_2)\exp(-\frac{1}{2}(u_1^2+u_2^2))du_1du_2,
$$
where $R=\{(u_1,u_2)|(a\cdot u_1+b\cdot u_2)>0\textrm{ and }(a\cdot u_1-b\cdot u_2)>0\}$ is the region where the two maximums are nonzero. After a change of coordinates $v_1=au_1$ and $v_2=bu_2$ we obtain:
$$
\int_{(v_1,v_2)\in R'}\frac{\left(v_1^2-v_2^2\right) e^{\frac{1}{2}
\left(-\frac{v_1^2}{a^2}-\frac{v_2^2}{b^2}\right)}}{2 \pi a b}dv_1dv_2
$$
where $R'=\{(v_1,v_2)|(v_1+v_2)>0\textrm{ and }(v_1-v_2)>0\}$. From here, change to polar coordinates $(r,t)$ and integrate over the radial coordinate to obtain:
$$
\int_{-\pi/4}^{\pi/4}\frac{a^3 b^3 \cos (2 t)}{\pi \left(a^2 \sin ^2(t)+b^2 \cos ^2(t)\right)^2}dt
$$
Using Mathematica for example, we obtain:
$$
\int_{-\pi/4}^{\pi/4}\frac{a^3 b^3 \cos (2 t)}{\pi \left(a^2 \sin ^2(t)+b^2 \cos ^2(t)\right)^2}dt = \frac{a b+(a-b) (a+b) \tan ^{-1}\left(\frac{a}{b}\right)}{\pi }
$$
From here, you can then express everything in terms of the inner product of $x$ and $y$ along with the norms of $x$ and $y$.
Of course, hopefully the result can be checked! May I ask where this integral came up? I bet there are others who can obtain something similar using some pretty slick symmetry arguments.
Cheers!
