[Closed] evenly distributed Points on Sphere

if we got evenly distribution on a sphere, what is next?

Sometimes it just takes starting from scratch… and someone that knows more than I do. As to whats next, I want to build a delaunay triangulation for the surface and eventually build a voronoi surface.

-Michael

-------------------------------------------------------------------------------
   -- Here is an implementation of the Saff and Kuijlaars algorithm for evenly distributed Points on Sphere
   -- Golden Section spiral
   
   delete objects--Lets clearout all the objects before we start
   
   n = 600
   inc = pi * (3.0 - sqrt(5.0))
   offs = 2.0 / N
   pointArray = #()
   for k=1 to n do (
       y = k * offs - 1.0 + (offs / 2)
       r = sqrt(1.0 - (y*y))
       phi = k * inc
       phiDeg = phi*57.2957795
       pnt = ([(cos(phiDeg)*r), (sin(phiDeg)*r),y])
       append pointArray pnt
   )
   geosphere radius:100 segments:10 wirecolor:gray
   for p in pointArray do (
       sphere pos:(p*100) radius:10.0 wirecolor:yellow
   )
   clearSelection()
   forceCompleteRedraw()
  -------------------------------------------------------------------------------

Here is a purely brute force way to approach it, slow and clunky but im sure someone here could make it better

https://vimeo.com/60932036

-Michael

	----------------------------------------------------------------------------------------
  	--	RandomPointsOnASphere
  	--
  	--	This is a brute force and hence slow way to get N# points on the surface
  	--	of the sphere but its fun to see the points being repulsed by their neighbors
  	--	
  	--	Here are the steps for the algorithm
  	--	1. N# paints are placed randomly within a cube
  	--	2. The points are then projected to the surface of the sphere
  	--	3. Each point is then checked against the others and the closest 
  	--		one is found and moved some amount to even out the distribution
  	--	4. The repulsion most likely moves the point off the surface of the sphere
  	--		so its important to re project the points back onto the surface
  	--	5. was and repeat until either the points are the minimal testing distance away 
  	--	   from one another or we run out of iterations
  	--
  	--	Michael Spaw
  	--	Morphographic.com
  	--	Started 03/01/13
  	--	ver 0.01 03/01/13			  
  	--
  	----------------------------------------------------------------------------------------
  
  	delete objects--Lets clearout all the objects before we start
  	aqua = (color 87 224 198)--Create a color to be used for the bounding shape
  	drawRad = 2.0 --the radius of the points
  	RepulsionRadius = 5
  	nPoints =300 --Total number of points created within the test
  	Iterations = 150 --Total attemps to move the points so they no longer collide
  	BoundingVolume = 100
  	halfVol = BoundingVolume*0.5 
  
  	initxpos = 50.0-drawrad --Inital x position for the circle minus some to keep it in the bounding rectangle
  	initypos = 50.0-drawrad --same for the y
  	initzpos = 50.0-drawrad --same for the y
  
  	delete objects--Lets clearout all the objects before we start
  
  	boundingbox= box length:BoundingVolume width:BoundingVolume height:BoundingVolume pos:[0,0,(-halfVol)] wirecolor:aqua  --a nice bounding rectangle for neatness
  	addmodifier boundingbox (Lattice Strut_Radius:0.1 Joint_Radius:0)
  	SphereProxy = GeoSphere radius:halfVol wirecolor:blue segments:8
  	centerPos = SphereProxy.center--set the center of the sphere
  
  	--Make n# points in the bounding area
  	points= #() 
  	for i = 1 to nPoints do(
  		p = (random [-initxpos,-initypos ,-initzpos] [initxpos,initypos,initzpos])
  		myPoint = GeoSphere pos:p radius:drawRad wirecolor:orange segments:3
  		 append points myPoint
  	)
  
  	--Now move all the points to the surface of the sphere
  	for i = 1 to points.count do (
  		pointVector = points[i].pos - centerPos
  		points[i].pos = normalize pointVector * halfVol + centerPos
  	)
  		
  	forceCompleteRedraw()
  	
  	
  	fn findClosestNodes points testPoint RepulsionRadius = 
  	(
  		closestNodes = #()
  		for i=1 to points.count do ( -- loop through the points to find any that are within the test radius of the point in question (obj)
  			if testPoint!=points[i] do ( --check to make sure we aren't comparing against ourself
  				(
  					currDist = distance testPoint.pos points[i].pos -- get the dist between comparison pairs
  					if (currDist < (RepulsionRadius*2)) do ( -- if the circles overlap, add to the return array
  						append closestNodes points[i]
  					)
  				)
  			)
  		)
  		closestNodes
  	)	
  
  		closestPoints =  findClosestNodes points points[1] RepulsionRadius
  
  	--Main relaxation loop
  	for j=1 to iterations do -- run for up to the total # of itterations
  	(
  		noChange = true --if tru we are going to have to do some work
  		for i = 1 to points.count do 
  			(
  			closestPoints =  findClosestNodes points points[i] RepulsionRadius
  			for k=1 to closestPoints.count do 
  				(
  					rad = RepulsionRadius*2
  					d = distance closestPoints[k].pos points[i].pos
  					if (d < rad) do 
  					(
  						noChange = false
  						currentpos = points[i].pos
  						v = (closestPoints[k].pos - points[i].pos) * -0.07
  						points[i].pos += v
  						newpos = points[i].pos
  						
  						--Now move all the points to the surface of the sphere
  						pointVector = points[i].pos - centerPos
  						points[i].pos = normalize pointVector * halfVol + centerPos
  					)
  				)
  			)
  		if noChange do 
  			(
  				format "Number of iterations: %
" j
  				exit
  			)
  			
  		format "Number of iterations: %
" j
  		forceCompleteRedraw()
  		windows.processPostedMessages()
  	)
  	
  	

Okay, next one in the row, been there, done that, once again long long time ago (not interactive, not polished, just a piece of code lost in my scripts folder; based on Saff and Kuiljaars’ paper as well), sure enough the next thread you post will be something about Voronoi partitioning

try destroyDialog fibSphere catch()
rollout fibSphere "Fibonacci Sphere"
(
	spinner spnRadius "Sphere Radius: " range:[0,1e9,100]
	spinner spnPointCount "Point count: " type:#integer range:[0,1e9,1000]
	radiobuttons rbSphereOrCircle "Node Type:" labels:#("Circle", "Sphere") offset:[0,10]
	button btnCreate "CREATE POINTS" offset:[0,15]

	fn sphericalToCartesian radius theta phi = 
		radius * [cos phi * sin theta, sin phi * sin theta, cos theta]

	on btnCreate pressed do
	(
		local radius = spnRadius.value
		local totalPoints = spnPointCount.value
		local custNode = if rbSphereOrCircle.state == 1 then circle radius:5 render_displayRenderMesh:true else sphere radius:5

		local pointsRoot = (totalPoints as double)^.5
		local a = 1 - (double 1) / (totalPoints - 3)
		local b = (double .5) * (totalPoints + 1) / (totalPoints - 3)
		local phiAng = radToDeg (- 2 * pi * ((double 1)/(double 1.61803399) - 1))
		local theta = 360
		local prevPhi = double 0
		local prevR = double 0

		local positions = for k = 2 to (totalPoints - 1) collect
		(		
			local k2 = a*k + b
			local h = -1 + (double 2) * (k2 - 1)/(totalPoints - 1)
			local r = (1 - h^2)^.5
			local theta = acos h
			local phi = mod (prevPhi + phiAng) 360

			prevPhi = phi
			sphericalToCartesian radius theta phi
		)

		custNode.radius = (distance positions[1] positions[2])/2 * 0.9

		for pos in positions do
			instance custNode pos:pos dir:pos wirecolor:red
		
		delete custNode
	)
)
createDialog fibSphere

Oh I think I like yours better Its nice to see this is well traveled ground. Thanks for sharing.
As far as voronoi partitioning do share as I haven’t gone there yet. I do however want to get to the point where i can build a voronoi based mesh on arbitrary sets of 3d points. not trivial from the looks of it.

Anything else you want to share would be great.

-Michael

Denis would you be willing to take a stab at doing a pure rejection based version of this. I’m sure it would be faster than just relaxing the points but I’m not sure what it would look like. My guess is you need the area of the sphere divided by the number of points to establish a rejection radius then add a new point when it doesn’t overlap any existing points. But I may be missing something. I still have to think its going to get slower prograssivle as more points are added but perhaps it would still cost less than moving all the points to a minimum energy state. god I wish I had taken math seriously as a kid

Thanks

-Michael

Here is a version that uses pure dart throwing, slower than molasses in the dead of winter for large point counts but it works. The distance between the points is calculated based on the area of the underlying sphere and the number of points you want. As you can imagine, the closer you get to that optimum distance the slower it will be to converge cause it takes a ton of attempts to find just the right spots to place a point.

I bet Denis could put this in a grid based acceleration structure and it would be a ton faster. I just need to find a bigger brain.

-Michael

 	----------------------------------------------------------------------------------------
 	--	RandomPointsOnASphere_Using_DartThrowing_05
 	--
 	--	This is a brute force and hence slow way to get N# points on the surface
 	--	of the sphere but its fun to see the points being repulsed by their neighbors
 	--	
 	--	Here are the steps for the algorithm
 	--	1. N# paints are placed randomly within a cube
 	--	2. The points are then projected to the surface of the sphere
 	--	3. Each point is then checked against the others and the closest 
 	--		one is found and moved some amout to even out the distribution
 	--	4. The repulsion most likely moves the point off the surface of the sphere
 	--		so its important to reproject the points back onto the sureface
 	--	5. was and repeat untill either the points are the minimul testing distance away 
 	--	   from one another or we run out of iterations
 	--
 	--	Michael Spaw
 	--	Started 03/06/13
 	--	ver 0.02 03/06/13			  
 	--
 	----------------------------------------------------------------------------------------
 	gc()
 	delete objects--Lets clearout all the objects before we start
 	aqua = (color 87 224 198)						-- Create a color to be used for the bounding shape
 	drawRad = 2.0 										-- The radius of the points
 	nPoints =300											-- Total number of points created within the test
 	Iterations = 300 										-- Total attemps to move the points so they no longer collide
 	TargetSphereRadius = 50.0 						-- Radius of the target sphere
 	BoundingVolume = TargetSphereRadius*2	-- Dimension of the bounding cube the inital points are scattered into
 	halfVol = BoundingVolume*0.5					-- Half of the bounding volume  
 	RepulsionRadius = 2*TargetSphereRadius/(sqrt nPoints)*0.84 -- This is the magic number that sets the repulsion radius based on the number of points and the target sphere
 	RepulsionImpulse = RepulsionRadius*0.015 	-- This is the amout each sphere will be moved during each step of the repulsion
 	RepulsionMultiplier = 0.5 							 --This controls how close to the optimal distance the points will be attempted to be placed 1 is the optimal and slower than snot. 
 																--A value less than one will arrive at the solution fasted but the result will be more points bunced together/Less evenly placed.
 
 	initxpos = 50.0-drawrad --Inital x position for the circle minus some to keep it in the bounding rectangle
 	initypos = 50.0-drawrad --same for the y
 	initzpos = 50.0-drawrad --same for the y
 
 	delete objects--Lets clearout all the objects before we start
 
 	boundingbox= box length:BoundingVolume width:BoundingVolume height:BoundingVolume pos:[0,0,(-halfVol)] wirecolor:aqua  --a nice bounding box for neatness
 	addmodifier boundingbox (Lattice Strut_Radius:0.1 Joint_Radius:0) --Lets make it look like its wirefame
 	SphereProxy = GeoSphere radius:TargetSphereRadius wirecolor:blue segments:8 --This is our standin sphere that the points will be placed on
 	centerPos = SphereProxy.center--set the center of the sphere
 	its =500
 	points= #() 
 	numNodes =0
 	timeStart = timeStamp()
 	
 	--Create a random point in the unit cube and test the vector against the sphere
 	fn getRandomSpherePoint = 
 	(
 		p = (random [-initxpos,-initypos ,-initzpos] [initxpos,initypos,initzpos]) -- Each point is placed within the bounding cube at a random position
 	)	
 
 	while points.count < nPoints do
 	(
 		Col = random green red
 		p = getRandomSpherePoint()
 		(
 		-- Lets make sure that the points are inside the sphere otherwise we 
 		-- will have way too many near the corners of the cube and not enough at the center of the faces	
 		if ((distance centerPos p) < TargetSphereRadius) then 
 			(
 				numToRemove = numNodes
 				for i=1 to its do
 				(
 					L = 9999999.9
 					for j=1 to points.count do
 					(
 						pointVector = p- centerPos --get the vector from the point to the center
 						p= normalize pointVector * halfVol + centerPos --Now move all the points to the surface of the sphere
 						d = distance p points[j]
 						if d < L do L = d
 					)
 					if L > (RepulsionRadius*2*RepulsionMultiplier) do 
 					(
 						GeoSphere pos:p radius:drawRad wirecolor:col segments:3
 						append points p
 						forceCompleteRedraw()
 						windows.processPostedMessages()
 					)
 					if points.count >= (nPoints + numToRemove) do
 					(
 						for i=1 to numToRemove do deleteitem points 1
 					)
 				)
 				for i=1 to numToRemove do deleteitem points 1
 			) 
 		)
 	)
 	
 	timeEnd = timeStamp()
 	forceCompleteRedraw()
 	windows.processPostedMessages()
 	format "Processing took % seconds
" ((timeEnd - timeStart) / 1000.0)			
 	format "Total # of points made % 
" points.count