1
$\begingroup$

I am looking for a way to sample points uniformly around a particular point on the unit sphere.

Working on spherical coordinates we can express any point on the unit sphere as $(\sin\theta \cos\phi, \sin\theta\sin\phi,\cos\theta)$. Say the initial point $r_0$ has coordinates $(\sin\theta_0 \cos\phi_0, \sin\theta_0\sin\phi_0,\cos\theta_0)$, I want to sample points uniformly over the solid angle of a cone around the initial point with an apex angle of $2\delta$ (see Figure below).

Cone around arbitrary direction

If the initial point is the north pole this is easy to do because $\phi$ can be taken uniformly as $2\pi \,\rm{rand}(0,1)$ and for $\theta$ the usual inverse transform method can be truncated to the desired range $(0,\delta]$ see wiki.

I don't know how to go around the general case though, I'd appreciate if anyone has some idea of how to attack this problem.

Thanks in advance!

Improve this question
$\endgroup$
4
  • 2
    $\begingroup$ The only thing missing was performing a rotation afterwards to the north pole case. Check math.stackexchange.com/questions/56784/… for more details $\endgroup$
    – David Leonardo Ramos
    Commented May 24, 2022 at 21:00
  • 1
    $\begingroup$ If we want the points to be uniformly distributed on the surface, it is not the polar angle $\phi$ that should be uniformly distributed but $\cos(\phi)$. $\endgroup$
    – Maxim Umansky
    Commented May 24, 2022 at 21:54
  • 1
    $\begingroup$ Details are given in Brannon's book on rotations: iopscience.iop.org/book/978-0-7503-1454-1/chapter/… $\endgroup$
    – Biswajit Banerjee
    Commented May 25, 2022 at 1:28
  • 1
    $\begingroup$ Uniformly sample over the sphere, and then reject the samples outside your desired region. $\endgroup$
    – ZR Han
    Commented May 27, 2022 at 0:55

1 Answer 1

Reset to default
1
$\begingroup$

So if I correctly understand, you know how to sample in the "vertical" cone, and you ask how to rotate?

Here is some pseudo-code to get the unit quaternion representing a rotation sending a line directed by a vector U to a line directed by a vector V:

getTheQuaternion <- function(U, V){
  U <- normalize(U)
  V <- normalize(V)
  d <- dotproduct(U, V)
  c <- sqrt(1 + d)
  r <- 1 / sqrt(2) / c
  W <- r * crossproduct(U, V)
  quaternion <- (c/sqrt(2), W[1], W[2], W[3])
  return(quaternion)
}

and some pseudo-code to get the rotation matrix corresponding to a unit quaternion:

quaternion2matrix <- function(q){
  (
    (1 - 2*q[3]^2 - 2*q[4]^2, 2*q[2]*q[3] - 2*q[4]*q[1], 2*q[2]*q[4] + 2*q[3]*q[1]),
    (2*q[2]*q[3] + 2*q[4]*q[1], 1 - 2*q[2]^2 - 2*q[4]^2, 2*q[3]*q[4] - 2*q[2]*q[1]),
    (2*q[2]*q[4] - 2*q[3]*q[1], 2*q[3]*q[4] + 2*q[2]*q[1], 1 - 2*q[2]^2 - 2*q[3]^2)
  )
}
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.