Adjustable Mesh Plot
From www.norsemathology.org
Contents |
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 cos(x + y) - (xy)2
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.
