how to generate new points as offset with gaussian distribution for some points in spherical coordinates in python

StackOverflow https://stackoverflow.com/questions/22439193

  •  15-06-2023
  •  | 
  •  

Question

I am working with some points in spherical coordinates. I need to generate new points as the error points for them and a kind of offset for the old points. The new point should be in a specific distance from the old one which distributing by gaussian distribution. The angle of new point compared to old one is not important.I am trying to generate new points for r direction. no matter what are phi and theta (Spherical coordinates)

To generate the new point distributing by gaussian function, I tried the numpy.rand.normal(mean,std,..). But It is generating 1D random points over mean value and this mean value is a real number. In my case I need an approach to specify the position of the old point and I have one given standard deviation for this distance from the original points.

Honesty, I dont have a copy of my code. It is on the university's server. But let's assume I have an array of size 100*3 including the spherical (or cartesian) coordinates of some points on a surface of a cylinder. In spherical case, the first column presents the radius value, the second column is theta and third one shows the phi for the points. now I want to generate random points from them using gaussian distribution. there is a given standard deviation for the gaussian distribution. The only important thing is that the new points generated by gaussian distribution are limited in r value. No matter the position of points in term of theta and phi. When I tried numpy.rand.normal(mean,std,..), this generate some random points over the mean value. It does not help me. I want new points over my old ones with the given STD.

any idea would be appreciated.

This is a code, similar to mine written By Ophion How to generate regular points on cylindrical surface

def make_cylinder(radius, length, nlength, alpha, nalpha, center, orientation):

#Create the length array
I = np.linspace(0, length, nlength)

#Create alpha array avoid duplication of endpoints
#Conditional should be changed to meet your requirements
if int(alpha) == 360:
    A = np.linspace(0, alpha, num=nalpha, endpoint=False)/180*np.pi
else:
    A = np.linspace(0, alpha, num=nalpha)/180*np.pi

#Calculate X and Y
X = radius * np.cos(A)
Y = radius * np.sin(A)

#Tile/repeat indices so all unique pairs are present
pz = np.tile(I, nalpha)
px = np.repeat(X, nlength)
py = np.repeat(Y, nlength)

points = np.vstack(( pz, px, py )).T

#Shift to center
shift = np.array(center) - np.mean(points, axis=0)
points += shift

#Orient tube to new vector

#Grabbed from an old unutbu answer
def rotation_matrix(axis,theta):
    a = np.cos(theta/2)
    b,c,d = -axis*np.sin(theta/2)
    return np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)],
                     [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)],
                     [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]])

ovec = orientation / np.linalg.norm(orientation)
cylvec = np.array([1,0,0])

if np.allclose(cylvec, ovec):
    return points

#Get orthogonal axis and rotation
oaxis = np.cross(ovec, cylvec)
rot = np.arccos(np.dot(ovec, cylvec))

R = rotation_matrix(oaxis, rot)
return points.dot(R)

now calling the function:

points = make_cylinder(3, 5, 5, 360, 10, [0,2,0], [1,0,0])

sigma = 0.5 # given STD ossfet_points = numpy.random.normal(np.mean(point[:,0]), sigma, size=(n,3))

Was it helpful?

Solution

If I'm not mistaken, you want random points on a spherical manifold with a gaussian distribution of distances from the center. If so, then you have the latter problem solved by sampling gaussian values of the radius using numpy.rand.normal

To get random spherical points is a little bit more tricky, but here's some code to do it (and a description of the math behind it at Wolfram MathWorld):

import numpy as np

num_points = 500
U = np.random.random(num_points)
V = np.random.random(num_points)

import math as m

def spherical_to_cartesian(vec):                                                  
    '''
    Convert spherical polar coordinates to cartesian coordinates:                       

    See the definition of spherical_cartesian_to_polar.                                 

    @param vec: A vector of the 3 polar coordinates (r, u, v)                           
    @return: (x, y, z)                                                                  
    ''' 
    (r, u, v) = vec                                                                     

    x = r * m.sin(u) * m.cos(v)                                                         
    y = r * m.sin(u) * m.sin(v)                                                         
    z = r * m.cos(u)                                                                    

    return [x, y, z]  

radius = 1.
points = np.array([spherical_to_cartesian([radius, 2 * np.pi * u, np.arccos(2*v - 1)]) for u,v in zip(U,V)])

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax = Axes3D(fig)

ax.plot(points[:,0], points[:,1], points[:,2], 'o')

Which will give you points like this:

enter image description here

Now if you want them to have normally distributed radii, you just need to substitute your randomly generated values in the list comprehension which uses the variable radius like this:

radii = np.random.normal(10, 3, 100)
points = np.array([spherical_to_cartesian([r, 2 * np.pi * u, np.arccos(2*v - 1)]) for r,u,v in zip(radii, U,V)])

Is this more or less what you're looking for?

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top