$search

tesselatedSphere.py

Go to the documentation of this file.

from Triangle import Triangle
#from Numeric import *
from numpy import *
from numpy.linalg import *
import vtk
import time
import os.path
from Cell import Cell

class TesselatedSphere:
        def __init__(self, levels, inter):
                self.triangles = list()
                self.icosahedron = list()
                self.polygons = list()
                self.max_levels = levels
                self.resolution = 0
                self.interpolateNormals = inter
                self.useLut = False
                self.convexHullPoints = dict()
                self.useNormalsOnSurface = True
                
        def resetPolygons(self):
                for cell in self.polygons[self.max_levels]:
                        cell.reset()
                
        def construct(self):
                
                #if file with lookup table exist, read it from file together with the
                #corresponding cells
                
                #create icosahedron vertices
                #define regular icosahedron vertex coordinates
                #http://en.wikipedia.org/wiki/Icosahedron

                tau = 0.8506508084 # t=(1+sqrt(5))/2, tau=t/sqrt(1+t^2) 
                one = 0.5257311121 # one=1/sqrt(1+t^2) , unit sphere 
                
                ZA = [tau, one, 0]
                ZB = [-tau, one, 0]
                ZC = [-tau, -one, 0]
                ZD = [tau, -one, 0]
                YA = [one, 0 , tau]
                YB = [one, 0 , -tau]
                YC = [-one, 0 , -tau] 
                YD = [-one, 0 , tau]
                XA = [0 , tau, one]
                XB = [0 , -tau, one]
                XC = [0 , -tau, -one]
                XD = [0 , tau, -one]

                self.icosahedron.append(Triangle([YA, XA, YD]))
                self.icosahedron.append(Triangle([YA, YD, XB]))
                self.icosahedron.append(Triangle([YB, YC, XD]))
                self.icosahedron.append(Triangle([YB, XC, YC]))
                self.icosahedron.append(Triangle([ZA, YA, ZD]))
                self.icosahedron.append(Triangle([ZA, ZD, YB]))
                self.icosahedron.append(Triangle([ZC, YD, ZB]))
                self.icosahedron.append(Triangle([ZC, ZB, YC]))
                self.icosahedron.append(Triangle([XA, ZA, XD]))
                self.icosahedron.append(Triangle([XA, XD, ZB]))
                self.icosahedron.append(Triangle([XB, XC, ZD]))
                self.icosahedron.append(Triangle([XB, ZC, XC]))
                self.icosahedron.append(Triangle([XA, YA, ZA]))
                self.icosahedron.append(Triangle([XD, ZA, YB]))
                self.icosahedron.append(Triangle([YA, XB, ZD]))
                self.icosahedron.append(Triangle([YB, ZD, XC]))
                self.icosahedron.append(Triangle([YD, XA, ZB]))
                self.icosahedron.append(Triangle([YC, ZB, XD]))
                self.icosahedron.append(Triangle([YD, ZC, XB]))
                self.icosahedron.append(Triangle([YC, XC, ZC]))

                self.polygons.append(self.icosahedron)
        
                #with the starting triangles of the icoseahedron, start to divide each triangle in four triangles
                #by a level of 2, each starting triangle will be splitted in four parts and these four parts again in four parts
                #more by level 2, therefore, giving 4x4x20=320 faces => (4^n)*20
                #level 3 => 1280, level 4 => 5120, level 5 => 20480

                for l in xrange(self.max_levels): 
                        #print l
                        #print len(self.polygons[l])
                        self.polygons.append(list())
                        for tri in self.polygons[l]:
                                #print tri
                                #subdivide triangle into 4 triangles
                                subdivision = tri.subdivide()
                                self.polygons[l+1].append(subdivision[0])
                                self.polygons[l+1].append(subdivision[1])
                                self.polygons[l+1].append(subdivision[2])
                                self.polygons[l+1].append(subdivision[3])
                                #print subdivision[0].normal, " ", subdivision[1].normal, " ", subdivision[2].normal , " ", subdivision[3].normal

                print "number of cells in tesselated sphere:", len(self.polygons[self.max_levels])
                
                self.resolution = 2.1993 / sqrt(len(self.polygons[self.max_levels])) * 180 / pi * sqrt(2)
                return

        def findCell(self, normal):
                #find best cell in icosahedron aproximation

                maxDot = -1000
                cell = 0

                for tri in self.icosahedron:
                        dotP = dot(tri.normal, normal)
                        if dotP > maxDot:
                                maxDot = dotP
                                cell = tri

                #RECURSIVE METHOD
                #refinedCell = cell.findCell(normal)
                #return refinedCell
                
                #ITERATIVE METHOD (its faster, almost 25% faster)
                #still very sloowww!!!! use loookup table

                while len(cell.childs) > 0:
                        maxDot = -1000
                        selCell = 0
        
                        for chi in cell.childs:
                                dotP = dot(normal,chi.normal)
                                if dotP > maxDot:
                                        maxDot = dotP
                                        selCell = chi
        
                        cell = selCell
                        
                return cell

        def supportingVertices(self, normalCenter, sphereCell, dist, factor):
                supportVerts = list()
                for ver in sphereCell.mappedNormals:
                        normal = array(normalCenter[1]) / linalg.norm(array(normalCenter[1]))
                        if abs(dot(normal,(array(ver[0]) - array(normalCenter[0])))) < (dist/factor):
                                supportVerts.append(ver)

                normalSum = array([0, 0, 0])
                centerSum = array([0, 0, 0])
                
                for ver in supportVerts:
                        normalSum = normalSum + ver[1]
                        centerSum = centerSum + ver[0]

                return (normalSum/len(supportVerts),supportVerts,centerSum/len(supportVerts))
        
        def getOrderedTuple(self, pt1, pt2):
                if pt1 > pt2:
                        return pt2,pt1
                return pt1,pt2
        
        def emptyCells(self):
                nonEmpty=0
                for cell in self.polygons[self.max_levels]:
                        if len(cell.mappedNormals) > 0:
                                nonEmpty += 1
                                        
                return len(self.polygons[self.max_levels]) - nonEmpty

        def SplitNormals180Degrees(self, normal, point):
                for nc in point:
                        angle = arccos( dot(nc, normal) / (linalg.norm(nc)*linalg.norm(nc1))) * 180 / math.pi
        
        def computeConvexHull(self, polyDataNormals):
                length = polyDataNormals.GetPointData().GetNumberOfTuples()
                
                print "number of normals:", length
                
                normals = polyDataNormals.GetPointData().GetNormals()
                f = open('coordinates.txt', 'w')
                f.write('3\n')
                f.write(str(length) + '\n')
                
                k=0
                indicesPoints = []
                while (k < length):
                        center = polyDataNormals.GetPoint(k)
                        normal = normals.GetTuple3(k)
                        
                        strCoordinates = str(center[0]) + " " + str(center[1]) + " " + str(center[2]) + "\n"
                        #print strCoordinates
                        f.write(strCoordinates)
                        indicesPoints.insert(k, center)
                        k+=1
                
                #call qconvex
                
                os.system("qconvex Qt o < coordinates.txt > convexHull.off")
                os.system("meshconv convexHull.off -c ply -tri")
                                
        def computeEGIFromVTKNormals(self, polyDataNormals, CoM):
                #print polyDataNormals
                #print polyDataNormals.GetCellData() #contains the normals of the cells (polygons -> triangles)
                #print polyDataNormals.GetPointData() #contains the normals of the vertices
                
                #some vertices will have more than one normal associated, if these associated normals
                #have an angle of 90 degrees between them, generate normals between them with a certain
                #discretization
                
                self.useNormalsOnSurface = 1
                if self.useNormalsOnSurface:
                        length = polyDataNormals.GetCellData().GetNumberOfTuples()
                        normals = polyDataNormals.GetCellData().GetNormals()
                        k=0
                        while (k < length):
                                triangle = polyDataNormals.GetCell(k)
                                pts = triangle.GetPoints()
                                points = pts.GetData()
                                pt1 = points.GetTuple3(0) 
                                pt2 = points.GetTuple3(1)
                                pt3 = points.GetTuple3(2)
                                center = [0,0,0]
                                normal = normals.GetTuple3(k)
        
                                triangle.TriangleCenter(pt1, pt2, pt3, center)
                                cell = self.findCell(normal)
                                cell.addInfo(center,normal, 1, k)
                                k+=1
                
                return self.polygons[self.max_levels]

        def computeCenterOfMass(self, polydata):
                #print polydata
                comx = 0
                comy = 0
                comz = 0
                
                cells = polydata.GetNumberOfCells()
                totalArea = 0
                for k in range(cells):
                        triangle = polydata.GetCell(k) #each cell represents a triangle in the mesh
                        pts = triangle.GetPoints()
                        points = pts.GetData()
                        pt1 = points.GetTuple3(0) 
                        pt2 = points.GetTuple3(1)
                        pt3 = points.GetTuple3(2)                                                       
                        center = [0,0,0]
                        triangle.TriangleCenter(pt1, pt2, pt3, center)
                        area = triangle.TriangleArea(pt1,pt2,pt3)
                        comx += center[0] * area
                        comy += center[1] * area
                        comz += center[2] * area
                        
                        totalArea += area
                        
                CoM = [comx/totalArea, comy/totalArea, comz/totalArea]  
                return CoM

---

vfh_recognition
Author(s): Aitor Aldoma
autogenerated on Tue Mar 5 15:13:23 2013