SamplingUniformSphere.py

#http://stackoverflow.com/questions/15880367/python-uniform-distribution-of-points-on-4-dimensional-sphere
dim = 3

norm = np.random.normal
normal_deviates = norm(size=(dim, N))

radius = np.sqrt((normal_deviates**2).sum(axis=0))
points = normal_deviates/radius
points = points.reshape(-1,3)
t = np.array( [ 0.0, 0.5, 0.5 ] )
c = np.array( [ 10, 10, 10 ] )

grid = np.mgrid[int( -0.5*imshape[0] ):int( 0.5*imshape[0]),int( -0.5*imshape[1] ):int( 0.5*imshape[1]),int( -0.5*imshape[2] ):int( 0.5*imshape[2]) ]

#x = np.linspace(int( -0.5*imshape[0] ), int( 0.5*imshape[0] ), imshape[0] ).astype(int)
#y = np.linspace(int( -0.5*imshape[1] ), int( 0.5*imshape[1] ), imshape[1] ).astype(int)
#z = np.linspace(int( -0.5*imshape[2] ), int( 0.5*imshape[2] ), imshape[2] ).astype(int)

#x1 = np.unique( ( x*m[0] + cmass[0] ).astype(int) )
#y1 = np.unique( ( y*m[1] + cmass[1] ).astype(int) )
#z1 = np.unique( ( z*m[2] + cmass[2] ).astype(int) )

p = c + t

a = np.array( [ (c[0],c[1]), (p[0], p[1]) ]  )
b = np.array( [ c[2], p[2] ] )

print np.linalg.lstsq( a, b )


t = grid * m[:,np.newaxis,np.newaxis,np.newaxis] +  c[:,np.newaxis,np.newaxis,np.newaxis]

#tvals = zip(*(x.flat for x in t))

#ray[] = 1

#raysfile = op.join( sim_dir, 'template/rays.nii.gz' )
#nib.save( nib.Nifti1Image( ray, bmaskimg.get_affine(), bmaskimg.get_header() ), raysfile )

places = np.where(t==grid, 1, 0)
z = places[0] * places[1]

print np.shape(places)
print np.shape(z)

raysfile = op.join( sim_dir, 'template/rays.nii.gz' )
nib.save( nib.Nifti1Image( z.astype(uint8), bmaskimg.get_affine() ), raysfile )

spots = np.zeros( shape=bmask.shape )
for cidx in saveidx:
    spots[cidx] = 1
nib.save( nib.Nifti1Image( spots, bmaskimg.get_affine(), bmaskimg.get_header() ), fixedskullspotsfile )