Option Explicit

' 3d Voronoi AKA Project Cell
' (c) Gabe Smedresman 2005
' All Rights Reserved.


''' global naming variables

Dim XX: XX = 0
Dim YY: YY = 1
Dim ZZ: ZZ = 2

GenerateVoronoiCells

Sub GenerateVoronoiCells()
  Randomize

	Dim arrBoundingVolume ' array of polysurfaces indicating bounding volume
	Dim arrPoints ' array of string references to voronoi points

	Dim i,j

  Dim startTime : startTime = Now

	' select objects

  arrBoundingVolume = Rhino.GetObjects("Select a Bounding Volume",16+8,vbTrue,vbTrue)
  If IsNull(arrBoundingVolume) Then Exit Sub
  HideObjects arrBoundingVolume
  
	' can select points and objects
	arrPoints = Rhino.GetObjects("Select Cell Points or Spheres",1+8+16,vbFalse,vbFalse)
	
	ShowObjects arrBoundingVolume
	
	If IsNull(arrPoints) Then Exit Sub

	' Check to make sure that all are points, or all are surfaces.
	Dim arrPoint
	For Each arrPoint in arrPoints
		If ((Rhino.IsSurface(arrPoint) <> Rhino.IsSurface(arrPoints(0))) Or (Rhino.IsPoint(arrPoint) <> Rhino.IsPoint(arrPoints(0)))) Then
			Rhino.Print("The selected objects must be either all points or all volumes.")
			Exit Sub
		End If
	Next
	
	' go through each point in the list, and create a cell, name it, color it, and update the progress report.
	
	Dim strCell, strName
	Dim arrCells()
	ReDim arrCells(UBound(arrPoints))
	
	Rhino.Print "Beginning cell divisions: " & UBound(arrPoints)+1 & " cells total."
	For i = 0 To UBound(arrPoints) ' for each point to make a cell around  ' one point at a time

	  strCell = GenerateCell(arrPoints(i),arrPoints,arrBoundingVolume)
	  Dim value : value = Int(Rnd()*255)
	  Rhino.objectColor strCell, RGB(value,value,255)
	  strName = "Cell #" & i
	  Rhino.objectname strCell, strName

    Dim strTime : strTime = GetTimeDescription(startTime, (i+1) * 1.0 / (UBound(arrPoints)+1) )
    
	  Rhino.Print i+1 & " of " & UBound(arrPoints)+1 & " cells (" & Int((i+1) * 100 / (UBound(arrPoints)+1)) & "%) completed. " & strTime
		arrCells(i) = strCell
	Next 'i
	
	'' Rhino.Print "Checking results."
	'' Dim testSuccess
	'' testSuccess = SanityCheckResults(arrCells, arrBoundingVolume)
	'' If testSuccess Then
	'' 	Rhino.Print("Results OK!")
	'' Else
	'' 	Rhino.Print("Sanity checks failed.")
	'' End If
	
	' hide source geometry to reveal the cells
	HideObjects arrPoints
	HideObjects arrBoundingVolume
	
  
End Sub

'' ----------------------------------------------------------------------
'' Check the results -- make sure that the generated cells when unioned
'' equal the bounding volume, and that none of them intersect each other.
'' ----------------------------------------------------------------------
Function SanityCheckResults(arrCells, arrBoundingVolume)
	Dim strCell, strCell2, arrBooleanResults
	' check that all cells don't intersect any other cell
	For Each strCell in arrCells
		For Each strCell2 in arrCells
			If strCell <> strCell2 Then
				Dim strCopy, strCopy2, centroid1, centroid2, center1, center2, scale, scaleArray
				strCopy = Rhino.CopyObject(strCell)
				strCopy2 = Rhino.CopyObject(strCell2)
				centroid1 = Rhino.SurfaceVolumeCentroid(strCopy)
				center1 = centroid1(0)
				centroid2 = Rhino.SurfaceVolumeCentroid(strCopy2)
				center2 = centroid2(0)
				scale = 0.9999
				scaleArray = Array(scale,scale,scale)
				Rhino.ScaleObject strCopy, center1, scaleArray
				Rhino.ScaleObject strCopy2, center2, scaleArray
				
				arrBooleanResults = Rhino.BooleanIntersection(array(strCopy), array(strCopy2), vbFalse)
				If IsArray(arrBooleanResults) Then
					Rhino.Print("Sanity check failure: Cells " & Rhino.ObjectName(strCell) & " and " & Rhino.ObjectName(strCell2) & " intersect.")
					Rhino.DeleteObjects(arrBooleanResults)
					Rhino.DeleteObject strCopy
					Rhino.DeleteObject strCopy2
					Rhino.objectColor strCell, RGB(255,0,0)
					Rhino.objectColor strCell2, RGB(255,0,0)
					SanityCheckResults = vbFalse
					Exit Function
				End If
				Rhino.DeleteObject strCopy
				Rhino.DeleteObject strCopy2
			End If
		Next
	Next
	SanityCheckResults = vbTrue
End Function

'' ----------------------------------------------------------------------
'' given the center, an array of points, and the bounding volume,
'' perform the operation, make sure all intermediate geometries
'' have been deleted, then return the voronoi cell.
'' ----------------------------------------------------------------------
Function GenerateCell(centerPoint,arrPoints,arrBoundingVolume)
  Dim arrBlocks
	arrBlocks = CreateBlocks(centerPoint, arrPoints)
	
	Dim strCell
  strCell = IntersectBlocks(arrBlocks,arrBoundingVolume)
  
  Dim m: For m = 0 To UBound(arrBlocks)
    If IsPolySurface(arrBlocks(m)) Then DeleteObject(arrBlocks(m))
  Next
    
  GenerateCell = strCell
End Function


'' ----------------------------------------------------------------------
'' using the array of point strings, and one point reference for the center
'' create a group of blocks, which, when intersected, will result in the 
'' voronoi volume. 
''
'' this is done by, for each point, generating a plane perpendicular to the 
'' line between the center and the test point, at the midpoint of that line.
'' this plane faces the center point and is equidistant from the two points
'' across the entire surface.
''
'' then extrude the plane towards the center point, theoretically an infinite
'' amount but really to the end of the point cloud. meaning, this extrusion
'' is closer to the center than it is to the test point.
''
'' intersect all of these volumes and you will have the area closest
'' to the center.
'' ----------------------------------------------------------------------
Function CreateBlocks(centerPoint, arrPoints)

  Dim arrBlocks
  ReDim arrBlocks(UBound(arrPoints))
	
	Dim newPlane
	Dim numBlocks : numBlocks = 0
	Dim i, interpolation
	Dim midpoint, normal
	
	Dim greatestdiagonalspread, geometryExtension
	greatestdiagonalspread = FindGreatestDiagonalSpread(arrPoints)
	geometryExtension = greatestdiagonalspread * 3
	  
	Dim centercoords,pointcoords
	centercoords = GetPointCoordinates(centerPoint)
  
	' ---- CREATE BOUNDING PLANES
	For i = 0 To UBound(arrPoints) ' for each point to consider boundary to

    If Not centerPoint = arrPoints(i) Then ' as long as the point isn't the central point
  
      pointcoords = GetPointCoordinates(arrPoints(i))
      
      interpolation = GetRadicalPlanePosition(centerPoint, arrPoints(i))
			midpoint = VectorInterpolate(centercoords, pointcoords, interpolation)

      ' normal goes from point inwards towards center
      normal = VectorUnitize(VectorSubtract(centercoords, pointcoords))
      
      Dim strPlane
      strPlane = CreatePlane(midpoint, normal, geometryExtension)
      
      If(Rhino.IsSurface(strPlane)) Then 'if it's a valid object, add to the array
	      'newPlane has been created

	      'strPath is now created	      
	      Dim strPath
	      strPath = Rhino.AddLine(midpoint, VectorAdd(midpoint,VectorScale(normal, geometryExtension)))
	      
	      Dim strExtrusion
	      strExtrusion = Rhino.ExtrudeSurface( strPlane, strPath )
	      Rhino.objectColor strExtrusion, RGB(128,128,128)
	      
        arrBlocks(numBlocks) = strExtrusion
	      numBlocks = numBlocks + 1
	      
	      If(Rhino.IsObject(strPath)) Then Rhino.DeleteObject(strPath)
	      If(Rhino.IsObject(strPlane)) Then Rhino.DeleteObject(strPlane)
	    End If 'end if it's a surface
	  End If 'end if considering different points
	Next 'i
	
	ReDim Preserve arrBlocks(numBlocks-1)
	'Rhino.Print(numBlocks & " blocks created.")
	
	CreateBlocks = arrBlocks
End Function

'' ----------------------------------------------------------------------
'' If the two objects are points, use the midpoint.
'' Otherwise, they are spheres.
'' If the spheres are equal size, return 0.5.
'' 0.0 would be at the center
'' sphere, and 1.0 would be at the distant sphere.
'' ----------------------------------------------------------------------
Function GetRadicalPlanePosition(strCenterSphere, strDistantSphere)
	Dim p
	If Rhino.IsPoint(strCenterSphere) Or Rhino.IsPoint(strDistantSphere) Then
		p = 0.5
	Else
		p = GetConstraint(strCenterSphere, strDistantSphere)
	End If
	GetRadicalPlanePosition = p
End Function

'' ----------------------------------------------------------------------
'' Given two spheres, find the radical plane position
'' between them, as a number between 0 (on center sphere)
'' and 1 (on distant sphere)
'' ----------------------------------------------------------------------
Function GetConstraint(strCenterSphere, strDistantSphere)
	Dim Ri, Rj, p
	Dim centerPoint, distantPoint, Dij

	
	Ri = GetObjectRadius(strCenterSphere)
	Rj = GetObjectRadius(strDistantSphere)
		
	centerPoint = GetPointCoordinates(strCenterSphere)
	distantPoint = GetPointCoordinates(strDistantSphere)
	Dij = VectorLength(VectorSubtract(centerPoint, distantPoint))
	
  ' first  calculate the raw distance
	p = ((Dij * Dij) + (Ri * Ri) - (Rj * Rj)) / (2 * Dij)
  ' then scale to 0.0-1.0
  p = p / Dij
  GetConstraint = p
End Function

'' ----------------------------------------------------------------------
'' Given a rhino object, in this case, either a point or a sphere, return
'' the coordinates of the center of that object.
'' ----------------------------------------------------------------------
Function GetPointCoordinates(strObject)
	Dim coords
	If Rhino.IsPoint(strObject) Then
		coords = Rhino.PointCoordinates(strObject)
	ElseIf Rhino.IsSurface(strObject) Then
		Dim arrBoundingBox
		arrBoundingBox = Rhino.BoundingBox(strObject)
		If IsArray(arrBoundingBox) Then
			coords = VectorMidpoint(arrBoundingBox(0), arrBoundingBox(6))
		Else
			Rhino.Print("Error: couldn't get the bounding box of the object.")
			coords = vbNull
		End If
	Else
		Rhino.Print("Error: tried to get the coordinates of an invalid object")
		coords = vbNull
	End If
	GetPointCoordinates = coords
End Function

'' ----------------------------------------------------------------------
'' Detect the 'radius' of any rhino object.  For spheres this should
'' return the radius of the sphere, but for the sake of robustness, and
'' because it's very difficult to determine what type of object an
'' arbitrary surface is, just return one half the max dimension of the
'' bounding box of the object.
'' ----------------------------------------------------------------------
Function GetObjectRadius(strObject)
	Dim radius
	
	If Rhino.IsPoint(strObject) Then
		Rhino.Print("Error: tried to get the radius of a point")
		radius = 0
	ElseIf Rhino.IsSurface(strObject) Then
		Dim arrBoundingBox, arrSpan
		arrBoundingBox = Rhino.BoundingBox(strObject)
		If IsArray(arrBoundingBox) Then
			arrSpan = VectorSubtract(arrBoundingBox(6), arrBoundingBox(0))
			radius = VectorMaxDimension(arrSpan) * 0.5
		Else
			Rhino.Print("Error: couldn't get the bounding box of the object.")
			radius = 0
		End If
	Else
		Rhino.Print("Error: tried to get the radius of a non-surface object")
		radius = 0
	End If
	GetObjectRadius = radius
End Function


'' ----------------------------------------------------------------------
'' Constructs a plane centered at the given point, facing
' the given normal, with a width and height of 'reach',
'' ----------------------------------------------------------------------
Function CreatePlane(point, normal, reach)
  Dim i
  Dim center: center = point
  
  ' make new coordinate system p,r,s: p is the line between 1 and 2, r is to the side, s is up(ish) from the line
  ' Dim p :  p  = VectorSubtract(arrPtTwo,arrPtOne) : p = VectorUnitize(p)
  Dim p : p = normal
  
  
  Dim up : up = Array(0,0,1)
  Dim r : r  = VectorCrossProduct(p,up) : If IsVectorZero(r) Then r = Array(0,1,0)
  r = VectorUnitize(r) '' points to the right
  Dim s :  s  = VectorCrossProduct(p,r) : s = VectorUnitize(s) '' points to perpendicular to p (forward) and r (side)
  
  ' now find four points, 1 up 2 left 3 down 4 right (looking from one to two)
  Dim arrCorners(3)
  arrCorners(0) = VectorAdd(center,VectorScale(s,reach))
  arrCorners(1) = VectorAdd(center,VectorScale(r,reach * -1))
  arrCorners(2) = VectorAdd(center,VectorScale(s,reach * -1))
  arrCorners(3) = VectorAdd(center,VectorScale(r,reach))
  
  Dim strPlane
  strPlane = Rhino.AddSrfPt( arrCorners )

  CreatePlane = strPlane
End Function

'' ----------------------------------------------------------------------
'' given the array of generated blocks, perform a boolean intersection
'' with the bounding volume and each block, one by one.
'' delete each source as you iterate. some blocks remain: if
'' the intersection fails (no areas intersect), the sources 
'' will not be deleted.
'' ----------------------------------------------------------------------
Function IntersectBlocks(arrBlocks, arrBoundingVolume)
  Dim i
  
  IntersectBlocks = Null
  i = 0
  Dim results, pendingresults
  Dim strBlock
  Dim group2

  Do
    strBlock = arrBlocks(i)
    group2 = array(strBlock)
    pendingresults = Rhino.BooleanIntersection(arrBoundingVolume,group2,vbFalse)
    i = i + 1
  Loop While Not IsArray(pendingresults)
  
  results = pendingresults
    
  Dim j
  For j = i To UBound(arrBlocks)
    strBlock = arrBlocks(j)
    group2 = array(strBlock)
    pendingresults = Rhino.BooleanIntersection(results,group2)
    If IsArray(pendingresults) Then results = pendingresults

  Next  
  
  IntersectBlocks = results(0)
End Function


'' ----------------------------------------------------------------------
'' find the maximum 3d diagonal length between one extreme corner
'' of a point cloud and the other
'' ----------------------------------------------------------------------
Function FindGreatestDiagonalSpread(arrPoints)
  If(Not IsArray(arrPoints)) Then FindGreatestDiagonalSpread = 0
  
  Dim max : max = GetPointCoordinates(arrPoints(0))
  Dim min : min = GetPointCoordinates(arrPoints(0))
  Dim i, pt, spread
  
  For i = 0 To UBound(arrPoints)
    pt = GetPointCoordinates(arrPoints(i))
    If(max(XX) < pt(XX)) Then max(XX) = pt(XX)
    If(max(YY) < pt(YY)) Then max(YY) = pt(YY)
    If(max(ZZ) < pt(ZZ)) Then max(ZZ) = pt(ZZ)
    If(min(XX) > pt(XX)) Then min(XX) = pt(XX)
    If(min(YY) > pt(YY)) Then min(YY) = pt(YY)
    If(min(ZZ) > pt(ZZ)) Then min(ZZ) = pt(ZZ)
  Next 'i
 
  FindGreatestDiagonalSpread = VectorLength(VectorSubtract(max,min))
End Function

'' ----------------------------------------------------------------------
'' given the starting time and fraction complete, estimate time to
'' completion and then generate a sentence of the format
'' "3 minutes 4 seconds elapsed: should be finished at 04:00 PM today."
'' ----------------------------------------------------------------------
Function GetTimeDescription(startTime, fractioncomplete)
  	Dim strDescription
  	strDescription = ""
	  Dim elapsedseconds : elapsedseconds = DateDiff("s",startTime,Now)
	  Dim m,h,s : s = elapsedseconds
	  m = Int(s/60) : s = s - m * 60
	  h = Int(m/60) : m = m - h * 60
	  If h = 1 Then strDescription = strDescription & " " & h & " hour"
	  If h > 1 Then strDescription = strDescription & " " & h & " hours"
	  If m = 1 Then strDescription = strDescription & " " & m & " minute"
	  If m > 1 Then strDescription = strDescription & " " & m & " minutes"
	  If s = 1 Then strDescription = strDescription & " " & s & " second"
	  If s > 1 Then strDescription = strDescription & " " & s & " seconds"
	  strDescription = strDescription & " elapsed:"
	  
	  Dim secondstogo : secondstogo = elapsedseconds / fractioncomplete * (1-fractioncomplete)
	  Dim ETA : ETA = DateAdd("s",secondstogo,Now)
	  
	  If(fractioncomplete < 1) Then
	  
	    strDescription = strDescription & " should be finished at " & FormatDateTime(ETA,4)	  
	    If(DatePart("d",ETA) = DatePart("d",Now)) Then
  	    strDescription = strDescription & " today."
	    Else
  	    strDescription = strDescription & " in " & DateDiff("d",Now,ETA) & " day(s)."
	    End If
	  Else
	    strDescription = strDescription & " finished " & Now & "!"
	  End If
	  
	  GetTimeDescription = strDescription
	  
End Function




''' ***************************************************************************
''' Rhinoscript Vector Functions
''' ***************************************************************************

''' ---------------------------------------------------------------------------
''' Make a vector from two 3D points
''' ---------------------------------------------------------------------------
Public Function VectorCreate(p1, p2)
  VectorCreate = Array(p2(0) - p1(0), p2(1) - p1(1), p2(2) - p1(2))
End Function

''' ---------------------------------------------------------------------------
''' Unitize a 3D vector
''' ---------------------------------------------------------------------------
Public Function VectorUnitize(v)

  VectorUnitize = Null
  Dim dist, x, y, z, x2, y2, z2

  x = v(XX) : y = v(YY) : z = v(ZZ)
  x2 = x * x : y2 = y * y : z2 = z * z

  dist = x2 + y2 + z2
  If dist <= 0.0 Then Exit Function
  dist = Sqr(dist)

  x = x / dist
  y = y / dist
  z = z / dist

  VectorUnitize = Array(x, y, z)
End Function

''' ---------------------------------------------------------------------------
''' Return the length of a 3D vector
''' ---------------------------------------------------------------------------
Public Function VectorLength(v)

  VectorLength = Null

  Dim dist, x, y, z, x2, y2, z2
  x = v(XX) : y = v(YY) : z = v(ZZ)
  x2 = x * x : y2 = y * y : z2 = z * z
  dist = x2 + y2 + z2

  VectorLength = Sqr(dist)
  
End Function

''' ---------------------------------------------------------------------------
''' Return the absolute value of the maximum dimension of this vector
''' ---------------------------------------------------------------------------
Public Function VectorMaxDimension(v)
	Dim maxDim
	
	If Abs(v(XX)) > Abs(v(YY)) Then		' X > Y
		If Abs(v(XX)) > Abs(v(ZZ)) Then	' X > Y and X > Z
			maxDim = Abs(v(XX))
		Else														' X > Y and Z >= X
			maxDim = Abs(v(ZZ))
		End If
	Else															' Y > X
		If Abs(v(YY)) > Abs(v(ZZ)) Then	' Y > X and Y > Z
			maxDim = Abs(v(YY))
		Else														' Y > X and Z >= Y
			maxDim = Abs(v(ZZ))
		End If
	End If
	VectorMaxDimension = maxDim
End Function

''' ---------------------------------------------------------------------------
''' Return the dot product of two 3D vectors
''' ---------------------------------------------------------------------------
Public Function VectorDotProduct(v1, v2)
  VectorDotProduct = v1(XX) * v2(XX) + v1(YY) * v2(YY) + v1(ZZ) * v2(ZZ)
End Function

''' ---------------------------------------------------------------------------
''' Return the cross product of two 3D vectors
''' ---------------------------------------------------------------------------
Public Function VectorCrossProduct(v1, v2)

  VectorCrossProduct = Null
  Dim x, y, z
  x = v1(YY) * v2(ZZ) - v1(ZZ) * v2(YY)
  y = v1(ZZ) * v2(XX) - v1(XX) * v2(ZZ)
  z = v1(XX) * v2(YY) - v1(YY) * v2(XX)
  VectorCrossProduct = Array(x, y, z)

End Function

 
''' ---------------------------------------------------------------------------
''' Add two 3D vectors
''' ---------------------------------------------------------------------------
Public Function VectorAdd(v1, v2)
  VectorAdd = Null
  VectorAdd = Array(v1(XX) + v2(XX), v1(YY) + v2(YY), v1(ZZ) + v2(ZZ))
End Function

''' ---------------------------------------------------------------------------
''' Subtract two 3D vectors
''' ---------------------------------------------------------------------------
Public Function VectorSubtract(v1, v2)
  VectorSubtract = Null
  VectorSubtract = Array(v1(XX) - v2(XX), v1(YY) - v2(YY), v1(ZZ) - v2(ZZ))
End Function

''' ---------------------------------------------------------------------------
''' Multiply two 3D vectors
''' ---------------------------------------------------------------------------
Public Function VectorMultiply(v1, v2)
  VectorMultiply = Null
  VectorMultiply = Array(v1(XX) * v2(XX), v1(YY) * v2(YY), v1(ZZ) * v2(ZZ))
End Function

''' ---------------------------------------------------------------------------
''' Scale a 3D vectors by a value
''' ---------------------------------------------------------------------------
Public Function VectorScale(v, d)
  VectorScale = Null
  VectorScale = Array(v(XX) * d, v(YY) * d, v(ZZ) * d)
End Function

''' ---------------------------------------------------------------------------
''' Compare two 3D vectors for equality
''' ---------------------------------------------------------------------------

Public Function VectorCompare(v1, v2)
  VectorCompare = vbFalse

  If v1(XX) = v2(XX) And v1(YY) = v2(YY) And v1(ZZ) = v2(ZZ) Then
    VectorCompare = vbTrue
  End If

End Function
 
''' ---------------------------------------------------------------------------
''' Negate a 3D vector
''' ---------------------------------------------------------------------------
Public Function VectorNegate(v)
  VectorNegate = Null
  VectorNegate = Array(-v(XX), -v(YY), -v(ZZ))
End Function

''' ---------------------------------------------------------------------------
''' Tiny vector test
''' ---------------------------------------------------------------------------

Public Function IsVectorTiny(v)

  IsVectorTiny = vbFalse
  Dim tol :  tol = 1.0e-12 ' ON_ZERO_TOLERANCE

  If (Abs(v(XX)) <= tol) And (Abs(v(YY)) <= tol) And (Abs(v(ZZ)) <= tol) Then
    IsVectorTiny = vbTrue
  End If

End Function

''' ---------------------------------------------------------------------------
''' Zero vector test
''' ---------------------------------------------------------------------------
Public Function IsVectorZero(v)

  IsVectorZero = vbFalse
  If (v(XX) = 0.0) And (v(YY) = 0.0) And (v(ZZ) = 0.0) Then IsVectorZero = vbTrue

End Function

''' My more specialized functions

''' ---------------------------------------------------------------------------
''' Find midpoint between two vectors
''' ---------------------------------------------------------------------------
Public Function VectorMidpoint(v1, v2)
  VectorMidpoint = Array((v1(XX) + v2(XX))/2, (v1(YY) + v2(YY))/2, (v1(ZZ) + v2(ZZ))/2)
End Function

''' ---------------------------------------------------------------------------
''' Interpolate between two vectors. If p = 0, return v1. If p = 1, return v2.
''' If p = 0.5, return the midpoint, and so on linearly between the two points.
''' ---------------------------------------------------------------------------
Public Function VectorInterpolate(v1, v2, p)
  VectorInterpolate = Array((v1(XX)*(1-p) + v2(XX)*p), (v1(YY)*(1-p) + v2(YY)*p), (v1(ZZ)*(1-p) + v2(ZZ)*p))
End Function

''' ---------------------------------------------------------------------------
''' Return the length of a 3D vector
''' ---------------------------------------------------------------------------
Public Function VectorSquaredDistance(v1,v2)

  VectorSquaredDistance = Null

  Dim dist, x, y, z, x2, y2, z2
  x = v1(XX) - v2(XX)
  y = v1(YY) - v2(YY)
  z = v1(ZZ) - v2(ZZ)
  x2 = x * x
  y2 = y * y
  z2 = z * z
  dist = x2 + y2 + z2

  VectorSquaredDistance = dist
  
End Function