The problem:
Given some distribution $P(x, y, z)$ over the surface of the sphere $S^2$, what is the distribution of lines $L$ generated by choosing two points $p_1, p_2 \sim P$ and returning the line $p_1 + \lambda(p_2 - p_1)$? Additionally, if I have some 2-dimensional surface $\mathcal{M} \subset \mathbb{R}^3$ defined by the equation $F(x, y, z) = 0$, what is the distribution on this surface of the intersection with lines $L$ chosen by the previous procedure?
I'm also interested in a generalization of this where $p_1$ and $p_2$ are chosen from different distributions $P_1, P_2$ on $S^2$.
Some more details on the problem: this question is related to the Crofton Formula https://en.wikipedia.org/wiki/Crofton_formula which connects the surface area of a surface in $\mathbb{R}^n$ to the number of intersections with 'randomly' chosen lines. In the paper https://orca.cardiff.ac.uk/id/eprint/1813/1/area.pdf they use this to estimate the surface area of volumes in $\mathbb{R}^3$ by intersecting it with random lines. I am trying to generalize this to distributions on the surface $\mathcal{M}$ which are not uniform.
Here is what I have tried so far:
Given a point $p = (x, y, z) \in \mathcal{M}$ on our surface, denote by $\hat{n}(p)$ the unit normal at $p$ given by $\frac{\nabla F}{|\nabla F|}$. Denote by $\hat{n}_S(s)$ the unit normal vector to the sphere at the point $s$, i.e. $\frac{s}{|s|}$. Let $s_1, s_2$ be points on $S^2$ such that the line $s_1 \to s_2$ intersects the point $p$. Denote by $Q(p)$ the density function on the surface $\mathcal{M}$ generated by this sampling procedure.
We have
$$ Q(p) = \int_{s_1 \in S^2} \int_{s_2 \in S^2} \left( \frac{(s_1 - s_2) \cdot \hat{n}(p)}{|(s_1 - s_2)| \cdot |\hat{n}(p)|} \right) \cdot \left( \frac{(s_1 - s_2) \cdot \hat{n}_S(s_1)}{|(s_1 - s_2)| \cdot |\hat{n}_S(s_1)|} \right) \cdot \left( \frac{|s_1 - p||s_2 - p|}{|s_1 - s_2|^2} \right) \cdot P(s_1) \cdot P(s_2) ds_1 ds_2 $$
The term
$$ \frac{(s_1 - s_2) \cdot \hat{n}(p)}{|(s_1 - s_2)| \cdot |\hat{n}(p)|} $$
represents the value $\cos(\theta)$ where $\theta$ is the angle between the line $L$ and the unit normal vector of the surface.
Similarly, the term
$$ \frac{(s_1 - s_2) \cdot \hat{n}_S(s_1)}{|(s_1 - s_2)| \cdot |\hat{n}_S(s_1)|} $$
represents the cosine of the angle between the line $L$ and the normal to the sphere at $s_1$ and $s_2$ (this angle is the same for both by symmetry.
Finally,
$$ \frac{|s_1 - p||s_2 - p|}{|s_1 - s_2|^2} $$
represents the ratio of the distances from $s_1$ to $p$ and $s_2$ to $p$. If we consider a right-angled triangle with sides $(s_1, s_2), (s_1, s_1 + ds_1), (s_1 + ds_1, s_2)$ and the congruent triangle $(p, s_2), (p, p + dp), (p + dp, s_2)$, we can see that the ratio $\frac{dp}{ds_1}$ is equal to
$$ \frac{|s_1 - p|}{|s_1 - s_2|} $$
By symmetry we get a similar term for $s_2$, and the product gives us the term above.
Why this doesn't seem to work:
I have tried generating a discrete approximation to this by taking a (roughly) uniform set $A$ of points on the surface of $S^2$ using a fibonacci lattice, and a uniformly distributed set of points on the surface of a torus.
For each point $p$ on the torus, and for each point $s_1 \in U$, I find the point $s_2$ which is closest to the intersection of the line $s_1 \to p$ with the sphere $S^2$. Using these pairs $(s_1, s_2)$ I then compute the integrand above and sum over all pairs to approximate the density at the point $p$.
This should generate roughly uniform values $Q(p)$ however the distribution has some strange correlations with the x, y, and z coordinates (see scatter plots of x/y/z vs. estimated $Q(p)$ below):
