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

       arrPoints = Rhino.GetObjects("Select Cell Points",1,vbFalse,vbFalse)

       ShowObjects arrBoundingVolume

       If IsNull(arrPoints) Then Exit Sub

       ' go through each point in the list, and create a cell, name it, color it, And update the progress report.

       Dim strCell, strName
       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(255,255,value)
         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

       Next 'i

       ' hide source geometry to reveal the cells
       HideObjects arrPoints
       HideObjects arrBoundingVolume


End Sub



'' ----------------------------------------------------------------------
'' 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
       Dim midpoint, normal

       Dim greatestdiagonalspread
       greatestdiagonalspread = FindGreatestDiagonalSpread(arrPoints)

       Dim centercoords,pointcoords
 centercoords = Rhino.PointCoordinates(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 = Rhino.PointCoordinates(arrPoints(i))
     Dim strPlane
     strPlane = CreateBisectingPlane(centercoords, pointcoords, greatestdiagonalspread)

     If(Rhino.IsSurface(strPlane)) Then 'if it's a valid object, add To the array
             'newPlane has been created

             midpoint = VectorMidpoint(centercoords, pointcoords)

             normal = VectorUnitize(VectorSubtract(centercoords, pointcoords))

             'strPath is now created

             Dim strPath
             strPath = Rhino.AddLine(midpoint, VectorAdd(midpoint,VectorScale(normal, greatestdiagonalspread)))

             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



'' ----------------------------------------------------------------------
'' Takes two point coordinate arrays, finds the midpoint
'' makes a coordinate system based on the line between these two points
'' and uses this coordinate system to make the plane that marks the 3d bisector
'' of that line.
'' IE this plane is the halfway boundary betweensc the two points
'' ----------------------------------------------------------------------

Function CreateBisectingPlane(arrPtOne, arrPtTwo, reach)
 Dim i
 Dim center
 center = VectorMidpoint(arrPtOne, arrPtTwo)

 ' 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 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 )

 CreateBisectingPlane = 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 = Rhino.PointCoordinates(arrPoints(0))
 Dim min : min = Rhino.PointCoordinates(arrPoints(0))
 Dim i, pt, spread

 For i = 0 To UBound(arrPoints)
   pt = Rhino.PointCoordinates(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 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 = Null
 VectorMidpoint = Array((v1(XX) + v2(XX))/2, (v1(YY) + v2(YY))/2,(v1(ZZ) + v2(ZZ))/2)
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