From 96404a99cb5edddce8951442fd4b35b1bd2e41b4 Mon Sep 17 00:00:00 2001 From: wukevin Date: Fri, 17 Oct 2014 16:19:28 -0700 Subject: [PATCH 1/4] Version 3.1 of Kevin's Grid-Building method. Makes a grid, and identifies each point as inside/outside/surface point. --- autopack/Grid Prototype/Grid3.1.py | 715 +++++++++++++++++++++++++++++ 1 file changed, 715 insertions(+) create mode 100644 autopack/Grid Prototype/Grid3.1.py diff --git a/autopack/Grid Prototype/Grid3.1.py b/autopack/Grid Prototype/Grid3.1.py new file mode 100644 index 0000000..ff8ac76 --- /dev/null +++ b/autopack/Grid Prototype/Grid3.1.py @@ -0,0 +1,715 @@ +# import sys +# sys.path.append('C:\Users\Kevin\AppData\Roaming\MAXON\CINEMA 4D R15 Student_0471F205\plugins\ePMV\mgl64\MGLToolsPckgs\AutoFill') + +import c4d +from math import sqrt, ceil +# import AutoFill +import autopack +from autopack import ingr_ui as UI +# helper = AutoFill.helper +helper = None +if helper is None: + import upy + helperClass = upy.getHelperClass() + helper = helperClass() + # AutoFill.helper = helper +import numpy +import heapq +""" +1. Build grid using radius of gyration of smallest guy we want to pack. + - SphereTreeMaker (Rmin is radius of gyration, Rmax is the encapsulating radius). +2. Walk through triangles (that make up A,B,C) that make up the containers, assign triangles + to closest grid points at the resolution of the grid. + for triangle in polyhedron: + if lengthOfAnySideLongerThanGridSpace: + subdivide + # Break any edge longer into segments max size of unit + # Marching cubes? + # Coarse projetion of our polyhedron onto our grid (at resolution of grid) + + +Plot computation time vs number of vertices/polygons in the original polyhedron + +""" + +class gridPoint: + def __init__(self,i,globalC,isPolyhedron): + self.index = int(i) + self.isOutside = None + self.minDistance = 99999 # Only store a number here if within certain distance from polyhedron + self.representsPolyhedron = isPolyhedron + self.closeFaces = [] + self.closestFaceIndex = 0 + self.testedEndpoint = None + self.allDistances = [] # Stores a tuple list of distances to all points. (point,distance) = (5,2.5) + self.globalCoord = numpy.array(globalC) # Stores the global coordinate associated with this point + +def visualizeAtomArray(name = 'Test', vertexCoordinates = [(0,0,0)], sphereRadius = 5): + """ + Given coordinates, visualize them in Cinema 4D using an Atom Array Object. + """ + atom = c4d.BaseObject(c4d.Oatomarray) + atom.SetName(name) + atom[1000] = 3 #radius cylinder + atom[1001] = sphereRadius #radius sphere + atom[1002] = 5 #subdivision + helper.addObjectToScene(helper.getCurrentScene(),atom,parent=None) + + parent = atom + coords = vertexCoordinates + nface = 0 + obj= c4d.PolygonObject(len(coords), nface) + obj.SetName(name + '_child') + cd4vertices = map(helper.FromVec,coords) + map(obj.SetPoint,range(len(coords)),cd4vertices) + helper.addObjectToScene(helper.getCurrentScene(),obj,parent=parent) + +def setupBBOX(object_target): + """ + Sets up a bounding box that encompasses a selected object in Cinema 4D. Returns the bounding box + object. + """ + mesh = object_target + + faces,vertices,vnormals,fn = helper.DecomposeMesh(mesh,edit=False,copy=False,tri=True,transform=True,fn=True) + + from autopack.Compartment import Compartment + o1 = Compartment(helper.getName(mesh),vertices, faces, vnormals,fnormals=fn) + o1.ref_obj = mesh + o1.number = 1 + + b = helper.Box("BBOX", cornerPoints=o1.bb) + bb=o1.bb + return o1 + +def makeGrid(lowerLeft,upperRight,spacing, inclusive = True): + """ + Given lower left back corner, an upper right front corner, and a grid spacing, construct a list of + coordinates that representsly an evenly spaced grid spanning the entire volume. Returned coordinates + preserve the right/lefthandedness of the input coordinates, making this function work equally well + in both 3D coordinate systems. + """ + pointCoords = [] + xSize,ySize,zSize = upperRight[0] - lowerLeft[0], upperRight[1] - lowerLeft[1], upperRight[2] - lowerLeft[2] + if inclusive: + nx,ny,nz = xSize / spacing + 1, ySize / spacing + 1, zSize / spacing + 1 + else: + nx,ny,nz = xSize / spacing, ySize / spacing, zSize / spacing + nx,ny,nz = round(nx),round(ny),round(nz) + x = numpy.linspace(lowerLeft[0],upperRight[0],nx) + y = numpy.linspace(lowerLeft[1],upperRight[1],ny) + z = numpy.linspace(lowerLeft[2],upperRight[2],nz) + for a in x: + for b in y: + for c in z: + pointCoords.append((a,b,c)) + # Order of above had to change to correctly produce a viable BBOX + return pointCoords, (nx,ny,nz) + +def getPointFrom3D(pt3d, gridSpacing, gridPtsPerEdge, BBOX): + """ + Starts from orthogonal grid. Has a lower left, upper right coordinate from bounding box. + Given lower left and known spacing of grid, we can provide any vector in 3D space and it will + algebraically return index of closest grid point. + + Given a RIGHT-HANDED fine coordinate, it will return the index of the closest point in the LEFT-HANDED coarse grid. + """ + x, y, z = pt3d # Continuous 3D point to be discretized + spacing1 = 1./gridSpacing # Grid spacing = diagonal of the voxel determined by smalled packing radius + NZ, NY, NX = gridPtsPerEdge # vector = [length, height, depth] of grid, units = gridPoints + + # Added to fix some weird edge behavior + #NZ, NY, NX = round(NZ),round(NY),round(NX) + #print(NX,NY,NZ) + OZ, OY, OX = BBOX.getBoundingBox()[0] # origin of Pack grid. OX and OZ were switched to compensate for coordinate systems + #print(OX,OY,OZ) + + # Algebra gives nearest gridPoint ID to pt3D + i = min( NX-1, max( 0, round((x-OX)*spacing1))) + j = min( NY-1, max( 0, round((y-OY)*spacing1))) + k = min( NZ-1, max( 0, round((z-OZ)*spacing1))) + result = k*NX*NY + j*NX + i + return int(result) + #return int(k*NX*NY + j*NX + i) + +def makeMarchingCube(gridSpacing,r): + """ + Create a numpy array that represents the precomputed distances to each point + for the cube of points surrounding our center point. + """ + def _pythagorean(*edgeLengths): + from math import sqrt + total = 0 + for length in edgeLengths: + total += length * length + distance = sqrt(total) + return distance + from math import ceil + pointsForRadius = ceil(r/gridSpacing) # Number of grid points required to represent our radius, rounded up + pointsInEdge = 2 * pointsForRadius + 1 # Number of points in one edge of our cube + + center = pointsForRadius # The index if our center point + cube = numpy.zeros(shape=(pointsInEdge,pointsInEdge,pointsInEdge)) + distX = numpy.zeros(shape=(pointsInEdge,pointsInEdge,pointsInEdge)) + distY = numpy.zeros(shape=(pointsInEdge,pointsInEdge,pointsInEdge)) + distZ = numpy.zeros(shape=(pointsInEdge,pointsInEdge,pointsInEdge)) + for a in range(pointsInEdge): + lenX = a - center + for b in range(pointsInEdge): + lenY = b - center + for c in range(pointsInEdge): + lenZ = c - center + cube[a][b][c] = _pythagorean(lenX,lenY,lenZ) * gridSpacing + distX[a][b][c] = lenX + distY[a][b][c] = lenY + distZ[a][b][c] = lenZ + return cube, distX, distY, distZ + +def f_dot_product(vector1,vector2): + """ + Return the dot product of two 3D vectors. + """ + dottedVectors = [vector1[i] * vector2[i] for i in range(len(vector1))] + return sum(dottedVectors) + +def vcross(v1,v2): + """ + Returns the cross-product of two 3D vectors. + """ + x1,y1,z1 = v1 + x2,y2,z2 = v2 + return (y1*z2-y2*z1, z1*x2-z2*x1, x1*y2-x2*y1) + +def vlen(v1): + """ + Returns the length of a 3D vector. + """ + x1,y1,z1 = v1 + return sqrt(x1*x1 + y1*y1 + z1*z1) + +def getCubeIndices(leftHandedCoord,gridPtsPerEdge,boundingBoxOrigin,gridSpacing,zippedNumbers): + """ + Given a left-handed coordinate, the lower left back corner of a volume, the number of grid points + per edge, the distance between grid points, and the precomputed distance cube, it returns the + indices of points that make up the cube surrounding a point. + """ + NZ,NY,NX = gridPtsPerEdge + OZ,OY,OX = boundingBoxOrigin + spacing1 = 1./gridSpacing + + zTemp,yTemp,xTemp = leftHandedCoord + i,j,k = round((xTemp-OX)*spacing1), round((yTemp-OY)*spacing1), round((zTemp-OZ)*spacing1) + cubeIndicies = [] + + for d,x,y,z in zippedNumbers: + newI, newJ, newK = i + x, j + y, k + z + if newI < 0 or newI > (NX-1) or newJ < 0 or newJ > (NY-1) or newK < 0 or newK > (NZ-1): + continue + desiredPointIndex = int(round(newK*NX*NY + newJ*NX + newI)) + cubeIndicies.append(desiredPointIndex) + return cubeIndicies + +def getOctahedronIndices(leftHandedCoord,gridPtsPerEdge,boundingBoxOrigin,gridSpacing): + """ + Similar to getCubeIndices, but isntead of returning indicies of the cube, it retunrs the indices + of the six points immediately surrounding the central point that make up an octahedron. + + This was originally written in an attempt to stremline the inside/outside testing pre-flood filling. + However, it turned out that the time it took to copute these 6 indices was slow enough to make + skipping a few inside/outside ray collision tests not worthwhile. Thus, this function is no longer in + use, but will stay here for booking purposes. + """ + NZ,NY,NX = gridPtsPerEdge + OZ,OY,OX = boundingBoxOrigin + spacing1 = 1./gridSpacing + + zTemp,yTemp,xTemp = leftHandedCoord + i,j,k = round((xTemp-OX)*spacing1), round((yTemp-OY)*spacing1), round((zTemp-OZ)*spacing1) + octahedronIndices = [] + zippedNumbers = [(0,0,1),(0,0,-1),(0,1,0),(0,-1,0),(1,0,0),(-1,0,0)] + for x,y,z in zippedNumbers: + newI, newJ, newK = i + x, j + y, k + z + if newI < 0 or newI > (NX-1) or newJ < 0 or newJ > (NY-1) or newK < 0 or newK > (NZ-1): + continue + desiredPointIndex = int(round(newK*NX*NY + newJ*NX + newI)) + octahedronIndices.append(desiredPointIndex) + return octahedronIndices + +def f_ray_intersect_polyhedron(pRayStartPos, pRayEndPos, faces,vertices, pTruncateToSegment): + """This function returns TRUE if a ray intersects a triangle. + It also calculates and returns the UV coordinates of said colision as part of the intersection test, + Makes sure that we are working with arrays + + * TAKES IN pRayStartPos as LEFT-HANDED COORDINATE (Z,Y,X) + * TAKES IN pRayEndPos AS LEFT-HANDED COORDINATE (Z,Y,X) + * faces is a list of vertex indices that defines the polyhedron + * vertices is the list of global coordinates that faces refers to + * pTruncateToSegment decides whether or not the segment terminates at the end position, or keeps on + going forever + """ + pRayStartPos = numpy.array(pRayStartPos) + pRayEndPos = numpy.array(pRayEndPos) + + vLineSlope = pRayEndPos - pRayStartPos #This line segment defines an infinite line to test for intersection + vPolyhedronPos = numpy.array((0,0,0)) + vTriPoints = vertices + vTriPolys = faces + vBackface = None + vBackfaceFinal = None + maxDistance = 99999 + vHitCount = 0 + + # Globalize the polyhedron + # Present in original coffee script, but unnecessary because DecomposeMesh gives global coords + #for i in range(len(vTriPoints)): + # vTriPoints[i] = vTriPoints[i] + vPolyhedronPos + + vEpsilon = 0.00001 + vBreakj = False + vCollidePos = None + + # Walk through each polygon in a polyhedron + for testingFace in vTriPolys: + # Loop through all the polygons in an input polyhedron + #vQuadrangle = 1 # Default says polygon is a quadrangle. + #vLoopLimit = 2 # Default k will loop through polygon assuming its a quad. + #if (vTriPolys[j+3] == vTriPolys[j+2]) # Test to see if quad is actually just a triangle. + # { + vQuadrangle = 0 # Current polygon is not a quad, it's a triangle. + vLoopLimit = 1 # Set k loop to only cycle one time. + + for k in range(vLoopLimit): + vTriPt0 = numpy.array(vTriPoints[testingFace[0]]) # Always get the first point of a quad/tri + vTriPt1 = numpy.array(vTriPoints[testingFace[1]]) # Get point 1 for a tri and a quad's first pass, but skip for a quad's second pass + vTriPt2 = numpy.array(vTriPoints[testingFace[2]]) # Get point 2 for a tri and a quad's first pass, but get point 3 only for a quad on its second pass. + + vE1 = vTriPt1 - vTriPt0 # Get the first edge as a vector. + vE2 = vTriPt2 - vTriPt0 # Get the second edge. + h = vcross(vLineSlope, vE2) + + a = f_dot_product(vE1,h) # Get the projection of h onto vE1. + if a > -vEpsilon and a < vEpsilon: + continue# If the ray is parallel to the plane then it does not intersect it, i.e, a = 0 +/- given rounding slope. + # If the polygon is a quadrangle, test the other triangle that comprises it. + + F = 1.0/a + s = pRayStartPos - vTriPt0 # Get the vector from the origin of the triangle to the ray's origin. + u = F * f_dot_product(s,h) + if u < 0.0 or u > 1.0: + continue + #/* Break if its outside of the triangle, but try the other triangle if in a quad. + #U is described as u = : start of vE1 = 0.0, to the end of vE1 = 1.0 as a percentage. + #If the value of the U coordinate is outside the range of values inside the triangle, + #then the ray has intersected the plane outside the triangle.*/ + + q = vcross(s, vE1) + v = F * f_dot_product(vLineSlope,q) + if v <0.0 or u+v > 1.0: + continue + #/* Break if outside of the triangles v range. + #If the value of the V coordinate is outside the range of values inside the triangle, + #then the ray has intersected the plane outside the triangle. + #U + V cannot exceed 1.0 or the point is not in the triangle. + #If you imagine the triangle as half a square this makes sense. U=1 V=1 would be in the + #lower left hand corner which would be in the second triangle making up the square.*/ + + vCollidePos = vTriPt0 + u*vE1 + v*vE2 # This is the global collision position. + assert len(vCollidePos) == 3 + + # The ray is hitting a triangle, now test to see if its a triangle hit by the ray. + vBackface = False + if f_dot_product(vLineSlope, vCollidePos - pRayStartPos) > 0: # This truncates our infinite line to a ray pointing from start THROUGH end positions. + vHitCount += 1 + #print(testingFace) + if pTruncateToSegment and vlen(vLineSlope) < vlen(vCollidePos - pRayStartPos): + print('broken') + break # This truncates our ray to a line segment from start to end positions. + + d = squaredTwoPointDistance(pRayStartPos,vCollidePos) + if d >= maxDistance: + continue + if a < vEpsilon: # Test to see if the triangle hit is a backface. + #set master grid to organelle->getname inside + vBackface = True + # This stuff is specific to our Point inside goals. + vBreakj = True # To see if a point is inside, I can stop at the first backface hit. + #break + else: + vBreakj = True + vBackfaceFinal = vBackface + maxDistance = d + #vBackfaceFinal = vBackfaces[distances.index(min(distances))] + return vHitCount,vBackfaceFinal + +def _findTriangleCenter(*args): + """ + Given three coordinates, it calculates the center of the triangle by averaging + the coordinates in each axes. + """ + p1,p2,p3 = args[0] + x = zip(p1,p2,p3) + result = [] + for elem in x: + result.append(sum(elem)/len(elem)) + return result + +def findPointsCenter(*args): + # Average down the column, such that we're averaging across all measurements in one dimension + center = numpy.mean(args[0], axis = 0) + return center + +def squaredTwoPointDistance(point1,point2): + """Computes and returns the distance between two points in 3D space""" + dists = [numpy.float64(point1[i] - point2[i]) for i in range(len(point1))] + distsSquare = [x**2 for x in dists] + dist = sum(distsSquare) + return dist + +def projectPolyhedronToGrid(gridSpacing): + startTime = time() + from math import ceil + # Get the object, and decompose it + object_target = helper.getCurrentSelection()[0] + faces,verticesLH,vnormals,faceNormals = helper.DecomposeMesh(object_target,edit=False,copy=False,tri=True,transform=True,fn=True) + # Setup a bounding box for our polyhedron + b = setupBBOX(object_target) + corners = b.getBoundingBox() + + # Create grid points to fill in our bbox with grid points at regular distances + points, gridPtsPerEdge = makeGrid(corners[0],corners[1],gridSpacing) # These are left-handed + + # Create a list of all our gridPoint objects + gridPoints = [] + i = 0 + for point in points: + gridPoints.append(gridPoint(i,point,isPolyhedron = False)) + i += 1 + assert len(gridPoints) == len(points) + + # Create grid points to fill in our bbox with grid points at regular distances + points, gridPtsPerEdge = makeGrid(corners[0],corners[1],gridSpacing) # These are left-handed + + # For every face in the polyhedron + allCoordinates = [] + for face in faces: + # Dist is a small function that calculates the distance between two points. + dist = lambda x,y: vlen([y[i] - x[i] for i in range(len(x))]) + # Get the vertex coordinates and conver to numpy arrays... + triCoords = [numpy.array(verticesLH[i]) for i in face] + allCoordinates.extend(triCoords) + # ...use them to define the u/v vectors + pos = triCoords[0] + u = triCoords[1] - pos + v = triCoords[2] - pos + # SOmetimes the hypotenuse isn't fully represented. To remedy this, we will use an addition w vector + w = triCoords[2] - triCoords[1] + + # If either u or v is greater than the grid spacing, then we need to subdivide it + # We will use ceil: if we have a u of length 16, and grid spacing of 5, then we want + # a u at 0, 5, 10, 15 which is [0, 1, 2, 3] * gridSpacing. We need ceil to make this happen. + + # This is a much more efficient solution than the previous method of artificially increasing + # the densith of the polygon mesh by 4, then mapping each of those small points to the grid. + # However, it is also leaks more, for (as of yet) unknown reasons. + + # Minimum is one because range(1) gives us [0] + uSubunits, vSubunits, wSubunits = 1, 1, 1 + if vlen(u) > gridSpacing: + uSubunits = ceil(vlen(u)/gridSpacing) + if vlen(v) > gridSpacing: + vSubunits = ceil(vlen(v)/gridSpacing) + if vlen(w) > gridSpacing: + wSubunits = ceil(vlen(w)/gridSpacing) + # Because we have observed leakage, maybe we want to try trying a denser interpolation, using numpy's linspace? + for uSub in range(uSubunits): + percentU = uSub * gridSpacing / vlen(u) + assert percentU < 1.0 # Make sure that we have not stepped outside of our original u vector + # h represents the height of the hypotenuse at this u. We cannot go past the hypotenuse, so this will be + # our upper bound. + h = percentU * u + (1 - percentU) * v + for vSub in range(vSubunits): + percentV = vSub * gridSpacing / vlen(v) + assert percentV < 1.0 # Make sure that we have not stepped oustide of our original v vector. + interpolatedPoint = percentU * u + percentV * v + # The original if: statement asks if the distance from the origin to the interpolated point is less than + # the distance from the origin to the hypotenuse point. + # if vlen(interpolatedPoint) < vlen(h): + # Wouldn't it be a better idea to measure distance to the u position instead? + if (vlen(interpolatedPoint - percentU * u) < vlen(h - percentU * u)): + allCoordinates.append(interpolatedPoint + pos) + else: + # Used to be a break, changed to continue (not 100% sure if it made too much difference) + continue + # Because the above only interpolates the face, and may not completely capture the hypotenuse, let's separately + # interpolate points on the hypotenuse (the w vector) + for wSub in range(wSubunits): + percentW = wSub * gridSpacing / vlen(w) + if percentW > 1.0: + percentW = 1.0 + interpolatedPoint = percentW * w + allCoordinates.append(interpolatedPoint + triCoords[1]) + visualizeAtomArray('fineCoords',allCoordinates,5) + projectedIndices = set() + for coord in allCoordinates: + projectedPointIndex = getPointFrom3D(coord[::-1],gridSpacing,gridPtsPerEdge,b) + projectedIndices.add(projectedPointIndex) + visualizeAtomArray('projectedCoords',[points[i] for i in projectedIndices],5) + print('Projecting polyhedron to grid took ' + str(time() - startTime) + ' seconds.') + return(projectedIndices) + +def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): + """ + Takes a polyhedron, and builds a grid. In this grid: + - Projects the polyhedron to the grid. + - Determines which points are inside/outside the polyhedron + - Determines point's distance to the polyhedron. + Usage: + Select desired object, run this program with desired grid spacing resolution. By default, + radius is 2 * gridSpacing. superFine provides the option doing a super leakproof test when + determining which points are inside or outside. This usually not necessary, because the + built-in algorithm has no known leakage cases. It is simply there as a safeguard. + """ + from time import time + if radius == None: + radius = gridSpacing + + startTime = time() + from math import ceil + # Get the object, and decompose it + object_target = helper.getCurrentSelection()[0] + faces,verticesLH,vnormals,faceNormals = helper.DecomposeMesh(object_target,edit=False,copy=False,tri=True,transform=True,fn=True) + # Setup a bounding box for our polyhedron + b = setupBBOX(object_target) + corners = b.getBoundingBox() + + # Create grid points to fill in our bbox with grid points at regular distances + points, gridPtsPerEdge = makeGrid(corners[0],corners[1],gridSpacing) # These are left-handed + + # Create a list of all our gridPoint objects + gridPoints = [] + i = 0 + for point in points: + gridPoints.append(gridPoint(i,point,isPolyhedron = False)) + i += 1 + assert len(gridPoints) == len(points) + + distanceCube,distX,distY,distZ = makeMarchingCube(gridSpacing,radius) + # Flatten and combine these arrays, so that it's easier to iterate over + distanceCubeF,distXF,distYF,distZF = distanceCube.flatten(),distX.flatten(),distY.flatten(),distZ.flatten() + zippedNumbers = zip(distanceCubeF,distXF,distYF,distZF) + + NZ,NY,NX = gridPtsPerEdge + OZ, OY, OX = b.getBoundingBox()[0] + spacing1 = 1./gridSpacing + # For every face in the polyhedron + allCoordinates = [] + pointsToTestInsideOutside = set() + # Walk through the faces + for face in faces: + # Dist is a small function that calculates the distance between two points. + dist = lambda x,y: vlen([y[i] - x[i] for i in range(len(x))]) + # Get the vertex coordinates and conver to numpy arrays... + triCoords = [numpy.array(verticesLH[i]) for i in face] + thisFaceFineCoords = list(triCoords) + allCoordinates.extend(triCoords) + # ...use them to define the u/v vectors + pos = triCoords[0] + u = triCoords[1] - pos + v = triCoords[2] - pos + # SOmetimes the hypotenuse isn't fully represented. To remedy this, we will use an addition w vector + w = triCoords[2] - triCoords[1] + + # If either u or v is greater than the grid spacing, then we need to subdivide it + # We will use ceil: if we have a u of length 16, and grid spacing of 5, then we want + # a u at 0, 5, 10, 15 which is [0, 1, 2, 3] * gridSpacing. We need ceil to make this happen. + + # This is a much more efficient solution than the previous method of artificially increasing + # the densith of the polygon mesh by 4, then mapping each of those small points to the grid. + # However, it is also leaks more, for (as of yet) unknown reasons. + + # It was discovered that by using the default gridspacing, some faces will produce leakage. One solution is to + # temporarily use a somewhat denser gridspacing to interpolate, and then project back onto our original spacing. + # We'll decrease the gridspacing by 25% (so that it's 75% of the original). This seems to fix the issue without + # any appreciable increase in computation time. + gridSpacingTempFine = gridSpacing * 3 / 4 + # Minimum is one because range(1) gives us [0] + uSubunits, vSubunits, wSubunits = 1, 1, 1 + if vlen(u) > gridSpacingTempFine: + uSubunits = ceil(vlen(u)/gridSpacingTempFine) + 1 + if vlen(v) > gridSpacingTempFine: + vSubunits = ceil(vlen(v)/gridSpacingTempFine) + 1 + if vlen(w) > gridSpacingTempFine: + wSubunits = ceil(vlen(w)/gridSpacingTempFine) + 1 + # Because we have observed leakage, maybe we want to try trying a denser interpolation, using numpy's linspace? + for uSub in range(uSubunits): + percentU = uSub * gridSpacingTempFine / vlen(u) + percentU = min(percentU, 1.0) # Make sure that we have not stepped outside of our original u vector + # h represents the height of the hypotenuse at this u. We cannot go past the hypotenuse, so this will be + # our upper bound. + h = percentU * u + (1 - percentU) * v + for vSub in range(vSubunits): + percentV = vSub * gridSpacingTempFine / vlen(v) + percentV = min(percentV, 1.0) # Make sure that we have not stepped oustide of our original v vector. + interpolatedPoint = percentU * u + percentV * v + # The original if: statement asks if the distance from the origin to the interpolated point is less than + # the distance from the origin to the hypotenuse point. + # if vlen(interpolatedPoint) < vlen(h): + # Wouldn't it be a better idea to measure distance to the u position instead? + if (vlen(interpolatedPoint - percentU * u) < vlen(h - percentU * u)): + allCoordinates.append(interpolatedPoint + pos) + thisFaceFineCoords.append(interpolatedPoint + pos) + else: + # Used to be a break, changed to continue (not 100% sure if it made too much difference) + break + # Because the above only interpolates the face, and may not completely capture the hypotenuse, let's separately + # interpolate points on the hypotenuse (the w vector) + for wSub in range(wSubunits): + percentW = wSub * gridSpacingTempFine / vlen(w) + percentW = min(percentW, 1.0) + interpolatedPoint = percentW * w + allCoordinates.append(interpolatedPoint + triCoords[1]) + thisFaceFineCoords.append(interpolatedPoint + triCoords[1]) + projectedIndices = set() + for coord in thisFaceFineCoords: + projectedPointIndex = getPointFrom3D(coord[::-1],gridSpacing,gridPtsPerEdge,b) + projectedIndices.add(projectedPointIndex) + + for P in projectedIndices: + # Get the point object corresponding to the index, and set its polyhedron attribute to true + g = gridPoints[P] + g.representsPolyhedron = True + # Get the coordinates of the point, and convert them to grid units + zTemp,yTemp,xTemp = g.globalCoord + i,j,k = round((xTemp-OX)*spacing1), round((yTemp-OY)*spacing1), round((zTemp-OZ)*spacing1) + # Let's step through our distance cube, and assign faces/closest distances to each + for d,x,y,z in zippedNumbers: + # Get the grid indices for the point we're considering, and pass if we're stepping oustide the boundaries + newI, newJ, newK = i + x, j + y, k + z + if newI < 0 or newI > (NX-1) or newJ < 0 or newJ > (NY-1) or newK < 0 or newK > (NZ-1): + continue + # Get the point index that this coordinate corresponds to. + desiredPointIndex = int(round(newK*NX*NY + newJ*NX + newI)) + desiredPoint = gridPoints[desiredPointIndex] + # Add the current face to the its list of closest faces + if face not in desiredPoint.closeFaces: + desiredPoint.closeFaces.append(face) + # Add the distance to the point's list of distances, and overwrite minimum distance if appropriate + desiredPoint.allDistances.append((v,d)) + if d < desiredPoint.minDistance: + desiredPoint.minDistance = d + desiredPoint.closestFaceIndex = len(desiredPoint.closeFaces) - 1 + # Later down the road, we want to test as few points as possible for inside/outside. Therefore, + # we will only test points that are + if abs(x) <= 1 and abs(y) <= 1 and abs(z) <= 1: + pointsToTestInsideOutside.add(desiredPointIndex) + visualizeAtomArray('fineCoords',allCoordinates,5) + projectedIndices = [x.index for x in gridPoints if x.representsPolyhedron == True] + # projectedIndices = set() + # for coord in allCoordinates: + # projectedPointIndex = getPointFrom3D(coord[::-1],gridSpacing,gridPtsPerEdge,b) + # projectedIndices.add(projectedPointIndex) + # visualizeAtomArray('projectedCoords',[points[i] for i in projectedIndices],5) + print('Projecting polyhedron to grid took ' + str(time() - startTime) + ' seconds.') + + # Let's start filling in inside outside. Here's the general algorithm: + # Walk through all the points in our grid. Once we encounter a point that has closest faces, + # then we know we need to test it for inside/outside. Once we test that for inside/outside, we + # fill in all previous points with that same inside outisde property. To account for the possible + # situation that there is a surface that is only partially bound by the bbox, then we need to + # reset the insideOutsideTracker every time we have a change in more than 1 of the 3 coordinates + # because that indicates we're starting a new row/column of points. + + # This tracks our inside/outside. + isOutsideTracker = None + # This tracks the points that we've iterated over which we do not know if inside/outside. Remember + # to reset every time we find an inside/outside. + emptyPointIndicies = [] + mismatchCounter = 0 + for g in gridPoints: + # Check if we've started a new line. If so, then we reset everything. + # This test should precede all other test, because we don't want old knowldge + # to carry over to the new line. + if g.index > 0: # We can't check the first element, so we can skip it. + coordDiff = g.globalCoord - gridPoints[g.index - 1].globalCoord + coordDiffNonzero = [x != 0 for x in coordDiff] + if sum(coordDiffNonzero) > 1: + # assert len(emptyPointIndicies) == 0 # When starting a new line, we shouldn't have any unknowns from the previous line + isOutsideTracker = None + emptyPointIndicies = [] + + # There's no point testing inside/outside for points that are on the surface. So if + # we get to one, we just skip over it. + if g.representsPolyhedron == True: + g.isOutside = None + continue + + if len(g.closeFaces) == 0: + # If it's not close to any faces, and we don't know if this row is inside/outside, then + # we have to wait till later to figure it out + if isOutsideTracker == None: + emptyPointIndicies.append(g.index) + # However, if we do know , we can just use the previous one to fill + else: + g.isOutside = isOutsideTracker + + # If there are close faces attached to it, then we need to test it for inside/outside. + else: + # Find centroid of all the vertices of all the close faces. This will be our endpoint + # when casting a ray for collision testing. + uniquePoints = [] + # This takes just the first face and projects to the center of it. + # [uniquePoints.append(x) for x in g.closeFaces[0] if x not in uniquePoints] + [uniquePoints.append(x) for x in g.closeFaces[g.closestFaceIndex] if x not in uniquePoints] + uniquePointsCoords = verticesLH[uniquePoints] + endPoint = findPointsCenter(uniquePointsCoords) + g.testedEndpoint = endPoint + + # Draw a ray to that point, and see if we hit a backface or not + numHits, thisBackFace = f_ray_intersect_polyhedron(g.globalCoord,endPoint,g.closeFaces,verticesLH,False) + + # We can check the face as well if we want to be super precise. If they dont' agree, we then check against the entire polyhedron. + if superFine == True: + if len(g.closeFaces) > 1: + uniquePoints2 = [] + [uniquePoints2.append(x) for x in g.closeFaces[1] if x not in uniquePoints2] + uniquePointsCoords2 = verticesLH[uniquePoints2] + endPoint2 = findPointsCenter(uniquePointsCoords2) + numHits2, thisBackFace2 = f_ray_intersect_polyhedron(g.globalCoord,endPoint2,g.closeFaces,verticesLH,False) + if len(g.closeFaces) == 1 or thisBackFace != thisBackFace2: + mismatchCounter += 1 + numHits, thisBackFace = f_ray_intersect_polyhedron(g.globalCoord,numpy.array([0.0,0.0,0.0]),faces,verticesLH,False) + + # Fill in inside outside attribute for this point, as pRayStartPos, pRayEndPos, faces, vertices, pTruncateToSegmentll as for any points before it + g.isOutside = not thisBackFace + isOutsideTracker = not thisBackFace + for i in emptyPointIndicies: + gridPoints[i].isOutside = isOutsideTracker + # Because we have filled in all the unknowns, we can reset that counter. + emptyPointIndicies = [] + + # Final pass through for quality/sanity checks. + for g in gridPoints: + if g.representsPolyhedron == True: + assert g.isOutside == None + continue + else: + if g.isOutside == None: + g.isOutside = True + + insidePoints = [g.index for g in gridPoints if g.isOutside == False] + outsidePoints = [g.index for g in gridPoints if g.isOutside == True] + surfacePoints = [g.index for g in gridPoints if g.representsPolyhedron == True] + + visualizeAtomArray('insidePoints',[points[i] for i in insidePoints],0.5) + visualizeAtomArray('outsidePoints',[points[i] for i in outsidePoints],0.5) + visualizeAtomArray('surfacePoints',[points[i] for i in surfacePoints],0.5) + print('Superfine was ' + str(superFine) + ' and there were ' + str(mismatchCounter) + ' mismatches.') + print('Grid construction took ' + str(time() - startTime) + ' seconds for ' + str(len(faces)) + ' faces and ' + str(len(gridPoints)) + ' points.') + + visualizeAtomArray('badPoint', [gridPoints[insidePoints[14088]].globalCoord]) + visualizeAtomArray('badPointEndPoint',[gridPoints[insidePoints[14088]].testedEndpoint]) + visualizeAtomArray('badPointCloseFaces',verticesLH[gridPoints[insidePoints[14088]].closeFaces]) + + return gridPoints + +projectPolyhedronToGrid2(5, superFine = False) From 3a182853e3efbcece9d9f1523ae6c122f6b0bf17 Mon Sep 17 00:00:00 2001 From: wukevin Date: Fri, 17 Oct 2014 17:24:39 -0700 Subject: [PATCH 2/4] Updated. Works perfectly on all structures that have no concavities. --- autopack/Grid Prototype/Grid3.1.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/autopack/Grid Prototype/Grid3.1.py b/autopack/Grid Prototype/Grid3.1.py index ff8ac76..f53ad31 100644 --- a/autopack/Grid Prototype/Grid3.1.py +++ b/autopack/Grid Prototype/Grid3.1.py @@ -507,7 +507,7 @@ def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): # Walk through the faces for face in faces: # Dist is a small function that calculates the distance between two points. - dist = lambda x,y: vlen([y[i] - x[i] for i in range(len(x))]) + # dist = lambda x,y: vlen([y[i] - x[i] for i in range(len(x))]) # Get the vertex coordinates and conver to numpy arrays... triCoords = [numpy.array(verticesLH[i]) for i in face] thisFaceFineCoords = list(triCoords) @@ -574,7 +574,7 @@ def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): projectedPointIndex = getPointFrom3D(coord[::-1],gridSpacing,gridPtsPerEdge,b) projectedIndices.add(projectedPointIndex) - for P in projectedIndices: + for P in list(projectedIndices): # Get the point object corresponding to the index, and set its polyhedron attribute to true g = gridPoints[P] g.representsPolyhedron = True @@ -706,9 +706,9 @@ def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): print('Superfine was ' + str(superFine) + ' and there were ' + str(mismatchCounter) + ' mismatches.') print('Grid construction took ' + str(time() - startTime) + ' seconds for ' + str(len(faces)) + ' faces and ' + str(len(gridPoints)) + ' points.') - visualizeAtomArray('badPoint', [gridPoints[insidePoints[14088]].globalCoord]) - visualizeAtomArray('badPointEndPoint',[gridPoints[insidePoints[14088]].testedEndpoint]) - visualizeAtomArray('badPointCloseFaces',verticesLH[gridPoints[insidePoints[14088]].closeFaces]) + # visualizeAtomArray('badPoint', [gridPoints[insidePoints[14067]].globalCoord]) + # visualizeAtomArray('badPointEndPoint',[gridPoints[insidePoints[14067]].testedEndpoint]) + # visualizeAtomArray('badPointCloseFaces',verticesLH[gridPoints[insidePoints[14067]].closeFaces]) return gridPoints From 5b800257e98d3c7556b78bef4ef8458079e4dc62 Mon Sep 17 00:00:00 2001 From: wukevin Date: Fri, 24 Oct 2014 14:42:17 -0700 Subject: [PATCH 3/4] Takes in a polyhedron in c4d, returns a grid with inside/outside/surface points annotated. Usage: projectPolyhedronToGrid2(x) where x is the grid spacing you want. Changes since 3.1: - Fixed leakage issues on concave sections - Code cleanup --- .../Grid Prototype/{Grid3.1.py => Grid3.2.py} | 243 ++++-------------- 1 file changed, 54 insertions(+), 189 deletions(-) rename autopack/Grid Prototype/{Grid3.1.py => Grid3.2.py} (72%) diff --git a/autopack/Grid Prototype/Grid3.1.py b/autopack/Grid Prototype/Grid3.2.py similarity index 72% rename from autopack/Grid Prototype/Grid3.1.py rename to autopack/Grid Prototype/Grid3.2.py index f53ad31..cbbf647 100644 --- a/autopack/Grid Prototype/Grid3.1.py +++ b/autopack/Grid Prototype/Grid3.2.py @@ -187,54 +187,6 @@ def vlen(v1): x1,y1,z1 = v1 return sqrt(x1*x1 + y1*y1 + z1*z1) -def getCubeIndices(leftHandedCoord,gridPtsPerEdge,boundingBoxOrigin,gridSpacing,zippedNumbers): - """ - Given a left-handed coordinate, the lower left back corner of a volume, the number of grid points - per edge, the distance between grid points, and the precomputed distance cube, it returns the - indices of points that make up the cube surrounding a point. - """ - NZ,NY,NX = gridPtsPerEdge - OZ,OY,OX = boundingBoxOrigin - spacing1 = 1./gridSpacing - - zTemp,yTemp,xTemp = leftHandedCoord - i,j,k = round((xTemp-OX)*spacing1), round((yTemp-OY)*spacing1), round((zTemp-OZ)*spacing1) - cubeIndicies = [] - - for d,x,y,z in zippedNumbers: - newI, newJ, newK = i + x, j + y, k + z - if newI < 0 or newI > (NX-1) or newJ < 0 or newJ > (NY-1) or newK < 0 or newK > (NZ-1): - continue - desiredPointIndex = int(round(newK*NX*NY + newJ*NX + newI)) - cubeIndicies.append(desiredPointIndex) - return cubeIndicies - -def getOctahedronIndices(leftHandedCoord,gridPtsPerEdge,boundingBoxOrigin,gridSpacing): - """ - Similar to getCubeIndices, but isntead of returning indicies of the cube, it retunrs the indices - of the six points immediately surrounding the central point that make up an octahedron. - - This was originally written in an attempt to stremline the inside/outside testing pre-flood filling. - However, it turned out that the time it took to copute these 6 indices was slow enough to make - skipping a few inside/outside ray collision tests not worthwhile. Thus, this function is no longer in - use, but will stay here for booking purposes. - """ - NZ,NY,NX = gridPtsPerEdge - OZ,OY,OX = boundingBoxOrigin - spacing1 = 1./gridSpacing - - zTemp,yTemp,xTemp = leftHandedCoord - i,j,k = round((xTemp-OX)*spacing1), round((yTemp-OY)*spacing1), round((zTemp-OZ)*spacing1) - octahedronIndices = [] - zippedNumbers = [(0,0,1),(0,0,-1),(0,1,0),(0,-1,0),(1,0,0),(-1,0,0)] - for x,y,z in zippedNumbers: - newI, newJ, newK = i + x, j + y, k + z - if newI < 0 or newI > (NX-1) or newJ < 0 or newJ > (NY-1) or newK < 0 or newK > (NZ-1): - continue - desiredPointIndex = int(round(newK*NX*NY + newJ*NX + newI)) - octahedronIndices.append(desiredPointIndex) - return octahedronIndices - def f_ray_intersect_polyhedron(pRayStartPos, pRayEndPos, faces,vertices, pTruncateToSegment): """This function returns TRUE if a ray intersects a triangle. It also calculates and returns the UV coordinates of said colision as part of the intersection test, @@ -341,18 +293,6 @@ def f_ray_intersect_polyhedron(pRayStartPos, pRayEndPos, faces,vertices, pTrunca #vBackfaceFinal = vBackfaces[distances.index(min(distances))] return vHitCount,vBackfaceFinal -def _findTriangleCenter(*args): - """ - Given three coordinates, it calculates the center of the triangle by averaging - the coordinates in each axes. - """ - p1,p2,p3 = args[0] - x = zip(p1,p2,p3) - result = [] - for elem in x: - result.append(sum(elem)/len(elem)) - return result - def findPointsCenter(*args): # Average down the column, such that we're averaging across all measurements in one dimension center = numpy.mean(args[0], axis = 0) @@ -365,99 +305,7 @@ def squaredTwoPointDistance(point1,point2): dist = sum(distsSquare) return dist -def projectPolyhedronToGrid(gridSpacing): - startTime = time() - from math import ceil - # Get the object, and decompose it - object_target = helper.getCurrentSelection()[0] - faces,verticesLH,vnormals,faceNormals = helper.DecomposeMesh(object_target,edit=False,copy=False,tri=True,transform=True,fn=True) - # Setup a bounding box for our polyhedron - b = setupBBOX(object_target) - corners = b.getBoundingBox() - - # Create grid points to fill in our bbox with grid points at regular distances - points, gridPtsPerEdge = makeGrid(corners[0],corners[1],gridSpacing) # These are left-handed - - # Create a list of all our gridPoint objects - gridPoints = [] - i = 0 - for point in points: - gridPoints.append(gridPoint(i,point,isPolyhedron = False)) - i += 1 - assert len(gridPoints) == len(points) - - # Create grid points to fill in our bbox with grid points at regular distances - points, gridPtsPerEdge = makeGrid(corners[0],corners[1],gridSpacing) # These are left-handed - - # For every face in the polyhedron - allCoordinates = [] - for face in faces: - # Dist is a small function that calculates the distance between two points. - dist = lambda x,y: vlen([y[i] - x[i] for i in range(len(x))]) - # Get the vertex coordinates and conver to numpy arrays... - triCoords = [numpy.array(verticesLH[i]) for i in face] - allCoordinates.extend(triCoords) - # ...use them to define the u/v vectors - pos = triCoords[0] - u = triCoords[1] - pos - v = triCoords[2] - pos - # SOmetimes the hypotenuse isn't fully represented. To remedy this, we will use an addition w vector - w = triCoords[2] - triCoords[1] - - # If either u or v is greater than the grid spacing, then we need to subdivide it - # We will use ceil: if we have a u of length 16, and grid spacing of 5, then we want - # a u at 0, 5, 10, 15 which is [0, 1, 2, 3] * gridSpacing. We need ceil to make this happen. - - # This is a much more efficient solution than the previous method of artificially increasing - # the densith of the polygon mesh by 4, then mapping each of those small points to the grid. - # However, it is also leaks more, for (as of yet) unknown reasons. - - # Minimum is one because range(1) gives us [0] - uSubunits, vSubunits, wSubunits = 1, 1, 1 - if vlen(u) > gridSpacing: - uSubunits = ceil(vlen(u)/gridSpacing) - if vlen(v) > gridSpacing: - vSubunits = ceil(vlen(v)/gridSpacing) - if vlen(w) > gridSpacing: - wSubunits = ceil(vlen(w)/gridSpacing) - # Because we have observed leakage, maybe we want to try trying a denser interpolation, using numpy's linspace? - for uSub in range(uSubunits): - percentU = uSub * gridSpacing / vlen(u) - assert percentU < 1.0 # Make sure that we have not stepped outside of our original u vector - # h represents the height of the hypotenuse at this u. We cannot go past the hypotenuse, so this will be - # our upper bound. - h = percentU * u + (1 - percentU) * v - for vSub in range(vSubunits): - percentV = vSub * gridSpacing / vlen(v) - assert percentV < 1.0 # Make sure that we have not stepped oustide of our original v vector. - interpolatedPoint = percentU * u + percentV * v - # The original if: statement asks if the distance from the origin to the interpolated point is less than - # the distance from the origin to the hypotenuse point. - # if vlen(interpolatedPoint) < vlen(h): - # Wouldn't it be a better idea to measure distance to the u position instead? - if (vlen(interpolatedPoint - percentU * u) < vlen(h - percentU * u)): - allCoordinates.append(interpolatedPoint + pos) - else: - # Used to be a break, changed to continue (not 100% sure if it made too much difference) - continue - # Because the above only interpolates the face, and may not completely capture the hypotenuse, let's separately - # interpolate points on the hypotenuse (the w vector) - for wSub in range(wSubunits): - percentW = wSub * gridSpacing / vlen(w) - if percentW > 1.0: - percentW = 1.0 - interpolatedPoint = percentW * w - allCoordinates.append(interpolatedPoint + triCoords[1]) - visualizeAtomArray('fineCoords',allCoordinates,5) - projectedIndices = set() - for coord in allCoordinates: - projectedPointIndex = getPointFrom3D(coord[::-1],gridSpacing,gridPtsPerEdge,b) - projectedIndices.add(projectedPointIndex) - visualizeAtomArray('projectedCoords',[points[i] for i in projectedIndices],5) - print('Projecting polyhedron to grid took ' + str(time() - startTime) + ' seconds.') - return(projectedIndices) - -def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): +def projectPolyhedronToGrid2(gridSpacing, radius = None): """ Takes a polyhedron, and builds a grid. In this grid: - Projects the polyhedron to the grid. @@ -485,7 +333,8 @@ def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): # Create grid points to fill in our bbox with grid points at regular distances points, gridPtsPerEdge = makeGrid(corners[0],corners[1],gridSpacing) # These are left-handed - # Create a list of all our gridPoint objects + # Create a gridPoint object for every single point we have in our grid. This pre-allocates all memory, so that future + # computations are faster. gridPoints = [] i = 0 for point in points: @@ -493,18 +342,19 @@ def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): i += 1 assert len(gridPoints) == len(points) + # Make a precomputed distance cube, where distanceCube,distX,distY,distZ = makeMarchingCube(gridSpacing,radius) - # Flatten and combine these arrays, so that it's easier to iterate over + # Flatten and combine these arrays. This makes it easier to iterate over. distanceCubeF,distXF,distYF,distZF = distanceCube.flatten(),distX.flatten(),distY.flatten(),distZ.flatten() zippedNumbers = zip(distanceCubeF,distXF,distYF,distZF) NZ,NY,NX = gridPtsPerEdge OZ, OY, OX = b.getBoundingBox()[0] - spacing1 = 1./gridSpacing - # For every face in the polyhedron - allCoordinates = [] + spacing1 = 1./gridSpacing # Inverse of the spacing. We compute this here, so we don't have to recompute it repeatedly + allCoordinates = [] # Tracker for all the fine coordiantes that we have interpolated for the faces of the polyhedron pointsToTestInsideOutside = set() - # Walk through the faces + # Walk through the faces, projecting each to the grid and marking immediate neighbors so we can test said + # neighbors for inside/outside later. for face in faces: # Dist is a small function that calculates the distance between two points. # dist = lambda x,y: vlen([y[i] - x[i] for i in range(len(x))]) @@ -512,11 +362,12 @@ def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): triCoords = [numpy.array(verticesLH[i]) for i in face] thisFaceFineCoords = list(triCoords) allCoordinates.extend(triCoords) - # ...use them to define the u/v vectors + # We will use these u/v vectors to interpolate points that reside on the face pos = triCoords[0] u = triCoords[1] - pos v = triCoords[2] - pos - # SOmetimes the hypotenuse isn't fully represented. To remedy this, we will use an addition w vector + # Smetimes the hypotenuse isn't fully represented. To remedy this, we will use an additional w vector + # to interpolate points on the hypotenuse w = triCoords[2] - triCoords[1] # If either u or v is greater than the grid spacing, then we need to subdivide it @@ -531,7 +382,8 @@ def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): # temporarily use a somewhat denser gridspacing to interpolate, and then project back onto our original spacing. # We'll decrease the gridspacing by 25% (so that it's 75% of the original). This seems to fix the issue without # any appreciable increase in computation time. - gridSpacingTempFine = gridSpacing * 3 / 4 + gridSpacingTempFine = gridSpacing / 3 + # Determine the number of grid spacing-sized points we can fit on each vector. # Minimum is one because range(1) gives us [0] uSubunits, vSubunits, wSubunits = 1, 1, 1 if vlen(u) > gridSpacingTempFine: @@ -544,7 +396,7 @@ def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): for uSub in range(uSubunits): percentU = uSub * gridSpacingTempFine / vlen(u) percentU = min(percentU, 1.0) # Make sure that we have not stepped outside of our original u vector - # h represents the height of the hypotenuse at this u. We cannot go past the hypotenuse, so this will be + # h represents the height of the hypotenuse at this u. Naturally, we cannot go past the hypotenuse, so this will be # our upper bound. h = percentU * u + (1 - percentU) * v for vSub in range(vSubunits): @@ -552,29 +404,35 @@ def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): percentV = min(percentV, 1.0) # Make sure that we have not stepped oustide of our original v vector. interpolatedPoint = percentU * u + percentV * v # The original if: statement asks if the distance from the origin to the interpolated point is less than - # the distance from the origin to the hypotenuse point. + # the distance from the origin to the hypotenuse point, as such: # if vlen(interpolatedPoint) < vlen(h): - # Wouldn't it be a better idea to measure distance to the u position instead? + # Wouldn't it be a better idea to measure distance to the u position instead? This is implemented below. if (vlen(interpolatedPoint - percentU * u) < vlen(h - percentU * u)): allCoordinates.append(interpolatedPoint + pos) thisFaceFineCoords.append(interpolatedPoint + pos) else: - # Used to be a break, changed to continue (not 100% sure if it made too much difference) break # Because the above only interpolates the face, and may not completely capture the hypotenuse, let's separately - # interpolate points on the hypotenuse (the w vector) + # interpolate points on the hypotenuse (the w vector). This way we can be 100% sure that our final solution is + # completely leak-proof for wSub in range(wSubunits): + # Apply the same proceudre we did above for u/v, just for w (for hypotenuse interpolation) percentW = wSub * gridSpacingTempFine / vlen(w) percentW = min(percentW, 1.0) interpolatedPoint = percentW * w allCoordinates.append(interpolatedPoint + triCoords[1]) thisFaceFineCoords.append(interpolatedPoint + triCoords[1]) projectedIndices = set() + # Once we have interpolated the face, let's project each interpolated point to the grid. for coord in thisFaceFineCoords: projectedPointIndex = getPointFrom3D(coord[::-1],gridSpacing,gridPtsPerEdge,b) projectedIndices.add(projectedPointIndex) + # Walk through each grid point that our face spans, gather its closest neighbors, and annotate the neighbors with + # minimum distance and closest faces, and flag them for testing inside/outside later. for P in list(projectedIndices): + # if P == 53144: + # print('We have projected to the point that causes issues!!!!!!!!!!!') # Get the point object corresponding to the index, and set its polyhedron attribute to true g = gridPoints[P] g.representsPolyhedron = True @@ -589,7 +447,11 @@ def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): continue # Get the point index that this coordinate corresponds to. desiredPointIndex = int(round(newK*NX*NY + newJ*NX + newI)) + # if desiredPointIndex == 54199: + # print('HIT with face ' + str(face) + ' and polygon point ' + str(P)) desiredPoint = gridPoints[desiredPointIndex] + if desiredPoint.representsPolyhedron == True: + continue # Add the current face to the its list of closest faces if face not in desiredPoint.closeFaces: desiredPoint.closeFaces.append(face) @@ -597,19 +459,19 @@ def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): desiredPoint.allDistances.append((v,d)) if d < desiredPoint.minDistance: desiredPoint.minDistance = d - desiredPoint.closestFaceIndex = len(desiredPoint.closeFaces) - 1 # Later down the road, we want to test as few points as possible for inside/outside. Therefore, # we will only test points that are - if abs(x) <= 1 and abs(y) <= 1 and abs(z) <= 1: - pointsToTestInsideOutside.add(desiredPointIndex) - visualizeAtomArray('fineCoords',allCoordinates,5) + # if abs(x) <= 1 and abs(y) <= 1 and abs(z) <= 1: + # pointsToTestInsideOutside.add(desiredPointIndex) + # visualizeAtomArray('fineCoords',allCoordinates,5) projectedIndices = [x.index for x in gridPoints if x.representsPolyhedron == True] # projectedIndices = set() # for coord in allCoordinates: # projectedPointIndex = getPointFrom3D(coord[::-1],gridSpacing,gridPtsPerEdge,b) # projectedIndices.add(projectedPointIndex) # visualizeAtomArray('projectedCoords',[points[i] for i in projectedIndices],5) - print('Projecting polyhedron to grid took ' + str(time() - startTime) + ' seconds.') + timeFinishProjection = time() + print('Projecting polyhedron to grid took ' + str(timeFinishProjection - startTime) + ' seconds.') # Let's start filling in inside outside. Here's the general algorithm: # Walk through all the points in our grid. Once we encounter a point that has closest faces, @@ -628,7 +490,7 @@ def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): for g in gridPoints: # Check if we've started a new line. If so, then we reset everything. # This test should precede all other test, because we don't want old knowldge - # to carry over to the new line. + # to carry over to the new line, since we don't know if the polygon is only partially encapsulated by the bounding box. if g.index > 0: # We can't check the first element, so we can skip it. coordDiff = g.globalCoord - gridPoints[g.index - 1].globalCoord coordDiffNonzero = [x != 0 for x in coordDiff] @@ -668,16 +530,16 @@ def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): numHits, thisBackFace = f_ray_intersect_polyhedron(g.globalCoord,endPoint,g.closeFaces,verticesLH,False) # We can check the face as well if we want to be super precise. If they dont' agree, we then check against the entire polyhedron. - if superFine == True: - if len(g.closeFaces) > 1: - uniquePoints2 = [] - [uniquePoints2.append(x) for x in g.closeFaces[1] if x not in uniquePoints2] - uniquePointsCoords2 = verticesLH[uniquePoints2] - endPoint2 = findPointsCenter(uniquePointsCoords2) - numHits2, thisBackFace2 = f_ray_intersect_polyhedron(g.globalCoord,endPoint2,g.closeFaces,verticesLH,False) - if len(g.closeFaces) == 1 or thisBackFace != thisBackFace2: - mismatchCounter += 1 - numHits, thisBackFace = f_ray_intersect_polyhedron(g.globalCoord,numpy.array([0.0,0.0,0.0]),faces,verticesLH,False) + # if superFine == True: + # if len(g.closeFaces) > 1: + # uniquePoints2 = [] + # [uniquePoints2.append(x) for x in g.closeFaces[1] if x not in uniquePoints2] + # uniquePointsCoords2 = verticesLH[uniquePoints2] + # endPoint2 = findPointsCenter(uniquePointsCoords2) + # numHits2, thisBackFace2 = f_ray_intersect_polyhedron(g.globalCoord,endPoint2,g.closeFaces,verticesLH,False) + # if len(g.closeFaces) == 1 or thisBackFace != thisBackFace2: + # mismatchCounter += 1 + # numHits, thisBackFace = f_ray_intersect_polyhedron(g.globalCoord,numpy.array([0.0,0.0,0.0]),faces,verticesLH,False) # Fill in inside outside attribute for this point, as pRayStartPos, pRayEndPos, faces, vertices, pTruncateToSegmentll as for any points before it g.isOutside = not thisBackFace @@ -695,7 +557,7 @@ def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): else: if g.isOutside == None: g.isOutside = True - + print('Flood filling grid inside/outside took ' + str(time() - timeFinishProjection) + ' seconds.') insidePoints = [g.index for g in gridPoints if g.isOutside == False] outsidePoints = [g.index for g in gridPoints if g.isOutside == True] surfacePoints = [g.index for g in gridPoints if g.representsPolyhedron == True] @@ -703,13 +565,16 @@ def projectPolyhedronToGrid2(gridSpacing, radius = None, superFine = False): visualizeAtomArray('insidePoints',[points[i] for i in insidePoints],0.5) visualizeAtomArray('outsidePoints',[points[i] for i in outsidePoints],0.5) visualizeAtomArray('surfacePoints',[points[i] for i in surfacePoints],0.5) - print('Superfine was ' + str(superFine) + ' and there were ' + str(mismatchCounter) + ' mismatches.') + # if superFine == True: + # print('Superfine was ' + str(superFine) + ' and there were ' + str(mismatchCounter) + ' mismatches.') print('Grid construction took ' + str(time() - startTime) + ' seconds for ' + str(len(faces)) + ' faces and ' + str(len(gridPoints)) + ' points.') - # visualizeAtomArray('badPoint', [gridPoints[insidePoints[14067]].globalCoord]) - # visualizeAtomArray('badPointEndPoint',[gridPoints[insidePoints[14067]].testedEndpoint]) - # visualizeAtomArray('badPointCloseFaces',verticesLH[gridPoints[insidePoints[14067]].closeFaces]) + # print('The bad point is at global index ' + str(insidePoints[15090])) + # visualizeAtomArray('badPoint', [gridPoints[insidePoints[15090]].globalCoord]) + # visualizeAtomArray('badPointEndPoint',[gridPoints[insidePoints[15090]].testedEndpoint]) + # visualizeAtomArray('badPointCloseFaces',verticesLH[gridPoints[insidePoints[15090]].closeFaces]) + # visualizeAtomArray('badPointMasterPoint', [gridPoints[53144].globalCoord]) return gridPoints -projectPolyhedronToGrid2(5, superFine = False) +projectPolyhedronToGrid2(5) From 9462241d873b46e091b6b72f9b26837befb5c6df Mon Sep 17 00:00:00 2001 From: wukevin Date: Thu, 30 Oct 2014 13:19:22 -0700 Subject: [PATCH 4/4] Addition of the inside/outside testing algorithm using the marching cube. --- autopack/Compartment.py | 588 ++++++++++++++++++++++++++++++++++++++++ autopack/Environment.py | 6 +- 2 files changed, 592 insertions(+), 2 deletions(-) diff --git a/autopack/Compartment.py b/autopack/Compartment.py index 16f8268..261ea3a 100644 --- a/autopack/Compartment.py +++ b/autopack/Compartment.py @@ -930,6 +930,8 @@ def BuildGrid(self, env): a, b = self.BuildGrid_jordan(env,ray=3) elif env.innerGridMethod == "pyray" and self.isOrthogonalBoudingBox!=1: # surfaces and interiors will be subtracted from it as normal! a, b = self.BuildGrid_pyray(env) + elif env.innerGridMethod == 'kevin' and self.isOrthogonalBoudingBox!=1: + a, b = self.BuildGrid_kevin(env) return a,b def prepare_buildgrid_box(self,env): @@ -981,9 +983,595 @@ def BuildGrid_box(self, env,vSurfaceArea): self.computeVolumeAndSetNbMol(env, self.surfacePoints, self.insidePoints,areas=vSurfaceArea) + def BuildGrid_kevin(self, env, superFine = False, visualize = False): + """ + Takes a polyhedron, and builds a grid. In this grid: + - Projects the polyhedron to the grid. + - Determines which points are inside/outside the polyhedron + - Determines point's distance to the polyhedron. + Usage: + Select desired object, run this program with desired grid spacing resolution. By default, + radius is 2 * gridSpacing. superFine provides the option doing a super leakproof test when + determining which points are inside or outside. This usually not necessary, because the + built-in algorithm has no known leakage cases. It is simply there as a safeguard. + """ + + def f_ray_intersect_polyhedron(pRayStartPos, pRayEndPos, faces,vertices, pTruncateToSegment): + """This function returns TRUE if a ray intersects a triangle. + It also calculates and returns the UV coordinates of said colision as part of the intersection test, + Makes sure that we are working with arrays + + * TAKES IN pRayStartPos as LEFT-HANDED COORDINATE (Z,Y,X) + * TAKES IN pRayEndPos AS LEFT-HANDED COORDINATE (Z,Y,X) + * faces is a list of vertex indices that defines the polyhedron + * vertices is the list of global coordinates that faces refers to + * pTruncateToSegment decides whether or not the segment terminates at the end position, or keeps on + going forever + """ + pRayStartPos = numpy.array(pRayStartPos) + pRayEndPos = numpy.array(pRayEndPos) + + vLineSlope = pRayEndPos - pRayStartPos #This line segment defines an infinite line to test for intersection + vPolyhedronPos = numpy.array((0,0,0)) + vTriPoints = vertices + vTriPolys = faces + vBackface = None + vBackfaceFinal = None + maxDistance = 99999 + vHitCount = 0 + + # Globalize the polyhedron + # Present in original coffee script, but unnecessary because DecomposeMesh gives global coords + #for i in range(len(vTriPoints)): + # vTriPoints[i] = vTriPoints[i] + vPolyhedronPos + + vEpsilon = 0.00001 + vBreakj = False + vCollidePos = None + + # Walk through each polygon in a polyhedron + for testingFace in vTriPolys: + # Loop through all the polygons in an input polyhedron + #vQuadrangle = 1 # Default says polygon is a quadrangle. + #vLoopLimit = 2 # Default k will loop through polygon assuming its a quad. + #if (vTriPolys[j+3] == vTriPolys[j+2]) # Test to see if quad is actually just a triangle. + # { + vQuadrangle = 0 # Current polygon is not a quad, it's a triangle. + vLoopLimit = 1 # Set k loop to only cycle one time. + + for k in range(vLoopLimit): + vTriPt0 = numpy.array(vTriPoints[testingFace[0]]) # Always get the first point of a quad/tri + vTriPt1 = numpy.array(vTriPoints[testingFace[1]]) # Get point 1 for a tri and a quad's first pass, but skip for a quad's second pass + vTriPt2 = numpy.array(vTriPoints[testingFace[2]]) # Get point 2 for a tri and a quad's first pass, but get point 3 only for a quad on its second pass. + + vE1 = vTriPt1 - vTriPt0 # Get the first edge as a vector. + vE2 = vTriPt2 - vTriPt0 # Get the second edge. + h = vcross(vLineSlope, vE2) + + a = f_dot_product(vE1,h) # Get the projection of h onto vE1. + if a > -vEpsilon and a < vEpsilon: + continue# If the ray is parallel to the plane then it does not intersect it, i.e, a = 0 +/- given rounding slope. + # If the polygon is a quadrangle, test the other triangle that comprises it. + + F = 1.0/a + s = pRayStartPos - vTriPt0 # Get the vector from the origin of the triangle to the ray's origin. + u = F * f_dot_product(s,h) + if u < 0.0 or u > 1.0: + continue + #/* Break if its outside of the triangle, but try the other triangle if in a quad. + #U is described as u = : start of vE1 = 0.0, to the end of vE1 = 1.0 as a percentage. + #If the value of the U coordinate is outside the range of values inside the triangle, + #then the ray has intersected the plane outside the triangle.*/ + + q = vcross(s, vE1) + v = F * f_dot_product(vLineSlope,q) + if v <0.0 or u+v > 1.0: + continue + #/* Break if outside of the triangles v range. + #If the value of the V coordinate is outside the range of values inside the triangle, + #then the ray has intersected the plane outside the triangle. + #U + V cannot exceed 1.0 or the point is not in the triangle. + #If you imagine the triangle as half a square this makes sense. U=1 V=1 would be in the + #lower left hand corner which would be in the second triangle making up the square.*/ + + vCollidePos = vTriPt0 + u*vE1 + v*vE2 # This is the global collision position. + assert len(vCollidePos) == 3 + + # The ray is hitting a triangle, now test to see if its a triangle hit by the ray. + vBackface = False + if f_dot_product(vLineSlope, vCollidePos - pRayStartPos) > 0: # This truncates our infinite line to a ray pointing from start THROUGH end positions. + vHitCount += 1 + #print(testingFace) + if pTruncateToSegment and vlen(vLineSlope) < vlen(vCollidePos - pRayStartPos): + print('broken') + break # This truncates our ray to a line segment from start to end positions. + + d = squaredTwoPointDistance(pRayStartPos,vCollidePos) + if d >= maxDistance: + continue + if a < vEpsilon: # Test to see if the triangle hit is a backface. + #set master grid to organelle->getname inside + vBackface = True + # This stuff is specific to our Point inside goals. + vBreakj = True # To see if a point is inside, I can stop at the first backface hit. + #break + else: + vBreakj = True + vBackfaceFinal = vBackface + maxDistance = d + #vBackfaceFinal = vBackfaces[distances.index(min(distances))] + return vHitCount,vBackfaceFinal + + def makeGrid(lowerLeft,upperRight,spacing, inclusive = True): + """ + Given lower left back corner, an upper right front corner, and a grid spacing, construct a list of + coordinates that representsly an evenly spaced grid spanning the entire volume. Returned coordinates + preserve the right/lefthandedness of the input coordinates, making this function work equally well + in both 3D coordinate systems. + """ + pointCoords = [] + xSize,ySize,zSize = upperRight[0] - lowerLeft[0], upperRight[1] - lowerLeft[1], upperRight[2] - lowerLeft[2] + if inclusive: + nx,ny,nz = xSize / spacing + 1, ySize / spacing + 1, zSize / spacing + 1 + else: + nx,ny,nz = xSize / spacing, ySize / spacing, zSize / spacing + nx,ny,nz = round(nx),round(ny),round(nz) + x = numpy.linspace(lowerLeft[0],upperRight[0],nx) + y = numpy.linspace(lowerLeft[1],upperRight[1],ny) + z = numpy.linspace(lowerLeft[2],upperRight[2],nz) + for a in x: + for b in y: + for c in z: + pointCoords.append((a,b,c)) + # Order of above had to change to correctly produce a viable BBOX + return pointCoords, (nx,ny,nz) + + def makeMarchingCube(gridSpacing,r): + """ + Create a numpy array that represents the precomputed distances to each point + for the cube of points surrounding our center point. + """ + def _pythagorean(*edgeLengths): + from math import sqrt + total = 0 + for length in edgeLengths: + total += length * length + distance = sqrt(total) + return distance + from math import ceil + pointsForRadius = ceil(r/gridSpacing) # Number of grid points required to represent our radius, rounded up + pointsInEdge = 2 * pointsForRadius + 1 # Number of points in one edge of our cube + + center = pointsForRadius # The index if our center point + cube = numpy.zeros(shape=(pointsInEdge,pointsInEdge,pointsInEdge)) + distX = numpy.zeros(shape=(pointsInEdge,pointsInEdge,pointsInEdge)) + distY = numpy.zeros(shape=(pointsInEdge,pointsInEdge,pointsInEdge)) + distZ = numpy.zeros(shape=(pointsInEdge,pointsInEdge,pointsInEdge)) + for a in range(pointsInEdge): + lenX = a - center + for b in range(pointsInEdge): + lenY = b - center + for c in range(pointsInEdge): + lenZ = c - center + cube[a][b][c] = _pythagorean(lenX,lenY,lenZ) * gridSpacing + distX[a][b][c] = lenX + distY[a][b][c] = lenY + distZ[a][b][c] = lenZ + return cube, distX, distY, distZ + + def getPointFrom3D(pt3d, gridSpacing, gridPtsPerEdge, bb): + """ + Starts from orthogonal grid. Has a lower left, upper right coordinate from bounding box. + Given lower left and known spacing of grid, we can provide any vector in 3D space and it will + algebraically return index of closest grid point. + + Given a RIGHT-HANDED fine coordinate, it will return the index of the closest point in the LEFT-HANDED coarse grid. + """ + x, y, z = pt3d # Continuous 3D point to be discretized + spacing1 = 1./gridSpacing # Grid spacing = diagonal of the voxel determined by smalled packing radius + NZ, NY, NX = gridPtsPerEdge # vector = [length, height, depth] of grid, units = gridPoints + + # Added to fix some weird edge behavior + #NZ, NY, NX = round(NZ),round(NY),round(NX) + #print(NX,NY,NZ) + OZ, OY, OX = bb[0] # origin of Pack grid. OX and OZ were switched to compensate for coordinate systems + #print(OX,OY,OZ) + + # Algebra gives nearest gridPoint ID to pt3D + i = min( NX-1, max( 0, round((x-OX)*spacing1))) + j = min( NY-1, max( 0, round((y-OY)*spacing1))) + k = min( NZ-1, max( 0, round((z-OZ)*spacing1))) + result = k*NX*NY + j*NX + i + return int(result) + + def squaredTwoPointDistance(point1,point2): + """Computes and returns the distance between two points in 3D space""" + dists = [numpy.float64(point1[i] - point2[i]) for i in range(len(point1))] + distsSquare = [x**2 for x in dists] + dist = sum(distsSquare) + return dist + + def findPointsCenter(*args): + # Average down the column, such that we're averaging across all measurements in one dimension + center = numpy.mean(args[0], axis = 0) + return center + + def vcross(v1,v2): + """ + Returns the cross-product of two 3D vectors. + """ + x1,y1,z1 = v1 + x2,y2,z2 = v2 + return (y1*z2-y2*z1, z1*x2-z2*x1, x1*y2-x2*y1) + + def vlen(v1): + """ + Returns the length of a 3D vector. + """ + x1,y1,z1 = v1 + return sqrt(x1*x1 + y1*y1 + z1*z1) + + def f_dot_product(vector1,vector2): + """ + Return the dot product of two 3D vectors. + """ + dottedVectors = [vector1[i] * vector2[i] for i in range(len(vector1))] + return sum(dottedVectors) + + def visualizeAtomArray(name = 'Test', vertexCoordinates = [(0,0,0)], sphereRadius = 5): + """ + Given coordinates, visualize them in Cinema 4D using an Atom Array Object. + """ + atom = c4d.BaseObject(c4d.Oatomarray) + atom.SetName(name) + atom[1000] = 3 #radius cylinder + atom[1001] = sphereRadius #radius sphere + atom[1002] = 5 #subdivision + helper.addObjectToScene(helper.getCurrentScene(),atom,parent=None) + + parent = atom + coords = vertexCoordinates + nface = 0 + obj= c4d.PolygonObject(len(coords), nface) + obj.SetName(name + '_child') + cd4vertices = map(helper.FromVec,coords) + map(obj.SetPoint,range(len(coords)),cd4vertices) + helper.addObjectToScene(helper.getCurrentScene(),obj,parent=parent) + + class gridPoint: + def __init__(self,i,globalC,isPolyhedron): + self.index = int(i) + self.isOutside = None + self.minDistance = 99999 # Only store a number here if within certain distance from polyhedron + self.representsPolyhedron = isPolyhedron + self.closeFaces = [] + self.closestFaceIndex = 0 + self.testedEndpoint = None + self.allDistances = [] # Stores a tuple list of distances to all points. (point,distance) = (5,2.5) + self.globalCoord = numpy.array(globalC) # Stores the global coordinate associated with this point + + if self.ghost : return + from time import time + # Copied from jordan + t0 = t1 = time() + if self.isBox : + self.overwriteSurfacePts = True + if self.overwriteSurfacePts: + self.ogsurfacePoints = self.vertices[:] + self.ogsurfacePointsNormals = self.vnormals[:] + else : + self.createSurfacePoints(maxl=env.grid.gridSpacing) + + # Inferred from reading Compartment.py + gridSpacing = env.grid.gridSpacing + + radius = gridSpacing + + startTime = time() + from math import ceil + # Get the object, and decompose it + # object_target = helper.getCurrentSelection()[0] + # faces,verticesLH,vnormals,faceNormals = helper.DecomposeMesh(object_target,edit=False,copy=False,tri=True,transform=True,fn=True) + faces = self.faces + verticesLH = self.vertices + vnormals = self.vnormals + # Setup a bounding box for our polyhedron + # b = setupBBOX(object_target) + # corners = b.getBoundingBox() + corners = self.bb + + # Create grid points to fill in our bbox with grid points at regular distances + points, gridPtsPerEdge = makeGrid(corners[0],corners[1],gridSpacing) # These are left-handed + + # Create a gridPoint object for every single point we have in our grid. This pre-allocates all memory, so that future + # computations are faster. + gridPoints = [] + i = 0 + for point in points: + gridPoints.append(gridPoint(i,point,isPolyhedron = False)) + i += 1 + assert len(gridPoints) == len(points) + + # Make a precomputed distance cube, where + distanceCube,distX,distY,distZ = makeMarchingCube(gridSpacing,radius) + # Flatten and combine these arrays. This makes it easier to iterate over. + distanceCubeF,distXF,distYF,distZF = distanceCube.flatten(),distX.flatten(),distY.flatten(),distZ.flatten() + zippedNumbers = zip(distanceCubeF,distXF,distYF,distZF) + + NZ,NY,NX = gridPtsPerEdge + OZ, OY, OX = self.bb[0] + spacing1 = 1./gridSpacing # Inverse of the spacing. We compute this here, so we don't have to recompute it repeatedly + allCoordinates = [] # Tracker for all the fine coordiantes that we have interpolated for the faces of the polyhedron + # pointsToTestInsideOutside = set() + # Walk through the faces, projecting each to the grid and marking immediate neighbors so we can test said + # neighbors for inside/outside later. + for face in faces: + # Dist is a small function that calculates the distance between two points. + # dist = lambda x,y: vlen([y[i] - x[i] for i in range(len(x))]) + # Get the vertex coordinates and conver to numpy arrays... + triCoords = [numpy.array(verticesLH[i]) for i in face] + thisFaceFineCoords = list(triCoords) + allCoordinates.extend(triCoords) + # We will use these u/v vectors to interpolate points that reside on the face + pos = triCoords[0] + u = triCoords[1] - pos + v = triCoords[2] - pos + # Smetimes the hypotenuse isn't fully represented. To remedy this, we will use an additional w vector + # to interpolate points on the hypotenuse + w = triCoords[2] - triCoords[1] + + # If either u or v is greater than the grid spacing, then we need to subdivide it + # We will use ceil: if we have a u of length 16, and grid spacing of 5, then we want + # a u at 0, 5, 10, 15 which is [0, 1, 2, 3] * gridSpacing. We need ceil to make this happen. + + # This is a much more efficient solution than the previous method of artificially increasing + # the densith of the polygon mesh by 4, then mapping each of those small points to the grid. + # However, it is also leaks more, for (as of yet) unknown reasons. + + # It was discovered that by using the default gridspacing, some faces will produce leakage. One solution is to + # temporarily use a somewhat denser gridspacing to interpolate, and then project back onto our original spacing. + # We'll decrease the gridspacing by 25% (so that it's 75% of the original). This seems to fix the issue without + # any appreciable increase in computation time. + gridSpacingTempFine = gridSpacing / 3 + # Determine the number of grid spacing-sized points we can fit on each vector. + # Minimum is one because range(1) gives us [0] + uSubunits, vSubunits, wSubunits = 1, 1, 1 + if vlen(u) > gridSpacingTempFine: + uSubunits = ceil(vlen(u)/gridSpacingTempFine) + 1 + if vlen(v) > gridSpacingTempFine: + vSubunits = ceil(vlen(v)/gridSpacingTempFine) + 1 + if vlen(w) > gridSpacingTempFine: + wSubunits = ceil(vlen(w)/gridSpacingTempFine) + 1 + # Because we have observed leakage, maybe we want to try trying a denser interpolation, using numpy's linspace? + for uSub in range(uSubunits): + percentU = uSub * gridSpacingTempFine / vlen(u) + percentU = min(percentU, 1.0) # Make sure that we have not stepped outside of our original u vector + # h represents the height of the hypotenuse at this u. Naturally, we cannot go past the hypotenuse, so this will be + # our upper bound. + h = percentU * u + (1 - percentU) * v + for vSub in range(vSubunits): + percentV = vSub * gridSpacingTempFine / vlen(v) + percentV = min(percentV, 1.0) # Make sure that we have not stepped oustide of our original v vector. + interpolatedPoint = percentU * u + percentV * v + # The original if: statement asks if the distance from the origin to the interpolated point is less than + # the distance from the origin to the hypotenuse point, as such: + # if vlen(interpolatedPoint) < vlen(h): + # Wouldn't it be a better idea to measure distance to the u position instead? This is implemented below. + if (vlen(interpolatedPoint - percentU * u) < vlen(h - percentU * u)): + allCoordinates.append(interpolatedPoint + pos) + thisFaceFineCoords.append(interpolatedPoint + pos) + else: + break + # Because the above only interpolates the face, and may not completely capture the hypotenuse, let's separately + # interpolate points on the hypotenuse (the w vector). This way we can be 100% sure that our final solution is + # completely leak-proof + for wSub in range(wSubunits): + # Apply the same proceudre we did above for u/v, just for w (for hypotenuse interpolation) + percentW = wSub * gridSpacingTempFine / vlen(w) + percentW = min(percentW, 1.0) + interpolatedPoint = percentW * w + allCoordinates.append(interpolatedPoint + triCoords[1]) + thisFaceFineCoords.append(interpolatedPoint + triCoords[1]) + projectedIndices = set() + # Once we have interpolated the face, let's project each interpolated point to the grid. + for coord in thisFaceFineCoords: + projectedPointIndex = getPointFrom3D(coord[::-1],gridSpacing,gridPtsPerEdge,self.bb) + projectedIndices.add(projectedPointIndex) + + # Walk through each grid point that our face spans, gather its closest neighbors, and annotate the neighbors with + # minimum distance and closest faces, and flag them for testing inside/outside later. + for P in list(projectedIndices): + # if P == 53144: + # print('We have projected to the point that causes issues!!!!!!!!!!!') + # Get the point object corresponding to the index, and set its polyhedron attribute to true + g = gridPoints[P] + g.representsPolyhedron = True + # Get the coordinates of the point, and convert them to grid units + zTemp,yTemp,xTemp = g.globalCoord + i,j,k = round((xTemp-OX)*spacing1), round((yTemp-OY)*spacing1), round((zTemp-OZ)*spacing1) + # Let's step through our distance cube, and assign faces/closest distances to each + for d,x,y,z in zippedNumbers: + # Get the grid indices for the point we're considering, and pass if we're stepping oustide the boundaries + newI, newJ, newK = i + x, j + y, k + z + if newI < 0 or newI > (NX-1) or newJ < 0 or newJ > (NY-1) or newK < 0 or newK > (NZ-1): + continue + # Get the point index that this coordinate corresponds to. + desiredPointIndex = int(round(newK*NX*NY + newJ*NX + newI)) + # if desiredPointIndex == 54199: + # print('HIT with face ' + str(face) + ' and polygon point ' + str(P)) + desiredPoint = gridPoints[desiredPointIndex] + if desiredPoint.representsPolyhedron == True: + continue + # Add the current face to the its list of closest faces + if face not in desiredPoint.closeFaces: + desiredPoint.closeFaces.append(face) + # Add the distance to the point's list of distances, and overwrite minimum distance if appropriate + desiredPoint.allDistances.append((v,d)) + if d < desiredPoint.minDistance: + desiredPoint.minDistance = d + # Later down the road, we want to test as few points as possible for inside/outside. Therefore, + # we will only test points that are + # if abs(x) <= 1 and abs(y) <= 1 and abs(z) <= 1: + # pointsToTestInsideOutside.add(desiredPointIndex) + # visualizeAtomArray('fineCoords',allCoordinates,5) + # projectedIndices = [x.index for x in gridPoints if x.representsPolyhedron == True] + # projectedIndices = set() + # for coord in allCoordinates: + # projectedPointIndex = getPointFrom3D(coord[::-1],gridSpacing,gridPtsPerEdge,b) + # projectedIndices.add(projectedPointIndex) + # visualizeAtomArray('projectedCoords',[points[i] for i in projectedIndices],5) + timeFinishProjection = time() + print('Projecting polyhedron to grid took ' + str(timeFinishProjection - startTime) + ' seconds.') + + # Let's start filling in inside outside. Here's the general algorithm: + # Walk through all the points in our grid. Once we encounter a point that has closest faces, + # then we know we need to test it for inside/outside. Once we test that for inside/outside, we + # fill in all previous points with that same inside outisde property. To account for the possible + # situation that there is a surface that is only partially bound by the bbox, then we need to + # reset the insideOutsideTracker every time we have a change in more than 1 of the 3 coordinates + # because that indicates we're starting a new row/column of points. + + # This tracks our inside/outside. + isOutsideTracker = None + # This tracks the points that we've iterated over which we do not know if inside/outside. Remember + # to reset every time we find an inside/outside. + emptyPointIndicies = [] + mismatchCounter = 0 + for g in gridPoints: + # Check if we've started a new line. If so, then we reset everything. + # This test should precede all other test, because we don't want old knowldge + # to carry over to the new line, since we don't know if the polygon is only partially encapsulated by the bounding box. + if g.index > 0: # We can't check the first element, so we can skip it. + coordDiff = g.globalCoord - gridPoints[g.index - 1].globalCoord + coordDiffNonzero = [x != 0 for x in coordDiff] + if sum(coordDiffNonzero) > 1: + # assert len(emptyPointIndicies) == 0 # When starting a new line, we shouldn't have any unknowns from the previous line + isOutsideTracker = None + emptyPointIndicies = [] + + # There's no point testing inside/outside for points that are on the surface. So if + # we get to one, we just skip over it. + if g.representsPolyhedron == True: + g.isOutside = None + continue + + if len(g.closeFaces) == 0: + # If it's not close to any faces, and we don't know if this row is inside/outside, then + # we have to wait till later to figure it out + if isOutsideTracker == None: + emptyPointIndicies.append(g.index) + # However, if we do know , we can just use the previous one to fill + else: + g.isOutside = isOutsideTracker + # If there are close faces attached to it, then we need to test it for inside/outside. + else: + # Find centroid of all the vertices of all the close faces. This will be our endpoint + # when casting a ray for collision testing. + uniquePoints = [] + # This takes just the first face and projects to the center of it. + # [uniquePoints.append(x) for x in g.closeFaces[0] if x not in uniquePoints] + [uniquePoints.append(x) for x in g.closeFaces[g.closestFaceIndex] if x not in uniquePoints] + uniquePointsCoords = verticesLH[uniquePoints] + endPoint = findPointsCenter(uniquePointsCoords) + g.testedEndpoint = endPoint + + # Draw a ray to that point, and see if we hit a backface or not + numHits, thisBackFace = f_ray_intersect_polyhedron(g.globalCoord,endPoint,g.closeFaces,verticesLH,False) + + # We can check the face as well if we want to be super precise. If they dont' agree, we then check against the entire polyhedron. + if superFine == True: + if len(g.closeFaces) > 1: + uniquePoints2 = [] + [uniquePoints2.append(x) for x in g.closeFaces[1] if x not in uniquePoints2] + uniquePointsCoords2 = verticesLH[uniquePoints2] + endPoint2 = findPointsCenter(uniquePointsCoords2) + numHits2, thisBackFace2 = f_ray_intersect_polyhedron(g.globalCoord,endPoint2,g.closeFaces,verticesLH,False) + if len(g.closeFaces) == 1 or thisBackFace != thisBackFace2: + mismatchCounter += 1 + numHits, thisBackFace = f_ray_intersect_polyhedron(g.globalCoord,numpy.array([0.0,0.0,0.0]),faces,verticesLH,False) + + # Fill in inside outside attribute for this point, as pRayStartPos, pRayEndPos, faces, vertices, pTruncateToSegmentll as for any points before it + g.isOutside = not thisBackFace + isOutsideTracker = not thisBackFace + for i in emptyPointIndicies: + gridPoints[i].isOutside = isOutsideTracker + # Because we have filled in all the unknowns, we can reset that counter. + emptyPointIndicies = [] + + # Final pass through for quality/sanity checks. + for g in gridPoints: + if g.representsPolyhedron == True: + assert g.isOutside == None + continue + else: + if g.isOutside == None: + g.isOutside = True + print('Flood filling grid inside/outside took ' + str(time() - timeFinishProjection) + ' seconds.') + insidePoints = [g.index for g in gridPoints if g.isOutside == False] + outsidePoints = [g.index for g in gridPoints if g.isOutside == True] + surfacePoints = [g.index for g in gridPoints if g.representsPolyhedron == True] + + if visualize == True: + visualizeAtomArray('insidePoints',[points[i] for i in insidePoints],0.5) + visualizeAtomArray('outsidePoints',[points[i] for i in outsidePoints],0.5) + visualizeAtomArray('surfacePoints',[points[i] for i in surfacePoints],0.5) + if superFine == True: + print('Superfine was ' + str(superFine) + ' and there were ' + str(mismatchCounter) + ' mismatches.') + print('Grid construction took ' + str(time() - startTime) + ' seconds for ' + str(len(faces)) + ' faces and ' + str(len(gridPoints)) + ' points.') + + # print('The bad point is at global index ' + str(insidePoints[15090])) + # visualizeAtomArray('badPoint', [gridPoints[insidePoints[15090]].globalCoord]) + # visualizeAtomArray('badPointEndPoint',[gridPoints[insidePoints[15090]].testedEndpoint]) + # visualizeAtomArray('badPointCloseFaces',verticesLH[gridPoints[insidePoints[15090]].closeFaces]) + # visualizeAtomArray('badPointMasterPoint', [gridPoints[53144].globalCoord]) + self.insidePoints = insidePoints + self.surfacePoints = surfacePoints + + # Pasted from BuildGrid_jordan + t1 = time() + nbGridPoints = len(env.grid.masterGridPositions) + + from scipy import spatial + srfPts = self.ogsurfacePoints + self.OGsrfPtsBht = bht = spatial.cKDTree(tuple(srfPts),leafsize=10) + surfPtsBB,surfPtsBBNorms = self.getSurfaceBB(srfPts,env) + srfPts = surfPtsBB + # I have no idea what idarray is..... + if autopack.verbose: + print ("compare length id distances",nbGridPoints,len(idarray),len(distances), + (nbGridPoints == len(idarray)), + (nbGridPoints == len(distances))) + ex = True#True if nbGridPoints == len(idarray) else False + surfacePoints,surfacePointsNormals = self.extendGridArrays(nbGridPoints, + srfPts,surfPtsBBNorms, + env,extended=ex) + + + insidePoints = insidePoints + if autopack.verbose: + print('time to extend arrays', time()-t1) + # print('Total time', time()-t0) + + self.insidePoints = insidePoints + self.surfacePoints = surfacePoints + self.surfacePointsCoords = surfPtsBB + self.surfacePointsNormals = surfacePointsNormals + if autopack.verbose: + print('%s surface pts, %d inside pts, %d tot grid pts, %d master grid'%( + len(self.surfacePoints), len(self.insidePoints), + nbGridPoints, len(env.grid.masterGridPositions))) + + normalList2,areas = self.getFaceNormals(self.vertices, faces,fillBB=env.fillBB) + vSurfaceArea = sum(areas) + self.computeVolumeAndSetNbMol(env, self.surfacePoints, + self.insidePoints,areas=vSurfaceArea) + return self.insidePoints, self.surfacePoints + + #Jordan Curve Theorem def BuildGrid_jordan(self, env,ray=1): """Build the compartment grid ie surface and inside point using jordan theorem and host raycast""" + # create surface points if self.ghost : return t0 = t1 = time() diff --git a/autopack/Environment.py b/autopack/Environment.py index 2263008..26f154a 100644 --- a/autopack/Environment.py +++ b/autopack/Environment.py @@ -1020,7 +1020,7 @@ def __init__(self,name="H"): "placeMethod": {"name":"placeMethod","value":"jitter","values":self.listPlaceMethod,"default":"placeMethod","type":"liste","description":" Overriding Packing Method = ","width":30}, "use_gradient":{"name":"use_gradient","value":False,"default":False,"type":"bool","description":"Use gradients if defined","width":150}, "gradients":{"name":"gradients","value":"","values":[],"default":"","type":"liste","description":"Gradients available","width":150}, - "innerGridMethod": {"name":"innerGridMethod","value":"bhtree","values":["bhtree","sdf","jordan","jordan3","pyray"],"default":"innerGridMethod","type":"liste","description":" Method to calculate the inner grid:","width":30}, + "innerGridMethod": {"name":"innerGridMethod","value":"bhtree","values":["bhtree","sdf","jordan","jordan3","pyray","kevin"],"default":"innerGridMethod","type":"liste","description":" Method to calculate the inner grid:","width":30}, "overwritePlaceMethod":{"name":"overwritePlaceMethod","value":False,"default":False,"type":"bool","description":"Overwrite per-ingredient packing method with Overriding Packing Method:","width":300}, "saveResult": {"name":"saveResult","value":False,"default":False,"type":"bool","description":"Save packing result to .apr file (enter full path below):","width":200}, "resultfile": {"name":"resultfile","value":"fillResult","default":"fillResult","type":"filename","description":"result filename","width":200}, @@ -2210,7 +2210,9 @@ def BuildGrids(self): elif self.innerGridMethod == "jordan3" and compartment.isOrthogonalBoudingBox!=1: # surfaces and interiors will be subtracted from it as normal! a, b = compartment.BuildGrid_jordan(self,ray=3) elif self.innerGridMethod == "pyray" and compartment.isOrthogonalBoudingBox!=1: # surfaces and interiors will be subtracted from it as normal! - a, b = compartment.BuildGrid_pyray(self) + a, b = compartment.BuildGrid_pyray(self) + elif self.innerGridMethod == "kevin" and compartment.isOrthogonalBoudingBox!=1: + a, b = compartment.BuildGrid_kevin(self) aInteriorGrids.append(a) print("I'm ruther in the loop") aSurfaceGrids.append(b)