Adjustable Mesh Plot

From www.norsemathology.org

Jump to: navigation, search

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.

Image: Grid.png

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.

Image: Vectors.png

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.

Image: Grid_sub_divided.png

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:

Image: Sage_plot.png

Resulting Image with the Mesh:

Image: Grid_big.png

Download

Click here to download the Sage worksheet file with the adjustable mesh algorithm.

Personal tools