randomly generate points within specific sphere coordinates in R

Question

I have specific x,y,z coordinates. I want to generate random points within a sphere given x as the center and x2 from another data frame as the edge of the radius (therefore the distance from x to x2 would be the length of the radius of the sphere).

I've seen a lot of discussion about how to do this appropriately mathematically (randomly distribute the points to avoid clustering) and was able to compile the easiest examples here and here for sample R code.

I also found this [R package sphereplot] (https://cran.r-project.org/web/packages/sphereplot/sphereplot.pdf) which might be easier, but am having a hard time understanding how to apply it.

These are all good starting points but using the sample code below I'm unsure how to apply it to specific starting points/spherical coordinates?

set.seed(101)
n <- 50
theta <- runif(n,0,2*pi)
u <- runif(n,-1,1)
x <- sqrt(1-u^2)*cos(theta)
y <- sqrt(1-u^2)*sin(theta)
z <- u

Using just one set/row of x,y,z coordinates from my data frame:

x = -0.0684486861
y= 0.0125857380
z= 0.0201056441

x2= -0.0684486861
y2 = 0.0125857380
z2= -0.0228805516

I want x,y,z to be the center of the sphere and the distance to x2,y2,z2 to be the radius length/edge of the sphere. Then generate random points from within the sphere with x,y,z as the center.

Eventually, I'm trying to do this with 100 spheres to compare if all the points in the second set of coordinates move in similar angles/directions in space.

Thanks for the guidance.

I'm having trouble understanding your question. A sphere exists in three dimensions. A circle exists in two. How are you planning to use a single `x` coordinate as "the center of the circle" or the center of a sphere? It takes two coordinate values for the former, three coordinate values for the latter. — pjs, Apr 07 '19 at 17:05
Well, I do have x,y,z coordinates per point -which I just added (I should mention these were generated from an ordination) but can't you just use the radius? Similar to calculating the volume of a sphere from the radius? — KNN, Apr 07 '19 at 17:43
No, you can't just use the radius because the radius defines a circle or sphere relative to a point, and it takes two or three coordinates to define a point in 2-space or 3-space, respectively. There are an infinite number of circles or spheres that have the same radius and x-coordinate. — pjs, Apr 07 '19 at 17:46
Thanks @pjs for helping me improve this. I was thinking about it wrong. I made some corrections that hopefully clarifies what I am trying to do. — KNN, Apr 07 '19 at 21:41
Do you want the points uniformly distributed through the sphere, i.e., the expected number of points in a given volume is directly proportional to the volume? — pjs, Apr 07 '19 at 21:52
No, I don't think so. I want to generate say 100 random points for each sphere regardless of the volume. — KNN, Apr 07 '19 at 22:57
Yes, but when you say you want random points they're going to have some distribution. Where do you want those points to fall? On the surface? Clustered nearer the center? Uniformly throughout? There's still some ambiguity in your question. — pjs, Apr 07 '19 at 23:51
I guess uniformly then. Sorry this is a little beyond my knowledge... — KNN, Apr 08 '19 at 01:54

Severin Pappadeux · Accepted Answer · 2019-04-10T02:53:01.920

Well, lets split problem in several subproblems.

First, is generating points uniformly distributed on a sphere (either volumetrically, or on surface), with center at (0,0,0), and given radius. Following http://mathworld.wolfram.com/SpherePointPicking.html, and quite close to code you shown,

rsphere <- function(n, r = 1.0, surface_only = FALSE) {
    phi       <- runif(n, 0.0, 2.0 * pi)
    cos_theta <- runif(n, -1.0, 1.0)
    sin_theta <- sqrt((1.0-cos_theta)*(1.0+cos_theta))
    radius <- r
    if (surface_only == FALSE) {
        radius <- r * runif(n, 0.0, 1.0)^(1.0/3.0)
    }

    x <- radius * sin_theta * cos(phi)
    y <- radius * sin_theta * sin(phi)
    z <- radius * cos_theta

    cbind(x, y, z)
}

set.seed(312345)
sphere_points <- rsphere(10000)

Second problem - move those points to the center at point X

rsphere <- function(n, r = 1.0, surface_only = FALSE, center=cbind(Xx, Xy, Xz)) {
    ....
    cbind(x+center[1], y+center[2], z+center[3])
}

Third problem - compute radius given center at (Xx, Xy, Xz) and surface point(Yx, Yy, Yz))

radius <- sqrt((Xx-Yx)**2+(Xy-Yy)**2+(Xz-Yz)**2)

Combine them all together for a full satisfaction. Ok, now that you provided values for center and radius, let's put it all together

rsphere <- function(n, r = 1.0, surface_only = FALSE, center=cbind(0.0, 0.0, 0.0)) {
    phi       <- runif(n, 0.0, 2.0 * pi)
    cos_theta <- runif(n, -1.0, 1.0)
    sin_theta <- sqrt((1.0-cos_theta)*(1.0+cos_theta))
    radius <- r
    if (surface_only == FALSE) {
        radius <- r * runif(n, 0.0, 1.0)^(1.0/3.0)
    }

    x <- radius * sin_theta * cos(phi)
    y <- radius * sin_theta * sin(phi)
    z <- radius * cos_theta

    # if radius is fixed, we could check it
    # rr = sqrt(x^2+y^2+z^2)
    # print(rr)

    cbind(x+center[1], y+center[2], z+center[3])
}

x1 = -0.0684486861
y1 = 0.0125857380
z1 = 0.0201056441

x2 = -0.0684486861
y2 = 0.0125857380
z2 = -0.0228805516

R = sqrt((x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2)
print(R)

set.seed(32345)
sphere_points <- rsphere(100000, R, FALSE, cbind(x1, y1, z1))

How it looks like?

UPDATE

Generated 10 points each on surface and in the volume and printed it, radius=2 looks ok to me

# 10 points uniform on surface, supposed to have fixed radius
sphere_points <- rsphere(10, 2, TRUE, cbind(x1, y1, z1))
for (k in 1:10) {
    rr <- sqrt((sphere_points[k,1]-x1)^2+(sphere_points[k,2]-y1)^2+(sphere_points[k,3]-z1)^2)
    print(rr)
}

# 10 points uniform in the sphere, supposed to have varying radius
sphere_points <- rsphere(10, 2, FALSE, cbind(x1, y1, z1))
for (k in 1:10) {
    rr <- sqrt((sphere_points[k,1]-x1)^2+(sphere_points[k,2]-y1)^2+(sphere_points[k,3]-z1)^2)
    print(rr)
}

got

[1] 2
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2

and

[1] 1.32571
[1] 1.505066
[1] 1.255023
[1] 1.82773
[1] 1.219957
[1] 1.641258
[1] 1.881937
[1] 1.083975
[1] 0.4745712
[1] 1.900066

what did you get?

Thank you. I am still uncertain how to apply this to my specific x,y,z coordinates instead of the 0,0,0? I can't follow any of the code in parts 2 or 3 actually. Though that might be just because I am really slow to learn R... — KNN, Apr 07 '19 at 18:12
@KNN I've updated the code together with all data you provided, please take a look — Severin Pappadeux, Apr 07 '19 at 22:23
Nice to see you used the cube root on radius to scale the results uniformly! — pjs, Apr 08 '19 at 18:20
An alternative way to generate the points on the surface of the sphere is to generate a vector of three independent N(0,1)'s, then normalize their values by the length of the vector (which can then be further scaled to be distributed into the interior of the sphere). This approach generalizes to arbitrary dimensional hyperspheres, while the trigonometric version gets really ugly really fast. — pjs, Apr 08 '19 at 18:26
Another variant is similar to [Marsaglia's polar method](https://en.wikipedia.org/wiki/Marsaglia_polar_method): generate vectors of three U(-1,1) coordinates, and reject while the norm is > 1. Then scale by the norm of the vector to place the coordinates on the surface of a sphere. In two or three dimensions this can be faster than generating normals (which is why Marsaglia did it), but the probability of acceptance for a vector decreases as the dimensionality increases, so the Gaussian approach is more general. — pjs, Apr 08 '19 at 18:32
Thank you very much @SeverinPappadeux (and @pjs) for your time and help! This works for me. Really, really appreciated. Just a follow-up question: If I only want to randomize points all at the same distance from the center point to my second point (x2,y2,z2) how do I fix that distance? — KNN, Apr 09 '19 at 00:44
@KNN You're welcome. You just set third parameter, `surface only`, to TRUE and that shall fix the radius. I just added few lines of code, you could uncomment them and verify that indeed radius is fixed and points are randomly distributed on the surface of the sphere — Severin Pappadeux, Apr 09 '19 at 01:08
@SeverinPappadeux I don't think that worked correctly regarding fixing the radius. Shouldn't all the generated points have the same euclidean distance from the center? I tested a few generated points to the center with `dist()` and `method="euclidean"` and the distances were different — KNN, Apr 09 '19 at 22:55
@KNN Well, they do, at least to me. Please check update. What did you get? — Severin Pappadeux, Apr 10 '19 at 02:53

randomly generate points within specific sphere coordinates in R

1 Answers1