I have not read the comment by B. Smith.

Instead i am going for an analytic solution:

The volume of a unit sphere is 4 pi/3

Imagine the function "fraction" f(z1,z2) where z1 and z2 are axial slices like latitude circles on the unit sphere between -0.5 and 0.5, and f is the factional volume between these latitudes with f(-.5,.5)=1.

Now do a double volume integral that considers every point in the volume against every other point. Then align the line made by these two points with a coordinate system that has the z-axis parallel. Find f(z1,z2) and keep averaging this fraction.

The double volume integral can be accomplished in cartesian coordinates. 

The values of z1 and z2 need not be the result of a coordinate transform, they can simply be distances from the sphere center, d1,d1. We will use z in a different context hereafter: the z position of one of the end points. d1,d2 can be arrived at by 1) noting the vector form of the line segment p1 = (x1,y1,z1) and p2 = (x2,y2,z2) and finding the parallel line passing through the origin, q= (0,0,0).

We find the orthogonal vector v from the line segment ends to the parallel axis line by subtracting off the parallel component: v-