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 def computeGradient(v, index,MS,x, normalPerFaces): grad=rg.Matrix(3,1) 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) return grad 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): grad=computeGradient(v,index, MS,x, normalPerFaces) hess=computeHessian(v,index,MS, x, normalPerFaces) hess.Invert(0.0) dx=(hess*grad) 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] if (math.sqrt(grad[0,0]*grad[0,0]+grad[1,0]*grad[1,0]+grad[2,0]*grad[2,0])< tol): 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] if f.IsQuad: _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) self.faceTableAdd(_f) self.halfEdgeAdd(self.faces[0]) 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): v.star.Add(h) 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) def halfEdgeAdd(self,f): 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): self.halfEdgeAdd(faces[1]); if (faces[1] is f): if (self._orientation[faces[0].N] == enum._UNKNOWN): self.halfEdgeAdd(faces[0]) else: faces=self.getfaces(J, I) if (faces.Count != 0): if (self._orientation[faces[0].N] == enum._UNKNOWN): self.halfEdgeAdd(faces[0]) def faceTableAdd2(self,i,j,f): self._faceTable[i].append((f,j)) def faceTableAdd(self,f): 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] self.faceTableAdd2(I, J, f); @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): newMesh.Vertices.Add(rg.Point3d(0,0,0)) 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, f.B, f.A) newMesh.Faces.AddFace(f.C + myMesh.Vertices.Count, f.B + myMesh.Vertices.Count, f.A + myMesh.Vertices.Count) else: newMesh.Faces.AddFace(f.D,f.C, f.B, f.A) 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 newMesh.Faces.AddFace(rg.MeshFace(d, c, b, a)) h = h.next if(h is boundary): break newMesh.UnifyNormals() a=newMesh
page revision: 10, last edited: 14 Oct 2015 16:51
