Adjustable Mesh Plot
Background
A typical implicit plot algorithm uses a fine grid of points to render the plot. For many graphs the grid is required to be so fine that the computational complexity becomes a major performance issue. The adjustable mesh algorithm is designed to intelligently determine the areas of the grid that need refinement, and only increase the resolution in those areas. This allows the algorithm to render the full detail of the plot without a lot of wasted computation.
The Data Structure
The adjustable mesh data structure is composed of two parts; faces and edges.
An edge is two points defining a line segment, which in addition stores references to the faces on each side of the edge.
A face is composed of three edges. Two edges that share a common vertex, and a hypotenuse that shares one vertex with each edge.
The Algorithm
The algorithm begins with an initial grid of coarse resolution, and scans the grid for areas where it needs to be refined.
The algorithm considers any area where the function is changing rapidly to be an area in need of refinement. How much the function is changing starts by looking at two faces that share a common hypotenuse. Using points at opposite corners of the square and the midpoint of the hypotenuse the algorithm creates to 3 dimensional vectors where the z component is the value of f(x,y) at each point. The two vectors are created by subtracting midpoint from opposite corner, and subtracting corner from midpoint. It then looks at the angle between the two vectors. If the angle is very close to 0 then the function isn't changing much in that area, but if the angle is above a predetermined threshold then the function is considered to be changing rapidly.
If the function is changing rapidly at some point in the grid the algorithm then refines the grid by adding a point at the midpoint of the square. The point is added at the middle of the square to avoid tears occurring in adjacent areas of differing resolutions.
The algorithm continues to go through all of the points in the grid for a max level of recursion specified when the object is created.
The code
Below is our implementation of the adjustable mesh algorithm in Sage.
PlanarEdge
The PlanarEdge class represents one edge in the grid.
''' Planar Edge: Holds information of a directed edge of a planar simplex. ''' class PlanarEdge: def __init__(self, pt1, pt2, face1, face2): self.__endPoints = [ pt1, pt2 ] self.__faces = [ face1, face2 ] def __eq__(self, rhs): return ( self.__endPoints[0] == rhs.__endPoints[0] and self.__endPoints[1] == rhs.__endPoints[1] ) def __str__(self): return "{ " + str(self.__endPoints[0]) + " , " + str(self.__endPoints[1]) + " }" ''' Returns true if pt is contained in this edge ''' def __contains__(self, pt): return ( self.__endPoints[0] == pt or self.__endPoints[1] == pt ) def getVertex(self, i): return self.__endPoints[i] def setVertex(self, i, pt): self.__endPoints[i] = pt def getFace(self, i): return self.__faces[i] def setFace(self, i, face): self.__faces[i] = face ''' Returns the midpoint of the edge ''' def getMidpoint(self): return ( (self.__endPoints[0][0] + self.__endPoints[1][0])/2, (self.__endPoints[0][1] + self.__endPoints[1][1])/2 ) ''' Adds a line representing the edge to the graphics object g with a line color of color ''' def drawEdge(self, g, color): x0,x1 = ( self.__endPoints[0][0], self.__endPoints[1][0] ) y0,y1 = ( self.__endPoints[0][1], self.__endPoints[1][1] ) l = line( [ (x0,y0),(x1,y1) ], rgbcolor=color ) g += l return g
PlanarFace
The PlanarFace class represents on face in the grid.
''' PlanarFace: face of a planar simplex ''' class PlanarFace: def __init__(self, leg1, leg2, hyp): self.__sides = [ leg1, leg2, hyp ] self.__isProcessed = False def __str__(self): return "{ " + str(self.__sides[0]) + " , " + str(self.__sides[1]) + " , " + str(self.__sides[2]) + " }" def getLegs(self): return [ self.__sides[0] , self.__sides[1] ] def getLeg(self,i): return self.__sides[i] def getHypotenuse(self): return self.__sides[2] def getEdges(self): return self.__sides def setLegs(self, leg1, leg2): self.__sides[0] = leg1 self.__sides[1] = leg2 def setHypotenuse(self, hyp): self.__sides[2] = hyp def getCommonVertex(self): if self.__sides[0].getVertex(0) == self.__sides[1].getVertex(0) or self.__sides[0].getVertex(0) == self.__sides[1].getVertex(1): return self.__sides[0].getVertex( 0 ) else: return self.__sides[0].getVertex( 1 ) def getVerticies(self): return [ self.getCommonVertex() , self.__sides[2].getVertex( 0 ), self.__sides[2].getVertex( 1 ) ] def getOrderedSides(self): ''' Returns the 3 points of the two short sides with the first 2 points in side1 and the last 2 points in side2. ''' if self.__sides[0].getVertex( 0 ) == self.__sides[1].getVertex( 0 ): return [ self.__sides[0].getVertex(1),self.__sides[0].getVertex(0),self.__sides[1].getVertex(1) ] elif self.__sides[0].getVertex( 0 ) == self.__sides[1].getVertex( 1 ): return [ self.__sides[0].getVertex(1), self.__sides[0].getVertex(0), self.__sides[1].getVertex(0) ] elif self.__sides[0].getVertex( 1 ) == self.__sides[1].getVertex( 1 ): return [ self.__sides[0].getVertex(0), self.__sides[0].getVertex(1), self.__sides[1].getVertex(0) ] else: return [ self.__sides[0].getVertex(0), self.__sides[0].getVertex(1), self.__sides[1].getVertex(1) ] def getOpposingFace(self): ''' Get the face that shares this face's long side. If there is no such face, returns None. ''' opposingFace = None if self.__sides[2].getFace(0) != self: opposingFace = self.__sides[2].getFace(0) else: opposingFace = self.__sides[2].getFace(1) return opposingFace def isProcessed(self): return self.__isProcessed def setProcessed(self, val): self.__isProcessed = val def toggleProcessed(self): self.__isProcessed = self.__isProcessed == 0 def split(self): extraFaces = [] opFace = self.getOpposingFace() if opFace != None and self.__sides[2] != opFace.getHypotenuse(): extraFaces += opFace.split() opFace = self.getOpposingFace() point = self.getOrderedSides() midpt = self.__sides[2].getMidpoint() newSide = [ PlanarEdge(point[0],midpt,None,None), PlanarEdge(midpt,point[2],None,None) ] extraFaces.append( self.splitFace( point[0],point[1],point[2],midpt,newSide, 0 ) ) if opFace != None: extraFaces.append( opFace.splitFace( point[0],opFace.getCommonVertex(),point[2],midpt,newSide,1 ) ) return extraFaces def splitFace(self,pt0,pt1,pt2,midpt,newSide,newFaceIndex): diagonalSide = PlanarEdge( pt1, midpt, self, None ) if pt0 in self.__sides[0]: oldSide = [ self.__sides[0], self.__sides[1] ] self.__sides[2] = self.__sides[0] self.__sides[0] = newSide[0] self.__sides[1] = diagonalSide else: oldSide = [ self.__sides[1], self.__sides[0] ] self.__sides[2] = self.__sides[1] self.__sides[0] = newSide[0] self.__sides[1] = diagonalSide newFace = PlanarFace( diagonalSide, newSide[1], oldSide[1] ) newSide[0].setFace( newFaceIndex, self ) newSide[1].setFace( newFaceIndex, newFace ) oldFace = oldSide[1].getFace(0) if oldFace == self: oldSide[1].setFace( 0, newFace ) else: oldSide[1].setFace( 1, newFace ) diagonalSide.setFace( 1, newFace ) return newFace
PlanarDomain
The PlanarDomain class generates the initial grid, and holds all of the faces contained in the grid.
''' PlanarDomain: Creates points to use in plotting surfaces as graphs ''' class PlanarDomain: def __init__(self,xmin,xmax,ymin,ymax,nx,ny): self.__xmin = xmin self.__xmax = xmax self.__ymin = ymin self.__ymax = ymax self.__nx = nx self.__ny = ny self.__dx = (xmax - xmin) / nx self.__dy = (ymax - ymin) / ny self.__facelist = self.__populateList() def getXDivisions(self): return self.__nx def getYDivisions(self): return self.__ny def getXChange(self): return self.__dx def getYChange(self): return self.__dy def getXmin(self): return self.__xmin def getXmax(self): return self.__xmax def getYmin(self): return self.__ymin def getYmax(self): return self.__ymax ''' Generates the initial grid based on xmin,xmax , ymin,ymax nx,ny ''' def __populateList(self): table = [ [None for i in range(0,self.__ny)] for j in range (0,self.__nx) ] table2 = [ [None for i in range(0,self.__ny)] for j in range (0,self.__nx) ] faceTable = [ table, table2 ] lowerLeft = ( self.__xmin, self.__ymin ) lowerRight = ( self.__xmin + self.__dx, self.__ymin ) upperLeft = ( self.__xmin, self.__ymin + self.__dy ) upperRight = ( self.__xmin + self.__dx, self.__ymin + self.__dy ) pts = [ upperLeft, upperRight ] base = PlanarEdge( lowerLeft, lowerRight, None, None ) right = PlanarEdge( lowerRight, upperRight, None, None ) top = PlanarEdge( upperLeft, upperRight, None, None ) left = PlanarEdge( lowerLeft, upperLeft, None, None ) diagonal = PlanarEdge( lowerLeft, upperRight, None, None ) topEdges = [ top ] faceTable[0][0][0] = PlanarFace( top, left, diagonal ) faceTable[1][0][0] = PlanarFace( right, base, diagonal ) diagonal.setFace( 0, faceTable[0][0][0] ) diagonal.setFace( 1, faceTable[1][0][0] ) base.setFace( 0, faceTable[1][0][0] ) right.setFace( 0, faceTable[1][0][0] ) top.setFace( 0, faceTable[0][0][0] ) left.setFace( 0, faceTable[0][0][0] ) for i in range(1, self.__nx): lowerLeft = lowerRight lowerRight = ( self.__xmin + (i + 1)*self.__dx, self.__ymin ) upperLeft = upperRight upperRight = ( self.__xmin + (i + 1)*self.__dx, self.__ymin + self.__dy ) pts.append( upperRight ) left = right base = PlanarEdge( lowerLeft, lowerRight, None, None ) right = PlanarEdge( lowerRight, upperRight, None, None ) top = PlanarEdge( upperLeft, upperRight, None, None ) topEdges.append( top ) diagonal = PlanarEdge( lowerLeft, upperRight, None, None ) faceTable[0][i][0] = PlanarFace( top, left, diagonal ) faceTable[1][i][0] = PlanarFace( right, base, diagonal ) diagonal.setFace( 0, faceTable[0][i][0] ) diagonal.setFace( 1, faceTable[1][i][0] ) base.setFace( 0, faceTable[0][i][0] ) right.setFace( 0, faceTable[1][i][0] ) top.setFace( 0, faceTable[0][i][0] ) left.setFace( 1, faceTable[0][i][0] ) for j in range(1, self.__ny): lowerLeft = pts[0] lowerRight = pts[1] upperLeft = ( self.__xmin, self.__ymin + (j+1)*self.__dy ) upperRight = ( self.__xmin + self.__dx, self.__ymin + (j+1)*self.__dy ) pts[0] = upperLeft pts[1] = upperRight left = PlanarEdge( lowerLeft, upperLeft, None, None ) base = topEdges[0] right = PlanarEdge( lowerRight, upperRight, None, None ) top = PlanarEdge( upperLeft, upperRight, None, None ) topEdges[0] = top diagonal = PlanarEdge( lowerLeft, upperRight, None, None ) faceTable[0][0][j] = PlanarFace( top, left, diagonal ) faceTable[1][0][j] = PlanarFace( right, base, diagonal ) diagonal.setFace( 0, faceTable[0][0][j] ) diagonal.setFace( 1, faceTable[1][0][j] ) base.setFace( 1, faceTable[1][0][j] ) top.setFace( 0, faceTable[0][0][j] ) left.setFace( 0, faceTable[0][0][j] ) for i in range(1,self.__nx): lowerLeft = lowerRight lowerRight = pts[i + 1] upperLeft = upperRight upperRight = ( self.__xmin + (i+1)*self.__dx, self.__ymin + (j+1)*self.__dy ) pts[i+1] = upperRight left = right base = topEdges[i] right = PlanarEdge( lowerRight, upperRight, None, None ) top = PlanarEdge( upperLeft, upperRight, None, None ) topEdges[i] = top diagonal = PlanarEdge( lowerLeft, upperRight, None, None ) faceTable[0][i][j] = PlanarFace( top, left, diagonal ) faceTable[1][i][j] = PlanarFace( right, base, diagonal ) diagonal.setFace( 0, faceTable[0][i][j] ) diagonal.setFace( 1, faceTable[1][i][j] ) base.setFace( 1, faceTable[1][i][j] ) right.setFace( 0, faceTable[1][i][j] ) top.setFace( 0, faceTable[0][i][j] ) left.setFace(1, faceTable[0][i][j] ) result = [] for i in range(0, self.__ny): for j in range(0,self.__nx): result.append( faceTable[0][i][j] ) result.append( faceTable[1][i][j] ) return result def getFaces(self): return self.__facelist def setFaces(self,faces): self.__facelist = faces def addFaces(self,faces): self.__facelist.append( faces )
AdjContourPlot
AdjContourPlot is the class used by end users to generate implicit plots.
''' AdjContourPlot: Generates an implicit or contour plot using an adjustable mesh algorithm ''' class AdjContourPlot: #Spefifies the threshold used in determing the rate of change MAX_CURVATURE = 0.2 def __init__(self,func,contours,xmin,xmax,ymin,ymax,nx,ny,showMesh,maxRecursion): self.__func = func self.__contours = contours self.__domain = PlanarDomain( xmin, xmax, ymin, ymax, nx, ny ) self.__showMesh = showMesh self.__maxRecursion = maxRecursion self.__domain.setFaces( self.__splitFaces() ) def __dot(self,v1,v2): s = sum( [ v1[x]*v2[x] for x in range(0,len(v1)) ] ) return s def __f(self, p ): return self.__func( float(p[0]), float(p[1]) ) def __angle(self,face): corner = face.getCommonVertex() midpoint = face.getHypotenuse().getMidpoint() dx = midpoint[0] - corner[0] dy = midpoint[1] - corner[1] opposingCorner = ( midpoint[0] + dx, midpoint[1] + dy ) heightCorner = self.__f( corner ) heightMidpoint = self.__f( midpoint ) heightOpposingCorner = self.__f( opposingCorner ) vectorCorner = ( dx, dy, heightMidpoint - heightCorner ) vectorMidpoint = ( dx, dy, heightOpposingCorner - heightMidpoint ) cosine = self.__dot( vectorCorner, vectorMidpoint ) / sqrt( RealNumber(self.__dot( vectorCorner, vectorCorner ),10) * RealNumber(self.__dot( vectorMidpoint,vectorMidpoint ),10) ) if RealNumber(cosine,10).n() > 1: cosine = 1 elif RealNumber(cosine,10).n() < -1: cosine = -1 return acos( cosine ) def __splitFaces(self): facesToCheck = self.__domain.getFaces() processedFaces = [] dx = self.__domain.getXChange() dy = self.__domain.getYChange() dLength = sqrt( dx^2 + dy^2 ) / 2 for i in range( 0, self.__maxRecursion ): partiallyProcessedFaces = [] for cur in facesToCheck: cur.setProcessed( False ) for k in range( 0, len( facesToCheck ) ): face = facesToCheck[k] if face.isProcessed() == False: if float(self.__angle( face )) > float(AdjContourPlot.MAX_CURVATURE * float(dLength)): newFaces = face.split() face.setProcessed( True ) partiallyProcessedFaces.append( face ) for cur in newFaces: cur.setProcessed( True ) partiallyProcessedFaces.append( cur ) else: processedFaces.append( face ) facesToCheck = partiallyProcessedFaces processedFaces += facesToCheck return processedFaces def __getContours(self): result = [] for cur in self.__domain.getFaces(): vertex = cur.getVerticies() height = [ self.__f( x ) for x in vertex ] for curContour in self.__contours: intersectionPoints = [] for i in range(0,3): iNext = (i+1)%3 if ( height[i] - curContour ) * ( height[iNext] - curContour ) <= 0 and height[i] != height[iNext]: t = (height[iNext]-curContour)/(height[iNext]-height[i]) intersectionPoints.append( ( vertex[iNext][0] + t*(vertex[i][0] - vertex[iNext][0] ), vertex[iNext][1] + t*(vertex[i][1] - vertex[iNext][1] ) ) ) if len( intersectionPoints ) == 2: result.append( PlanarEdge( intersectionPoints[0], intersectionPoints[1], None, None ) ) return result def plot(self): g = Graphics() xmin = self.__domain.getXmin() xmax = self.__domain.getXmax() ymin = self.__domain.getYmin() ymax = self.__domain.getYmax() if self.__showMesh: color = (1,0,0) faces = self.__domain.getFaces() for cur in faces: edges = cur.getEdges() for curEdge in edges: g = curEdge.drawEdge( g, (0,0,1) ) edges = self.__getContours() color = (0,0,0) for cur in edges: g = cur.drawEdge( g, color ) return g
Usage
The constructor for AdjContourPlot takes the following parameters:
- func: A lambda or symbolic function of two variables
- contours: A list of contours to plot
- xmin: The minimum x value to include in the plot
- xmax: The maximum x value to include in the plot
- ymin: The minimum y value to include in the plot
- ymax: The maximum y value to include in the plot
- nx: The number of initial x divisions in the grid
- ny: The number of initial y divisions in the grid
- showMesh: A boolean specifies if the mesh should be displayed in the plot
- maxRecursion: The maximum number of times the algorithm will refine the grid
Example
Plot of
Sage Code:
function = lambda x,y : cos(x+y) - (x*y)^2 myPlot = AdjContourPlot( function, [0], -3, 3, -3, 3, 20, 20, False, 2 ) myPlot.plot()
Resulting Image:
Resulting Image with the Mesh:
Download
Click here to download the Sage worksheet file with the adjustable mesh algorithm.