Thickness

This small script is what I wrote as the prototype for the sprout.
This script contains class definitions of half-edge data structure (which itself would be useful).

• GH file containing this Python script.

http://mikity.wikidot.com/local--files/detail/Thickeness.gh

• An example mesh surface (containing a hole). Require Rhinoceros.

http://mikity.wikidot.com/local--files/detail/testMesh.3dm

```import Rhino.Geometry as rg
import rhinoscriptsyntax as rs
import math
for e in v.onering:
N = normalPerFaces[e.owner.N]
cx = x[index * 3 + 0]
cy = x[index * 3 + 1]
cz = x[index * 3 + 2]

grad[0,0] = grad[0,0] + 2 * N.X * N.X * cx + 2 * N.X * N.Y * cy + 2 * N.X * N.Z * cz - 2 * N.X
grad[1,0] = grad[1,0] + 2 * N.Y * N.Y * cy + 2 * N.Y * N.Z * cz + 2 * N.Y * N.X * cx - 2 * N.Y
grad[2,0] = grad[2,0] + 2 * N.Z * N.Z * cz + 2 * N.Z * N.X * cx + 2 * N.Z * N.Y * cy - 2 * N.Z

dot = cx * N.X + cy * N.Y + cz * N.Z
norm = N.X*N.X+N.Y*N.Y+N.Z*N.Z

grad[0,0] = grad[0,0] + 0.01*(2 * cx - 4 * dot * N.X + norm * 2 * dot * N.X)
grad[1,0] = grad[1,0] + 0.01*(2 * cy - 4 * dot * N.Y + norm * 2 * dot * N.Y)
grad[2,0] = grad[2,0] + 0.01*(2 * cz - 4 * dot * N.Z + norm * 2 * dot * N.Z)

def computeHessian(v,index,MS,x,normalPerFaces):
hess=rg.Matrix(3,3)
for i in range(3):
for j in range(3):
hess[i,j]=0
for e in v.onering:
N = normalPerFaces[e.owner.N]
cx = x[index * 3 + 0]
cy = x[index * 3 + 1]
cz = x[index * 3 + 2]

hess[0,0]= hess[0,0]+2 * N.X * N.X
hess[1,0]= hess[1,0]+2 * N.X * N.Y
hess[2,0]= hess[2,0]+2 * N.X * N.Z
hess[0,1]= hess[0,1]+2 * N.Y * N.X
hess[1,1]= hess[1,1]+2 * N.Y * N.Y
hess[2,1]= hess[2,1]+2 * N.Y * N.Z
hess[0,2]= hess[0,2]+2 * N.Z * N.X
hess[1,2]= hess[1,2]+2 * N.Z * N.Y
hess[2,2]= hess[2,2]+2 * N.Z * N.Z

dot = cx * N.X + cy * N.Y + cz * N.Z
norm = N.X*N.X+N.Y*N.Y+N.Z*N.Z

hess[0,0] = hess[0,0]+0.01*(2 - 4 * N.X * N.X + norm * 2 * N.X * N.X)
hess[1,0] = hess[1,0]+0.01*(-4 * N.Y * N.X + norm * 2 * N.Y * N.X)
hess[2,0] = hess[2,0]+0.01*(-4 * N.Z * N.X + norm * 2 * N.Z * N.X)
hess[1,1] = hess[1,1]+0.01*(2 - 4 * N.Y * N.Y + norm * 2 * N.Y * N.Y)
hess[2,1] = hess[2,1]+0.01*(-4 * N.Z * N.Y + norm * 2 * N.Z * N.Y)
hess[0,1] = hess[0,1]+0.01*(-4 * N.X * N.Y + norm * 2 * N.X * N.Y)
hess[2,2] = hess[2,2]+0.01*(2 - 4 * N.Z * N.Z + norm * 2 * N.Z * N.Z)
hess[0,2] = hess[0,2]+0.01*(-4 * N.Y * N.Z + norm * 2 * N.X * N.Z)
hess[1,2] = hess[1,2]+0.01*(-4 * N.Y * N.Z + norm * 2 * N.Y * N.Z)

return hess
def computeNormal(input):
(t,myMesh,MS)=input
numvar = len(MS.vertices) * 3
x = range(numvar) #initialize
nodes = myMesh.Vertices
normalPerFaces = range(MS.nFaces())
#initial guess
for (v,index) in zip(MS.vertices,range(len(MS.vertices))):
N=rg.Vector3d(0,0,0)
for e in v.onering:
#compute normal
P = nodes[e.P.N]
Q = nodes[e.next.P.N]
R = nodes[e.next.next.P.N]
b = P-Q
c = R-Q
n = rg.Vector3d.CrossProduct(b,c)
n.Unitize()
normalPerFaces[e.owner.N] = n
N = N+n
N.Unitize()
x[index * 3 + 0] = N.X
x[index * 3 + 1] = N.Y
x[index * 3 + 2] = N.Z
tol=0.00000001
index=0
for (v,index) in zip(MS.vertices,range(len(MS.vertices))):
for i in range(1):
hess=computeHessian(v,index,MS, x, normalPerFaces)
hess.Invert(0.0)
x[index * 3 + 0] = x[index*3+0]-dx[0,0]
x[index * 3 + 1] = x[index*3+1]-dx[1,0]
x[index * 3 + 2] = x[index*3+2]-dx[2,0]

break
output={}
for (v,index) in zip(MS.vertices,range(len(MS.vertices))):
P=nodes[v.N]
N=rg.Vector3f(x[index * 3 + 0],x[index * 3 + 1],x[index * 3 + 2])
output[v]=(P,N)
return (t,MS,output)

class _enum:
_NULL=-100
_POSITIVE=0
_UNKNOWN=-1
_NEGATIVE=2
enum=_enum()
class face:
def __init__(self,_N, _corner):
self.corner=_corner
self.N = _N;
class halfedge:
def isNaked(self):
global enum
if self.pair == enum._NULL:
return True
else:
return False

def isBoundary(self):
global enum
if self.owner == enum._NULL:
return True
else:
return False
def __init__(self,_P):
global enum
self.P = _P;
self.pair=enum._NULL
if _P.hf_begin == enum._NULL:
_P.hf_begin = self

class vertex:
star = []
onering = []

def __init__(self,_N):
global enum
self.hf_begin=enum._NULL
self.N = _N;
def isInner(self):
return self.onering[0] == self.onering[len(onering) - 1].next
def isBoundary(self):
return self.onering[0] != self.onering[len(onering) - 1].next;
class MeshStructure:
boundaryStart=[]
vertices = []
faces = []
halfedges = []
innerVertices = []
outerVertices = []
_halfedgeTable = []
_faceTable=[] #items are also lists
def edges(self):
res = []
for e in self.halfedges:
if (e.P.N < e.next.P.N):
res.append(e)
return res
def nVertices(self):
return len(self.vertices)
def nFaces(self):
return len(self.faces)
def __init__(self):
self.Clear()

def getfaces(self,I,J):
faces=[]
ff = self._faceTable[I]
for (_f,_J) in ff:
if _J==J:
faces.append(_f)
return faces
def gethalfedges(self,I, J):
lhalfedges=[]
ff = self._halfedgeTable[I]
for f in ff:
if f.next.P.N == J:
lhalfedges.append(f)
return lhalfedges

def Construct(self,val):
global enum
_nVertices = val.Vertices.Count
_nFaces = val.Faces.Count
self._orientation =[]
for i in range(_nFaces):
self._orientation.append(enum._UNKNOWN)
self._faceTable = []
self._halfedgeTable = []
for i in range(_nVertices):
self._faceTable.append([])
self._halfedgeTable.append([])

for i in range(_nVertices):
_v = vertex(i)
self.vertices.append(_v)

for i in range(_nFaces):
f = val.Faces[i]
_f=face(i,(f.A,f.B,f.C,f.D))
else:
_f = face(i, (f.A, f.B, f.C))
self.faces.append(_f)

for h in self.halfedges:
i = h.P.N
j = h.next.P.N;
lhalfedges=self.gethalfedges(i, j)
self._halfedgeTable[i].append(h)

for h in self.halfedges:
i = h.P.N
j = h.next.P.N
lhalfedges=self.gethalfedges(j, i)
if (len(lhalfedges)==0):
h.pair = enum._NULL
else:
h.pair = lhalfedges[0]

for v in self.vertices:
h = v.hf_begin
while(True):
if (h.prev.isNaked()):
while (not h.isNaked()):
h = h.pair.next
v.hf_begin = h
break
h = h.prev.pair
if h is v.hf_begin:
break

boundary_complements = []
nBoundary = 0
for v in self.vertices:
h = v.hf_begin
if (h.isNaked()):
flag = True
for i in range(nBoundary):
if h in boundary_complements[i]:
flag = False
break
if flag:
boundary_complements.append([]);
while(True):
boundary_complements[nBoundary].append(h)
h = h.next.P.hf_begin
if h is v.hf_begin:
break
nBoundary=nBoundary+1

self.boundaryStart = []
for i in range(nBoundary):
self.boundaryStart.append(enum._NULL)
for boundary_complement in boundary_complements:
boundary = []
for i in range(len(boundary_complement)):
boundary.append(halfedge(boundary_complement[i].next.P))
boundary[i].pair = boundary_complement[i]
boundary_complement[i].pair = boundary[i]
self.halfedges.extend(boundary)
self.boundaryStart[boundary_complements.IndexOf(boundary_complement)] = boundary[0]
for  i in range(len(boundary)):
boundary[i].owner = enum._NULL
if (i != 0):
boundary[i].next = boundary_complement[i - 1].pair
else:
boundary[i].next = boundary_complement[len(boundary_complement) - 1].pair
if (i != boundary.Count - 1):
boundary[i].prev = boundary_complement[i + 1].pair
else:
boundary[i].prev = boundary_complement[0].pair
for e in self.halfedges:
if (e.isNaked()):
print("error")
for v in self.vertices:
h = v.hf_begin
v.star=[]
while(True):
if (h.isBoundary()):
break
h = h.prev.pair
if h is v.hf_begin:
break
for v in self.vertices:
h = v.hf_begin
v.onering=[]
while(True):
while(True):
h = h.next
v.onering.append(h)
if h.next.next.P is v:
break
if h.next.pair.isBoundary():
break
h = h.next.pair
if h is v.hf_begin:
break
self.innerVertices=[]
self.outerVertices=[]
for v in self.vertices:
if v.hf_begin.pair.isBoundary():
self.outerVertices.append(v)
else:
self.innerVertices.append(v)

global enum
_o = enum._UNKNOWN
for i in range(len(f.corner)):
I = f.corner[i]
if (i == len(f.corner) - 1):
J=f.corner[0]
else:
J=f.corner[i + 1]
faces=self.getfaces(I, J)
if (len(faces) == 2):
if (faces[0] is f):
if (_orientation[faces[1].N] != enum._UNKNOWN):
if _orientation[faces[1].N] == enum._POSITIVE:
_o=enum._NEGATIVE
else:
_o=enum._POSITIVE
if (faces[1] is f):
if (_orientation[faces[0].N] != enum._UNKNOWN):
if _orientation[faces[0].N] == enum._POSITIVE:
_o=enum._NEGATIVE
else:
_o=enum._POSITIVE
else:
faces=self.getfaces(J, I)
if len(faces) != 0:
if (self._orientation[faces[0].N] != enum._UNKNOWN):
_o = self._orientation[faces[0].N]
if _o==enum._UNKNOWN:
self._orientation[f.N] = enum._POSITIVE
else:
self._orientation[f.N] = _o
if (len(f.corner) == 3 and self._orientation[f.N] == enum._POSITIVE):
he1 = halfedge(self.vertices[f.corner[0]])
he2 = halfedge(self.vertices[f.corner[1]])
he3 = halfedge(self.vertices[f.corner[2]])
self.halfedges.append(he1)
self.halfedges.append(he2)
self.halfedges.append(he3)
he1.prev = he3
he1.next = he2
he1.owner = f
he2.prev = he1
he2.next = he3
he2.owner = f
he3.prev = he2
he3.next = he1
he3.owner = f
f.firsthalfedge = he1

if (len(f.corner) == 3 and self._orientation[f.N] == enum._NEGATIVE):
he1 = halfedge(self.vertices[f.corner[2]])
he2 = halfedge(self.vertices[f.corner[1]])
he3 = halfedge(self.vertices[f.corner[0]])
self.halfedges.append(he1)
self.halfedges.append(he2)
self.halfedges.append(he3)
he1.prev = he3
he1.next = he2
he1.owner = f
he2.prev = he1
he2.next = he3
he2.owner = f
he3.prev = he2
he3.next = he1
he3.owner = f
f.firsthalfedge = he1
if (len(f.corner) == 4 and self._orientation[f.N] == enum._POSITIVE):
he1 = halfedge(self.vertices[f.corner[0]])
he2 = halfedge(self.vertices[f.corner[1]])
he3 = halfedge(self.vertices[f.corner[2]])
he4 = halfedge(self.vertices[f.corner[3]])
self.halfedges.append(he1)
self.halfedges.append(he2)
self.halfedges.append(he3)
self.halfedges.append(he4)
he1.prev = he4
he1.next = he2
he1.owner = f
he2.prev = he1
he2.next = he3
he2.owner = f
he3.prev = he2
he3.next = he4
he3.owner = f
he4.prev = he3
he4.next = he1
he4.owner = f
f.firsthalfedge = he1
if (len(f.corner) == 4 and self._orientation[f.N] == enum._NEGATIVE):
he1 = halfedge(self.vertices[f.corner[3]])
he2 = halfedge(self.vertices[f.corner[2]])
he3 = halfedge(self.vertices[f.corner[1]])
he4 = halfedge(self.vertices[f.corner[0]])
self.halfedges.append(he1)
self.halfedges.append(he2)
self.halfedges.append(he3)
self.halfedges.append(he4)
he1.prev = he4
he1.next = he2
he1.owner = f
he2.prev = he1
he2.next = he3
he2.owner = f
he3.prev = he2
he3.next = he4
he3.owner = f
he4.prev = he3
he4.next = he1
he4.owner = f
f.firsthalfedge = he1
for i in range(len(f.corner)):
I = f.corner[i]
if(i == len(f.corner) - 1):
J=f.corner[0]
else:
J=f.corner[i + 1]
faces=self.getfaces(I, J)
if (len(faces) == 2):
if (faces[0] is f):
if (self._orientation[faces[1].N] == enum._UNKNOWN):
if (faces[1] is f):
if (self._orientation[faces[0].N] == enum._UNKNOWN):
else:
faces=self.getfaces(J, I)
if (faces.Count != 0):
if (self._orientation[faces[0].N] == enum._UNKNOWN):
self._faceTable[i].append((f,j))
for i in range(len(f.corner)):
I = f.corner[i]
if (i == len(f.corner)-1):
J = f.corner[0]
else:
J = f.corner[i + 1]
@staticmethod
def CreateFrom(val):
ret = MeshStructure()
ret.Construct(val)
return ret
def Clear(self):
global enum
self.vertices=[]
self.faces=[]
self.halfedges=[]
self.innerVertices=[]
self.outerVertices=[]
self.boundaryStart = enum._NULL

myMesh=M
thickness = t
MS = MeshStructure.CreateFrom(myMesh)
(t,MS,output)=computeNormal((t,myMesh,MS))

listNormal = []
newMesh=rg.Mesh()
for i in range(myMesh.Vertices.Count*2):

for v in MS.vertices:
(P,N)=output[v]
P1=P+rg.Vector3f.Multiply(thickness/2.0,N)
P2=P+rg.Vector3f.Multiply(thickness/2.0,rg.Vector3f(-N.X,-N.Y,-N.Z))
newMesh.Vertices[v.N]=P1
listNormal.append(rg.Line(P1,P2))
newMesh.Vertices[myMesh.Vertices.Count+v.N]=P2

for i in range(myMesh.Faces.Count):
f = myMesh.Faces[i]
if (f.IsTriangle):
newMesh.Faces.AddFace(f.C + myMesh.Vertices.Count, f.B + myMesh.Vertices.Count, f.A + myMesh.Vertices.Count)
else:
newMesh.Faces.AddFace(f.D+myMesh.Vertices.Count,f.C + myMesh.Vertices.Count, f.B + myMesh.Vertices.Count, f.A + myMesh.Vertices.Count)

for boundary in MS.boundaryStart:
h = boundary
while(True):
a = h.P.N
b = h.next.P.N
c = b + myMesh.Vertices.Count
d = a + myMesh.Vertices.Count