path: root/Build/source/utils/asymptote/base/smoothcontour3.asy

// Copyright 2015 Charles Staats III
// 
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
// 
//     http://www.apache.org/licenses/LICENSE-2.0
// 
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

// smoothcontour3
// An Asymptote module for drawing smooth implicitly defined surfaces
// author: Charles Staats III
// charles dot staats dot iii at gmail dot com

import graph_settings;  // for nmesh
import three;

/***********************************************/
/******** LINEAR ALGEBRA ROUTINES **************/
/********     LEAST-SQUARES       **************/
/***********************************************/

// Apply a matrix to a vector.
real[] apply(real[][] matrix, real[] v) {
  real[] ans = new real[matrix.length];
  for (int r = 0; r < matrix.length; ++r) {
    ans[r] = 0;
    for (int c = 0; c < v.length; ++c) {
      ans[r] += matrix[r][c] * v[c];
    }
  }
  return ans;
}

// Apply the transpose of a matrix to a vector,
// without actually forming the transpose.
real[] applytranspose(real[][] matrix, real[] v) {
  real[] ans = new real[matrix[0].length];
  for (int r = 0; r < ans.length; ++r) ans[r] = 0;
  for (int c = 0; c < matrix.length; ++c) {
    for (int r = 0; r < ans.length; ++r) {
      ans[r] += matrix[c][r] * v[c];
    }
  }
  return ans;
}

// For a matrix A, returns the matrix product
// (A transposed) * A.
// The transpose of A is never actually formed.
real[][] AtA(real[][] matrix) {
  real[][] toreturn = new real[matrix[0].length][matrix[0].length];
  for (int i = 0; i < toreturn.length; ++i) {
    for (int j = 0; j < toreturn.length; ++j) {
      toreturn [i][j] = 0;
    }
  }
  for (int k = 0; k < matrix.length; ++k) {
    for (int i = 0; i < toreturn.length; ++i) {
      for (int j = 0; j < toreturn.length; ++j) {
	toreturn[i][j] += matrix[k][i] * matrix[k][j];
      }
    }
  }
  return toreturn;
}

// Assuming A is a matrix with independent columns, returns
// the unique vector y minimizing |Ay - b|^2 (the L2 norm).
// If the columns of A are not linearly independent,
// throws an error (if warn == true) or returns an empty array
// (if warn == false).
real[] leastsquares(real[][] A, real[] b, bool warn = true) {
  real[] solution = solve(AtA(A), applytranspose(A, b), warn=false);
  if (solution.length == 0 && warn)
    abort("Cannot compute least-squares approximation for " +
	  "a matrix with linearly dependent columns.");
  return solution;
}

/***********************************************/
/******** CREATING BEZIER PATCHES **************/
/******** WITH SPECIFIED NORMALS  **************/
/***********************************************/

// The weight given to minimizing the sum of squares of
// the mixed partials at the corners of the bezier patch.
// If this weight is zero, the result is undefined in
// places and can be rather wild even where it is
// defined.
// The struct is used to as a namespace.
struct pathwithnormals_settings {
  static real wildnessweight = 1e-3;
}
private from pathwithnormals_settings unravel wildnessweight;

// The Bernstein basis polynomials of degree 3:
real B03(real t) { return (1-t)^3; }
real B13(real t) { return 3*t*(1-t)^2; }
real B23(real t) { return 3*t^2*(1-t); }
real B33(real t) { return t^3; }

private typedef real function(real);
function[] bernstein = new function[] {B03, B13, B23, B33};

// This function attempts to produce a Bezier patch
// with the specified boundary path and normal directions.
// For instance, the patch should be normal to
// u0normals[0] at (0, 0.25),
// normal to u0normals[1] at (0, 0.5), and
// normal to u0normals[2] at (0, 0.75).
// The actual normal (as computed by the patch.normal() function)
// may be parallel to the specified normal, antiparallel, or
// even zero.
//
// A small amount of deviation is allowed in order to stabilize
// the algorithm (by keeping the mixed partials at the corners from
// growing too large).
//
// Note that the specified normals are projected to be orthogonal to
// the specified boundary path. However, the entries in the array
// remain intact.
patch patchwithnormals(path3 external, triple[] u0normals, triple[] u1normals,
		       triple[] v0normals, triple[] v1normals)
{
  assert(cyclic(external));
  assert(length(external) == 4);
  assert(u0normals.length == 3);
  assert(u1normals.length == 3);
  assert(v0normals.length == 3);
  assert(v1normals.length == 3);
  
  triple[][] controlpoints = new triple[4][4];
  controlpoints[0][0] = point(external,0);
  controlpoints[1][0] = postcontrol(external,0);
  controlpoints[2][0] = precontrol(external,1);
  controlpoints[3][0] = point(external,1);
  controlpoints[3][1] = postcontrol(external,1);
  controlpoints[3][2] = precontrol(external,2);
  controlpoints[3][3] = point(external,2);
  controlpoints[2][3] = postcontrol(external,2);
  controlpoints[1][3] = precontrol(external,3);
  controlpoints[0][3] = point(external,3);
  controlpoints[0][2] = postcontrol(external,3);
  controlpoints[0][1] = precontrol(external, 4);

  real[][] matrix = new real[24][12];
  for (int i = 0; i < matrix.length; ++i)
    for (int j = 0; j < matrix[i].length; ++j)
      matrix[i][j] = 0;
  real[] rightvector = new real[24];
  for (int i = 0; i < rightvector.length; ++i)
    rightvector[i] = 0;

  void addtocoeff(int i, int j, int count, triple coeffs) {
    if (1 <= i && i <= 2 && 1 <= j && j <= 2) {
      int position = 3 * (2 * (i-1) + (j-1));
      matrix[count][position] += coeffs.x;
      matrix[count][position+1] += coeffs.y;
      matrix[count][position+2] += coeffs.z;
    } else {
      rightvector[count] -= dot(controlpoints[i][j], coeffs);
    }
  }

  void addtocoeff(int i, int j, int count, real coeff) {
    if (1 <= i && i <= 2 && 1 <= j && j <= 2) {
      int position = 3 * (2 * (i-1) + (j-1));
      matrix[count][position] += coeff;
      matrix[count+1][position+1] += coeff;
      matrix[count+2][position+2] += coeff;
    } else {
      rightvector[count] -= controlpoints[i][j].x * coeff;
      rightvector[count+1] -= controlpoints[i][j].y * coeff;
      rightvector[count+2] -= controlpoints[i][j].z * coeff;
    }
  }

  int count = 0;

  void apply_u0(int j, real a, triple n) {
    real factor = 3 * bernstein[j](a);
    addtocoeff(0,j,count,-factor*n);
    addtocoeff(1,j,count,factor*n);
  }
  void apply_u0(real a, triple n) {
    triple tangent = dir(external, 4-a);
    n -= dot(n,tangent)*tangent;
    n = unit(n);
    for (int j = 0; j < 4; ++j) {
      apply_u0(j,a,n);
    }
    ++count;
  }
  apply_u0(0.25, u0normals[0]);
  apply_u0(0.5, u0normals[1]);
  apply_u0(0.75, u0normals[2]);

  void apply_u1(int j, real a, triple n) {
    real factor = 3 * bernstein[j](a);
    addtocoeff(3,j,count,factor*n);
    addtocoeff(2,j,count,-factor*n);
  }
  void apply_u1(real a, triple n) {
    triple tangent = dir(external, 1+a);
    n -= dot(n,tangent)*tangent;
    n = unit(n);
    for (int j = 0; j < 4; ++j)
      apply_u1(j,a,n);
    ++count;
  }
  apply_u1(0.25, u1normals[0]);
  apply_u1(0.5, u1normals[1]);
  apply_u1(0.75, u1normals[2]);

  void apply_v0(int i, real a, triple n) {
    real factor = 3 * bernstein[i](a);
    addtocoeff(i,0,count,-factor*n);
    addtocoeff(i,1,count,factor*n);
  }
  void apply_v0(real a, triple n) {
    triple tangent = dir(external, a);
    n -= dot(n,tangent) * tangent;
    n = unit(n);
    for (int i = 0; i < 4; ++i)
      apply_v0(i,a,n);
    ++count;
  }
  apply_v0(0.25, v0normals[0]);
  apply_v0(0.5, v0normals[1]);
  apply_v0(0.75, v0normals[2]);

  void apply_v1(int i, real a, triple n) {
    real factor = 3 * bernstein[i](a);
    addtocoeff(i,3,count,factor*n);
    addtocoeff(i,2,count,-factor*n);
  }
  void apply_v1(real a, triple n) {
    triple tangent = dir(external, 3-a);
    n -= dot(n,tangent)*tangent;
    n = unit(n);
    for (int i = 0; i < 4; ++i)
      apply_v1(i,a,n);
    ++count;
  }
  apply_v1(0.25, v1normals[0]);
  apply_v1(0.5, v1normals[1]);
  apply_v1(0.75, v1normals[2]);

  addtocoeff(0,0,count,9*wildnessweight);
  addtocoeff(1,1,count,9*wildnessweight);
  addtocoeff(0,1,count,-9*wildnessweight);
  addtocoeff(1,0,count,-9*wildnessweight);
  count+=3;
  addtocoeff(3,3,count,9*wildnessweight);
  addtocoeff(2,2,count,9*wildnessweight);
  addtocoeff(3,2,count,-9*wildnessweight);
  addtocoeff(2,3,count,-9*wildnessweight);
  count+=3;
  addtocoeff(0,3,count,9*wildnessweight);
  addtocoeff(1,2,count,9*wildnessweight);
  addtocoeff(1,3,count,-9*wildnessweight);
  addtocoeff(0,2,count,-9*wildnessweight);
  count += 3;
  addtocoeff(3,0,count,9*wildnessweight);
  addtocoeff(2,1,count,9*wildnessweight);
  addtocoeff(3,1,count,-9*wildnessweight);
  addtocoeff(2,0,count,-9*wildnessweight);
  count += 3;

  real[] solution = leastsquares(matrix, rightvector, warn=false);
  if (solution.length == 0) { // if the matrix was singular
    write("Warning: unable to solve matrix for specifying normals "
	  + "on bezier patch. Using standard method.");
    return patch(external);
  }
  
  for (int i = 1; i <= 2; ++i) {
    for (int j = 1; j <= 2; ++j) {
      int position = 3 * (2 * (i-1) + (j-1));
      controlpoints[i][j] = (solution[position],
			     solution[position+1],
			     solution[position+2]);
    }
  }

  return patch(controlpoints);
}

// A wrapper for the previous function when the normal direction
// is given as a function of direction. The wrapper can also
// accommodate cyclic boundary paths of between one and four
// segments, although the results are best by far when there
// are four segments.
patch patchwithnormals(path3 external, triple normalat(triple)) {
  assert(cyclic(external));
  assert(1 <= length(external) && length(external) <= 4);
  while (length(external) < 4) external = external -- cycle;
  triple[] u0normals = new triple[3];
  triple[] u1normals = new triple[3];
  triple[] v0normals = new triple[3];
  triple[] v1normals = new triple[3];
  for (int i = 1; i <= 3; ++i) {
    v0normals[i-1] = unit(normalat(point(external, i/4)));
    u1normals[i-1] = unit(normalat(point(external, 1 + i/4)));
    v1normals[i-1] = unit(normalat(point(external, 3 - i/4)));
    u0normals[i-1] = unit(normalat(point(external, 4 - i/4)));
  }
  return patchwithnormals(external, u0normals, u1normals, v0normals, v1normals);
}

/***********************************************/
/*********** ROOT-FINDER UTILITY ***************/
/***********************************************/

// Namespace
struct rootfinder_settings {
  static real roottolerance = 1e-4;
}

// Find a root for the specified continuous (but not
// necessarily differentiable) function. Whatever
// value t is returned, it is guaranteed that either
// t is within tolerance of a sign change, or 
// abs(f(t)) <= 0.1 tolerance.
// An error is thrown if fa and fb are both positive
// or both negative.
//
// In the current implementation, binary search is interleaved
// with a modified version of linear interpolation.
real findroot(real f(real), real a, real b,
	       real tolerance = rootfinder_settings.roottolerance,
	       real fa = f(a), real fb = f(b))
{
  if (fa == 0) return a;
  if (fb == 0) return b;
  real g(real);
  if (fa < 0) {
    assert(fb > 0);
    g = f;
  } else {
    assert(fb < 0);
    fa = -fa;
    fb = -fb;
    g = new real(real t) { return -f(t); };
  }

  real t = a;
  real ft = fa;

  while (b - a > tolerance && abs(ft) > 0.1*tolerance) {
    t = a + (b - a) / 2;
    ft = g(t);
    if (ft == 0) return t;
    else if (ft > 0) {
      b = t;
      fb = ft;
    } else if (ft < 0) {
      a = t;
      fa = ft;
    }
    
    // linear interpolation
    t = a - (b - a) / (fb - fa) * fa;

    // If the interpolated value is close to one edge of
    // the interval, move it farther away from the edge in
    // an effort to catch the root in the middle.
    if (t - a < (b-a)/8) t = a + 2*(t-a);
    else if (b - t < (b-a)/8) t = b - 2*(b-t);

    assert(t >= a && t <= b);
    
    ft = g(t);
    if (ft == 0) return t;
    else if (ft > 0) {
      b = t;
      fb = ft;
    } else if (ft < 0) {
      a = t;
      fa = ft;
    }

  }
  return a - (b - a) / (fb - fa) * fa;
}

/***********************************************/
/********* DUAL CUBE GRAPH UTILITY *************/
/***********************************************/

// Suppose a plane intersects a (hollow) cube, and
// does not intersect any vertices. Then its intersection
// with cube forms a cycle. The goal of the code below
// is to reconstruct the order of the cycle
// given only an unordered list of which edges the plane
// intersects.
//
// Basically, the question is this: If we know the points
// in which a more-or-less planar surface intersects the
// edges of cube, how do we connect those points?
//
// When I wrote the code, I was thinking in terms of the
// dual graph of a cube, in which "vertices" are really
// faces of the cube and "edges" connect those "vertices."

// An enum for the different "vertices" (i.e. faces)
// available. NULL_VERTEX is primarily intended as a
// return value to indicate the absence of a desired
// vertex.
private int NULL_VERTEX = -1;
private int XHIGH = 0;
private int XLOW = 1;
private int YHIGH = 2;
private int YLOW = 3;
private int ZHIGH = 4;
private int ZLOW = 5;

// An unordered set of nonnegative integers.
// Since the intent is to use
// only the six values from the enum above, no effort
// was made to use scalable algorithms.
struct intset {
  private bool[] ints = new bool[0];
  private int size = 0;

  bool contains(int item) {
    assert(item >= 0);
    if (item >= ints.length) return false;
    return ints[item];
  }

  // Returns true if the item was added (i.e., was
  // not already present).
  bool add(int item) {
    assert(item >= 0);
    while (item >= ints.length) ints.push(false);
    if (ints[item]) return false;
    ints[item] = true;
    ++size;
    return true;
  }

  int[] elements() {
    int[] toreturn;
    for (int i = 0; i < ints.length; ++i) {
      if (ints[i]) toreturn.push(i);
    }
    return toreturn;
  }

  int size() { return size; }
}

// A map from integers to sets of integers. Again, no
// attempt is made to use scalable data structures.
struct int_to_intset {
  int[] keys = new int[0];
  intset[] values = new intset[0];

  void add(int key, int value) {
    for (int i = 0; i < keys.length; ++i) {
      if (keys[i] == key) {
	values[i].add(value);
	return;
      }
    }
    keys.push(key);
    intset newset;
    values.push(newset);
    newset.add(value);
  }

  private int indexOf(int key) {
    for (int i = 0; i < keys.length; ++i) {
      if (keys[i] == key) return i;
    }
    return -1;
  }

  int[] get(int key) {
    int i = indexOf(key);
    if (i < 0) return new int[0];
    else return values[i].elements();
  }

  int numvalues(int key) {
    int i = indexOf(key);
    if (i < 0) return 0;
    else return values[i].size();
  }

  int numkeys() {
    return keys.length;
  }
}

// A struct intended to represent an undirected edge between
// two "vertices."
struct edge {
  int start;
  int end;
  void operator init(int a, int b) {
    start = a;
    end = b;
  }
  bool bordersvertex(int v) { return start == v || end == v; }
}

string operator cast(edge e) {
  int a, b;
  if (e.start <= e.end) {a = e.start; b = e.end;}
  else {a = e.end; b = e.start; }
  return (string)a + " <-> " + (string)b;
}

bool operator == (edge a, edge b) {
  if (a.start == b.start && a.end == b.end) return true;
  if (a.start == b.end && a.end == b.start) return true;
  return false;
}

string operator cast(edge[] edges) {
  string toreturn = "{ ";
  for (int i = 0; i < edges.length; ++i) {
    toreturn += edges[i];
    if (i < edges.length-1) toreturn += ", ";
  }
  return toreturn + " }";
}

// Finally, the function that strings together a list of edges
// into a cycle. It makes assumptions that hold true if the
// list of edges did in fact come from a plane intersection
// containing no vertices of the cube. For instance, such a
// plane can contain at most two noncollinear points of any
// one face; consequently, no face can border more than two of
// the selected edges.
//
// If the underlying assumptions prove to be false, the function
// returns null.
int[] makecircle(edge[] edges) {
  if (edges.length == 0) return new int[0];
  int_to_intset graph;
  for (edge e : edges) {
    graph.add(e.start, e.end);
    graph.add(e.end, e.start);
  }
  int currentvertex = edges[0].start;
  int startvertex = currentvertex;
  int lastvertex = NULL_VERTEX;
  int[] toreturn = new int[0];
  do {
    toreturn.push(currentvertex);
    int[] adjacentvertices = graph.get(currentvertex);
    if (adjacentvertices.length != 2) return null;
    for (int v : adjacentvertices) {
      if (v != lastvertex) {
	lastvertex = currentvertex;
	currentvertex = v;
	break;
      }
    }
  } while (currentvertex != startvertex);
  if (toreturn.length != graph.numkeys()) return null;
  toreturn.cyclic = true;
  return toreturn;
}

/***********************************************/
/********** PATHS BETWEEN POINTS ***************/
/***********************************************/
// Construct paths between two points with additional
// constraints; for instance, the path must be orthogonal
// to a certain vector at each of the endpoints, must
// lie within a specified plane or a specified face
// of a rectangular solid,....

// A vector (typically a normal vector) at a specified position.
struct positionedvector {
  triple position;
  triple direction;
  void operator init(triple position, triple direction) {
    this.position = position;
    this.direction = direction;
  }
}

string operator cast(positionedvector vv) {
  return "position: " + (string)(vv.position) + " vector: " + (string)vv.direction;
}

// The angle, in degrees, between two vectors.
real angledegrees(triple a, triple b) {
  real lengthprod = abs(a) * abs(b);
  if (lengthprod == 0) return 0;
  return aCos(dot(a,b) / lengthprod);
}

// A path (single curved segment) between two points. At each point
// is specified a vector orthogonal to the path.
path3 pathbetween(positionedvector v1, positionedvector v2) {
  triple n1 = unit(v1.direction);
  triple n2 = unit(v2.direction);

  triple p1 = v1.position;
  triple p2 = v2.position;
  triple delta = p2-p1;

  triple dir1 = delta - dot(delta, n1)*n1;
  triple dir2 = delta - dot(delta, n2)*n2;
  return p1 {dir1} .. {dir2} p2;
}

// Assuming v1 and v2 are linearly independent, returns an array {a, b}
// such that a v1 + b v2 is the orthogonal projection of toproject onto
// the span of v1 and v2. If v1 and v2 are dependent, returns an empty array
// (if warn==false) or throws an error (if warn==true).
real[] projecttospan_findcoeffs(triple toproject, triple v1, triple v2,
				   bool warn=false) {
  real[][] matrix = {{v1.x, v2.x},
		     {v1.y, v2.y},
		     {v1.z, v2.z}};
  real[] desiredanswer = {toproject.x, toproject.y, toproject.z};
  return leastsquares(matrix, desiredanswer, warn=warn);
}

// Project the triple toproject into the span of a and b, but restrict
// to the quarter-plane of linear combinations a v1 + b v2 such that
// a >= mincoeff and b >= mincoeff. If v1 and v2 are linearly dependent,
// return a random (positive) linear combination.
triple projecttospan(triple toproject, triple v1, triple v2,
		       real mincoeff = 0.05) {
  real[] coeffs = projecttospan_findcoeffs(toproject, v1, v2, warn=false);
  real a, b;
  if (coeffs.length == 0) {
    a = mincoeff + unitrand();
    b = mincoeff + unitrand();
  } else {
    a = max(coeffs[0], mincoeff);
    b = max(coeffs[1], mincoeff);
  }
  return a*v1 + b*v2;
}

// A path between two specified vertices of a cyclic path. The
// path tangent at each endpoint is guaranteed to lie within the
// quarter-plane spanned by positive linear combinations of the
// tangents of the two outgoing paths at that endpoint.
path3 pathbetween(path3 edgecycle, int vertex1, int vertex2) {
  triple point1 = point(edgecycle, vertex1);
  triple point2 = point(edgecycle, vertex2);
  
  triple v1 = -dir(edgecycle, vertex1, sign=-1);
  triple v2 =  dir(edgecycle, vertex1, sign= 1);
  triple direction1 = projecttospan(unit(point2-point1), v1, v2);

  v1 = -dir(edgecycle, vertex2, sign=-1);
  v2 =  dir(edgecycle, vertex2, sign= 1);
  triple direction2 = projecttospan(unit(point1-point2), v1, v2);

  return point1 {direction1} .. {-direction2} point2;
}

// This function applies a heuristic to choose two "opposite"
// vertices (separated by three segments) of edgecycle, which
// is required to be a cyclic path consisting of 5 or 6 segments.
// The two chosen vertices are pushed to savevertices.
//
// The function returns a path between the two chosen vertices. The
// path tangent at each endpoint is guaranteed to lie within the
// quarter-plane spanned by positive linear combinations of the
// tangents of the two outgoing paths at that endpoint.
path3 bisector(path3 edgecycle, int[] savevertices) {
  real mincoeff = 0.05;
  assert(cyclic(edgecycle));
  int n = length(edgecycle);
  assert(n >= 5 && n <= 6);
  triple[] forwarddirections = sequence(new triple(int i) {
      return dir(edgecycle, i, sign=1);
    }, n);
  forwarddirections.cyclic = true;
  triple[] backwarddirections = sequence(new triple(int i) {
      return -dir(edgecycle, i, sign=-1);
    }, n);
  backwarddirections.cyclic = true;
  real[] angles = sequence(new real(int i) {
      return angledegrees(forwarddirections[i], backwarddirections[i]);
    }, n);
  angles.cyclic = true;
  int lastindex = (n == 5 ? 4 : 2);
  real maxgoodness = 0;
  int chosenindex = -1;
  triple directionout, directionin;
  for (int i = 0; i <= lastindex; ++i) {
    int opposite = i + 3;
    triple vec = unit(point(edgecycle, opposite) - point(edgecycle, i));
    real[] coeffsbegin = projecttospan_findcoeffs(vec, forwarddirections[i],
						    backwarddirections[i]);
    if (coeffsbegin.length == 0) continue;
    coeffsbegin[0] = max(coeffsbegin[0], mincoeff);
    coeffsbegin[1] = max(coeffsbegin[1], mincoeff);

    real[] coeffsend = projecttospan_findcoeffs(-vec, forwarddirections[opposite],
						    backwarddirections[opposite]);
    if (coeffsend.length == 0) continue;
    coeffsend[0] = max(coeffsend[0], mincoeff);
    coeffsend[1] = max(coeffsend[1], mincoeff);

    real goodness = angles[i] * angles[opposite] * coeffsbegin[0] * coeffsend[0]
      * coeffsbegin[1] * coeffsend[1];
    if (goodness > maxgoodness) {
      maxgoodness = goodness;
      directionout = coeffsbegin[0] * forwarddirections[i] +
	coeffsbegin[1] * backwarddirections[i];
      directionin = -(coeffsend[0] * forwarddirections[opposite] +
		       coeffsend[1] * backwarddirections[opposite]);
      chosenindex = i;
    }
  }
  if (chosenindex == -1) {
    savevertices.push(0);
    savevertices.push(3);
    return pathbetween(edgecycle, 0, 3);
  } else {
    savevertices.push(chosenindex);
    savevertices.push(chosenindex+3);
    return point(edgecycle, chosenindex) {directionout} ..
      {directionin} point(edgecycle, chosenindex + 3);
  }
}

// A path between two specified points (with specified normals) that lies
// within a specified face of a rectangular solid.
path3 pathinface(positionedvector v1, positionedvector v2,
		    triple facenorm, triple edge1normout, triple edge2normout)
{
  triple dir1 = cross(v1.direction, facenorm);
  real dotprod = dot(dir1, edge1normout);
  if (dotprod > 0) dir1 = -dir1;
  // Believe it or not, this "tiebreaker" is actually relevant at times,
  // for instance, when graphing the cone x^2 + y^2 = z^2 over the region
  // -1 <= x,y,z <= 1.
  else if (dotprod == 0 && dot(dir1, v2.position - v1.position) < 0) dir1 = -dir1;

  triple dir2 = cross(v2.direction, facenorm);
  dotprod = dot(dir2, edge2normout);
  if (dotprod < 0) dir2 = -dir2;
  else if (dotprod == 0 && dot(dir2, v2.position - v1.position) < 0) dir2 = -dir2;

  return v1.position {dir1} .. {dir2} v2.position;
}

triple normalout(int face) {
  if (face == XHIGH) return X;
  else if (face == YHIGH) return Y;
  else if (face == ZHIGH) return Z;
  else if (face == XLOW) return -X;
  else if (face == YLOW) return -Y;
  else if (face == ZLOW) return -Z;
  else return O;
}

// A path between two specified points (with specified normals) that lies
// within a specified face of a rectangular solid.
path3 pathinface(positionedvector v1, positionedvector v2,
		    int face, int edge1face, int edge2face) {
  return pathinface(v1, v2, normalout(face), normalout(edge1face),
		      normalout(edge2face));
}

/***********************************************/
/******** DRAWING IMPLICIT SURFACES ************/
/***********************************************/

// Quadrilateralization:
// Produce a surface (array of *nondegenerate* Bezier patches) with a
// specified three-segment boundary. The surface should approximate the
// zero locus of the specified f with its specified gradient.
//
// If it is not possible to produce the desired result without leaving the
// specified rectangular region, returns a length-zero array.
//
// Dividing a triangle into smaller quadrilaterals this way is opposite
// the usual trend in mathematics. However, the pathwithnormals algorithm
// does a poor job of choosing a good surface when the boundary path does
// not consist of four positive-length segments.
patch[] triangletoquads(path3 external, real f(triple), triple grad(triple),
			  triple a, triple b) {
  static real epsilon = 1e-3;
  assert(length(external) == 3);
  assert(cyclic(external));

  triple c0 = point(external, 0);
  triple c1 = point(external, 1);
  triple c2 = point(external, 2);

  triple center = (c0 + c1 + c2) / 3;
  triple n = unit(cross(c1-c0, c2-c0));

  real g(real t) { return f(center + t*n); }

  real tmin = -realMax, tmax = realMax;
  void absorb(real t) {
    if (t < 0) tmin = max(t,tmin);
    else tmax = min(t,tmax);
  }
  if (n.x != 0) {
    absorb((a.x - center.x) / n.x);
    absorb((b.x - center.x) / n.x);
  }
  if (n.y != 0) {
    absorb((a.y - center.y) / n.y);
    absorb((b.y - center.y) / n.y);
  }
  if (n.z != 0) {
    absorb((a.z - center.z) / n.z);
    absorb((b.z - center.z) / n.z);
  }

  real fa = g(tmin);
  real fb = g(tmax);
  if ((fa > 0 && fb > 0) || (fa < 0 && fb < 0)) {
    return new patch[0];
  } else {
    real t = findroot(g, tmin, tmax, fa=fa, fb=fb);
    center += t * n;
  }

  n = unit(grad(center));

  triple m0 = point(external, 0.5);
  positionedvector m0 = positionedvector(m0, unit(grad(m0)));
  triple m1 = point(external, 1.5);
  positionedvector m1 = positionedvector(m1, unit(grad(m1)));
  triple m2 = point(external, 2.5);
  positionedvector m2 = positionedvector(m2, unit(grad(m2)));
  positionedvector c = positionedvector(center, unit(grad(center)));

  path3 pathto_m0 = pathbetween(c, m0);
  path3 pathto_m1 = pathbetween(c, m1);
  path3 pathto_m2 = pathbetween(c, m2);

  path3 quad0 = subpath(external, 0, 0.5)
    & reverse(pathto_m0)
    & pathto_m2
    & subpath(external, -0.5, 0)
    & cycle;
  path3 quad1 = subpath(external, 1, 1.5)
    & reverse(pathto_m1)
    & pathto_m0
    & subpath(external, 0.5, 1)
    & cycle;
  path3 quad2 = subpath(external, 2, 2.5)
    & reverse(pathto_m2)
    & pathto_m1
    & subpath(external, 1.5, 2)
    & cycle;

  return new patch[] {patchwithnormals(quad0, grad),
      patchwithnormals(quad1, grad),
      patchwithnormals(quad2, grad)};   
}

// Returns true if the point is "nonsingular" (in the sense that the magnitude
// of the gradient is not too small) AND very close to the zero locus of f
// (assuming f is locally linear).
bool check_fpt_zero(triple testpoint, real f(triple), triple grad(triple)) {
  real testval = f(testpoint);
  real slope = abs(grad(testpoint));
  static real tolerance = 2*rootfinder_settings.roottolerance;
  return !(slope > tolerance && abs(testval) / slope > tolerance);
}

// Returns true if pt lies within the rectangular solid with
// opposite corners at a and b.
bool checkptincube(triple pt, triple a, triple b) {
  real xmin = a.x;
  real xmax = b.x;
  real ymin = a.y;
  real ymax = b.y;
  real zmin = a.z;
  real zmax = b.z;
  if (xmin > xmax) { real t = xmax; xmax=xmin; xmin=t; }
  if (ymin > ymax) { real t = ymax; ymax=ymin; ymin=t; }
  if (zmin > zmax) { real t = zmax; zmax=zmin; zmin=t; }

  return ((xmin <= pt.x) && (pt.x <= xmax) &&
	  (ymin <= pt.y) && (pt.y <= ymax) &&
	  (zmin <= pt.z) && (pt.z <= zmax));

}

// A convenience function for combining the previous two tests.
bool checkpt(triple testpt, real f(triple), triple grad(triple),
	      triple a, triple b) {
  return checkptincube(testpt, a, b) &&
    check_fpt_zero(testpt, f, grad);
}

// Attempts to fill in the boundary cycle with a collection of
// patches to approximate smoothly the zero locus of f. If unable to
// do so while satisfying certain checks, returns null.
// This is distinct from returning an empty
// array, which merely indicates that the boundary cycle is too small
// to be worth filling in.
patch[] quadpatches(path3 edgecycle, positionedvector[] corners,
		     real f(triple), triple grad(triple),
		     triple a, triple b) {
  assert(corners.cyclic);

  // The tolerance for considering two points "essentially identical."
  static real tolerance = 2.5 * rootfinder_settings.roottolerance;
  
  // If there are two neighboring vertices that are essentially identical,
  // unify them into one.
  for (int i = 0; i < corners.length; ++i) {
    if (abs(corners[i].position - corners[i+1].position) < tolerance) {
      if (corners.length == 2) return new patch[0];
      corners.delete(i);
      edgecycle = subpath(edgecycle, 0, i)
	& subpath(edgecycle, i+1, length(edgecycle))
	& cycle;
      --i;
      assert(length(edgecycle) == corners.length);
    }
  }

  static real areatolerance = tolerance^2;

  assert(corners.length >= 2);
  if (corners.length == 2) {
    // If the area is too small, just ignore it; otherwise, subdivide.
    real area0 = abs(cross(-dir(edgecycle, 0, sign=-1, normalize=false),
			   dir(edgecycle, 0, sign=1, normalize=false)));
    real area1 = abs(cross(-dir(edgecycle, 1, sign=-1, normalize=false),
			   dir(edgecycle, 1, sign=1, normalize=false)));
    if (area0 < areatolerance && area1 < areatolerance) return new patch[0];
    else return null;
  }
  if (length(edgecycle) > 6) abort("too many edges: not possible.");

  for (int i = 0; i < length(edgecycle); ++i) {
    if (angledegrees(dir(edgecycle,i,sign=1),
		     dir(edgecycle,i+1,sign=-1)) > 80) {
      return null;
    }
  }
  
  if (length(edgecycle) == 3) {
    patch[] toreturn = triangletoquads(edgecycle, f, grad, a, b);
    if (toreturn.length == 0) return null;
    else return toreturn;
  }
  if (length(edgecycle) == 4) {
    return new patch[] {patchwithnormals(edgecycle, grad)};
  }

  int[] bisectorindices;
  path3 middleguide = bisector(edgecycle, bisectorindices);
  
  triple testpoint = point(middleguide, 0.5);
  if (!checkpt(testpoint, f, grad, a, b)) {
    return null;
  }
  
  patch[] toreturn = null;
  path3 firstpatch = subpath(edgecycle, bisectorindices[0], bisectorindices[1])
    & reverse(middleguide) & cycle;
  if (length(edgecycle) == 5) {
    path3 secondpatch = middleguide
	& subpath(edgecycle, bisectorindices[1], 5+bisectorindices[0]) & cycle;
    toreturn = triangletoquads(secondpatch, f, grad, a, b);
    if (toreturn.length == 0) return null;
    toreturn.push(patchwithnormals(firstpatch, grad));
  } else {
    // now length(edgecycle) == 6
    path3 secondpatch = middleguide
      & subpath(edgecycle, bisectorindices[1], 6+bisectorindices[0])
      & cycle;
    toreturn = new patch[] {patchwithnormals(firstpatch, grad),
	patchwithnormals(secondpatch, grad)};
  }
  return toreturn;
}

// Numerical gradient of a function
typedef triple vectorfunction(triple);
vectorfunction nGrad(real f(triple)) {
  static real epsilon = 1e-3;
  return new triple(triple v) {
    return ( (f(v + epsilon*X) - f(v - epsilon*X)) / (2 epsilon),
			 (f(v + epsilon*Y) - f(v - epsilon*Y)) / (2 epsilon),
			 (f(v + epsilon*Z) - f(v - epsilon*Z)) / (2 epsilon) );
  };
}

// A point together with a value at that location.
struct evaluatedpoint {
  triple pt;
  real value;
  void operator init(triple pt, real value) {
    this.pt = pt;
    this.value = value;
  }
}

triple operator cast(evaluatedpoint p) { return p.pt; }

// Compute the values of a function at every vertex of an nx by ny by nz
// array of rectangular solids.
evaluatedpoint[][][] make3dgrid(triple a, triple b, int nx, int ny, int nz,
				   real f(triple), bool allowzero = false)
{
  evaluatedpoint[][][] toreturn = new evaluatedpoint[nx+1][ny+1][nz+1];
  for (int i = 0; i <= nx; ++i) {
    for (int j = 0; j <= ny; ++j) {
      for (int k = 0; k <= nz; ++k) {
	triple pt = (interp(a.x, b.x, i/nx),
		     interp(a.y, b.y, j/ny),
		     interp(a.z, b.z, k/nz));
	real value = f(pt);
	if (value == 0 && !allowzero) value = 1e-5;
	toreturn[i][j][k] = evaluatedpoint(pt, value);
      }
    }
  }
  return toreturn;
}

// The following utilities make, for instance, slice(A, i, j, k, l)
// equivalent to what A[i:j][k:l] ought to mean for two- and three-
// -dimensional arrays of evaluatedpoints and of positionedvectors.
typedef evaluatedpoint T;
T[][] slice(T[][] a, int start1, int end1, int start2, int end2) {
  T[][] toreturn = new T[end1-start1][];
  for (int i = start1; i < end1; ++i) {
    toreturn[i-start1] = a[i][start2:end2];
  }
  return toreturn;
}
T[][][] slice(T[][][] a, int start1, int end1,
	      int start2, int end2,
	      int start3, int end3) {
  T[][][] toreturn = new T[end1-start1][][];
  for (int i = start1; i < end1; ++i) {
    toreturn[i-start1] = slice(a[i], start2, end2, start3, end3);
  }
  return toreturn;
}
typedef positionedvector T;
T[][] slice(T[][] a, int start1, int end1, int start2, int end2) {
  T[][] toreturn = new T[end1-start1][];
  for (int i = start1; i < end1; ++i) {
    toreturn[i-start1] = a[i][start2:end2];
  }
  return toreturn;
}
T[][][] slice(T[][][] a, int start1, int end1,
	      int start2, int end2,
	      int start3, int end3) {
  T[][][] toreturn = new T[end1-start1][][];
  for (int i = start1; i < end1; ++i) {
    toreturn[i-start1] = slice(a[i], start2, end2, start3, end3);
  }
  return toreturn;
}

// An object of class gridwithzeros stores the values of a function at each vertex
// of a three-dimensional grid, together with zeros of the function along edges
// of the grid and the gradient of the function at each such zero.
struct gridwithzeros {
  int nx, ny, nz;
  evaluatedpoint[][][] corners;
  positionedvector[][][] xdirzeros;
  positionedvector[][][] ydirzeros;
  positionedvector[][][] zdirzeros;
  triple grad(triple);
  real f(triple);
  int maxdepth;

  // Populate the edges with zeros that have a sign change and are not already
  // populated.
  void fillzeros() {
    for (int j = 0; j < ny+1; ++j) {
      for (int k = 0; k < nz+1; ++k) {
	real y = corners[0][j][k].pt.y;
	real z = corners[0][j][k].pt.z;
	real f_along_x(real t) { return f((t, y, z)); }
	for (int i = 0; i < nx; ++i) {
	  if (xdirzeros[i][j][k] != null) continue;
	  evaluatedpoint start = corners[i][j][k];
	  evaluatedpoint end = corners[i+1][j][k];
	  if ((start.value > 0 && end.value > 0) || (start.value < 0 && end.value < 0))
	    xdirzeros[i][j][k] = null;
	  else {
	    triple root = (0,y,z);
	    root += X * findroot(f_along_x, start.pt.x, end.pt.x,
				    fa=start.value, fb=end.value);
	    triple normal = grad(root);
	    xdirzeros[i][j][k] = positionedvector(root, normal);
	  }
	}
      }
    }

    for (int i = 0; i < nx+1; ++i) {
      for (int k = 0; k < nz+1; ++k) {
	real x = corners[i][0][k].pt.x;
	real z = corners[i][0][k].pt.z;
	real f_along_y(real t) { return f((x, t, z)); }
	for (int j = 0; j < ny; ++j) {
	  if (ydirzeros[i][j][k] != null) continue;
	  evaluatedpoint start = corners[i][j][k];
	  evaluatedpoint end = corners[i][j+1][k];
	  if ((start.value > 0 && end.value > 0) || (start.value < 0 && end.value < 0))
	    ydirzeros[i][j][k] = null;
	  else {
	    triple root = (x,0,z);
	    root += Y * findroot(f_along_y, start.pt.y, end.pt.y,
				    fa=start.value, fb=end.value);
	    triple normal = grad(root);
	    ydirzeros[i][j][k] = positionedvector(root, normal);
	  }
	}
      }
    }

    for (int i = 0; i < nx+1; ++i) {
      for (int j = 0; j < ny+1; ++j) {
	real x = corners[i][j][0].pt.x;
	real y = corners[i][j][0].pt.y;
	real f_along_z(real t) { return f((x, y, t)); }
	for (int k = 0; k < nz; ++k) {
	  if (zdirzeros[i][j][k] != null) continue;
	  evaluatedpoint start = corners[i][j][k];
	  evaluatedpoint end = corners[i][j][k+1];
	  if ((start.value > 0 && end.value > 0) || (start.value < 0 && end.value < 0))
	    zdirzeros[i][j][k] = null;
	  else {
	    triple root = (x,y,0);
	    root += Z * findroot(f_along_z, start.pt.z, end.pt.z,
			       fa=start.value, fb=end.value);
	    triple normal = grad(root);
	    zdirzeros[i][j][k] = positionedvector(root, normal);
	  }
	}
      }
    }
  }

  // Fill in the grid vertices and the zeros along edges. Each cube starts at
  // depth one and the depth increases each time it subdivides; maxdepth is the
  // maximum subdivision depth. When a cube at maxdepth cannot be resolved to
  // patches, it is left empty.
  void operator init(int nx, int ny, int nz,
		     real f(triple), triple a, triple b,
		     int maxdepth = 6) {
    this.nx = nx;
    this.ny = ny;
    this.nz = nz;
    grad = nGrad(f);
    this.f = f;
    this.maxdepth = maxdepth;
    corners = make3dgrid(a, b, nx, ny, nz, f);
    xdirzeros = new positionedvector[nx][ny+1][nz+1];
    ydirzeros = new positionedvector[nx+1][ny][nz+1];
    zdirzeros = new positionedvector[nx+1][ny+1][nz];

    for (int i = 0; i <= nx; ++i) {
      for (int j = 0; j <= ny; ++j) {
	for (int k = 0; k <= nz; ++k) {
	  if (i < nx) xdirzeros[i][j][k] = null;
	  if (j < ny) ydirzeros[i][j][k] = null;
	  if (k < nz) zdirzeros[i][j][k] = null;
	}
      }
    }

    fillzeros();
  }

  // Doubles nx, ny, and nz by halving the sizes of the cubes along the x, y, and z
  // directions (resulting in 8 times as many cubes). Already existing data about
  // function values and zeros is copied; vertices and edges with no such pre-existing
  // data are populated.
  //
  // Returns true if subdivide succeeded, false if it failed (because maxdepth
  // was exceeded).
  bool subdivide() {
    if (maxdepth <= 1) {
      return false;
    }
    --maxdepth;
    triple a = corners[0][0][0];
    triple b = corners[nx][ny][nz];
    nx *= 2;
    ny *= 2;
    nz *= 2;
    evaluatedpoint[][][] oldcorners = corners;
    corners = new evaluatedpoint[nx+1][ny+1][nz+1];
    for (int i = 0; i <= nx; ++i) {
      for (int j = 0; j <= ny; ++j) {
	for (int k = 0; k <= nz; ++k) {
	  if (i % 2 == 0 && j % 2 == 0 && k % 2 == 0) {
	    corners[i][j][k] = oldcorners[quotient(i,2)][quotient(j,2)][quotient(k,2)];
	  } else {
	    triple pt = (interp(a.x, b.x, i/nx),
			 interp(a.y, b.y, j/ny),
			 interp(a.z, b.z, k/nz));
	    real value = f(pt);
	    if (value == 0) value = 1e-5;
	    corners[i][j][k] = evaluatedpoint(pt, value);
	  }
	}
      }
    }

    positionedvector[][][] oldxdir = xdirzeros;
    xdirzeros = new positionedvector[nx][ny+1][nz+1];
    for (int i = 0; i < nx; ++i) {
      for (int j = 0; j < ny + 1; ++j) {
	for (int k = 0; k < nz + 1; ++k) {
	  if (j % 2 != 0 || k % 2 != 0) {
	    xdirzeros[i][j][k] = null;
	  } else {
	    positionedvector zero = oldxdir[quotient(i,2)][quotient(j,2)][quotient(k,2)];
	    if (zero == null) {
	      xdirzeros[i][j][k] = null;
	      continue;
	    }
	    real x = zero.position.x;
	    if (x > interp(a.x, b.x, i/nx) && x < interp(a.x, b.x, (i+1)/nx)) {
	      xdirzeros[i][j][k] = zero;
	    } else {
	      xdirzeros[i][j][k] = null;
	    }
	  }
	}
      }
    }

    positionedvector[][][] oldydir = ydirzeros;
    ydirzeros = new positionedvector[nx+1][ny][nz+1];
    for (int i = 0; i < nx+1; ++i) {
      for (int j = 0; j < ny; ++j) {
	for (int k = 0; k < nz + 1; ++k) {
	  if (i % 2 != 0 || k % 2 != 0) {
	    ydirzeros[i][j][k] = null;
	  } else {
	    positionedvector zero = oldydir[quotient(i,2)][quotient(j,2)][quotient(k,2)];
	    if (zero == null) {
	      ydirzeros[i][j][k] = null;
	      continue;
	    }
	    real y = zero.position.y;
	    if (y > interp(a.y, b.y, j/ny) && y < interp(a.y, b.y, (j+1)/ny)) {
	      ydirzeros[i][j][k] = zero;
	    } else {
	      ydirzeros[i][j][k] = null;
	    }
	  }
	}
      }
    }

    positionedvector[][][] oldzdir = zdirzeros;
    zdirzeros = new positionedvector[nx+1][ny+1][nz];
    for (int i = 0; i < nx + 1; ++i) {
      for (int j = 0; j < ny + 1; ++j) {
	for (int k = 0; k < nz; ++k) {
	  if (i % 2 != 0 || j % 2 != 0) {
	    zdirzeros[i][j][k] = null;
	  } else {
	    positionedvector zero = oldzdir[quotient(i,2)][quotient(j,2)][quotient(k,2)];
	    if (zero == null) {
	      zdirzeros[i][j][k] = null;
	      continue;
	    }
	    real z = zero.position.z;
	    if (z > interp(a.z, b.z, k/nz) && z < interp(a.z, b.z, (k+1)/nz)) {
	      zdirzeros[i][j][k] = zero;
	    } else {
	      zdirzeros[i][j][k] = null;
	    }
	  }
	}
      }
    }

    fillzeros();
    return true;
  }

  // Forward declaration of the draw method, which will be called by drawcube().
  patch[] draw(bool[] reportactive = null);

  // Construct the patches, assuming that we are working
  // with a single cube (nx = ny = nz = 1). This method will subdivide the
  // cube if necessary. The parameter reportactive should be an array of
  // length 6. Setting an entry to true indicates that the surface abuts the
  // corresponding face (according to the earlier enum), and thus that the
  // algorithm should be sure that something is drawn in the cube sharing
  // that face--even if all the vertices of that cube have the same sign.
  patch[] drawcube(bool[] reportactive = null) {
    // First, determine which edges (if any) actually have zeros on them.
    edge[] zeroedges = new edge[0];
    positionedvector[] zeros = new positionedvector[0];

    int currentface, nextface;

    void pushifnonnull(positionedvector v) {
      if (v != null) {
	zeroedges.push(edge(currentface, nextface));
	zeros.push(v);
      }
    }
    positionedvector findzero(int face1, int face2) {
      edge e = edge(face1, face2);
      for (int i = 0; i < zeroedges.length; ++i) {
	if (zeroedges[i] == e) return zeros[i];
      }
      return null;
    }
    
    currentface = XLOW;
    nextface = YHIGH;
    pushifnonnull(zdirzeros[0][1][0]);
    nextface = YLOW;
    pushifnonnull(zdirzeros[0][0][0]);
    nextface = ZHIGH;
    pushifnonnull(ydirzeros[0][0][1]);
    nextface = ZLOW;
    pushifnonnull(ydirzeros[0][0][0]);

    currentface = XHIGH;
    nextface = YHIGH;
    pushifnonnull(zdirzeros[1][1][0]);
    nextface = YLOW;
    pushifnonnull(zdirzeros[1][0][0]);
    nextface = ZHIGH;
    pushifnonnull(ydirzeros[1][0][1]);
    nextface = ZLOW;
    pushifnonnull(ydirzeros[1][0][0]);

    currentface = YHIGH;
    nextface = ZHIGH;
    pushifnonnull(xdirzeros[0][1][1]);
    currentface = ZHIGH;
    nextface = YLOW;
    pushifnonnull(xdirzeros[0][0][1]);
    currentface = YLOW;
    nextface = ZLOW;
    pushifnonnull(xdirzeros[0][0][0]);
    currentface = ZLOW;
    nextface = YHIGH;
    pushifnonnull(xdirzeros[0][1][0]);

    //Now, string those edges together to make a circle.

    patch[] subdividecube() {
      if (!subdivide()) {
	return new patch[0];
      }
      return draw(reportactive);
    }
    if (zeroedges.length < 3) {
      return subdividecube();
    }
    int[] faceorder = makecircle(zeroedges);
    if (alias(faceorder,null)) {
      return subdividecube();
    }
    positionedvector[] patchcorners = new positionedvector[0];
    for (int i = 0; i < faceorder.length; ++i) {
      patchcorners.push(findzero(faceorder[i], faceorder[i+1]));
    }
    patchcorners.cyclic = true;

    //Now, produce the cyclic path around the edges.
    path3 edgecycle;
    for (int i = 0; i < faceorder.length; ++i) {
      path3 currentpath = pathinface(patchcorners[i], patchcorners[i+1],
					faceorder[i+1], faceorder[i],
					faceorder[i+2]);
      triple testpoint = point(currentpath, 0.5);
      if (!checkpt(testpoint, f, grad, corners[0][0][0], corners[1][1][1])) {
	return subdividecube();
      }
       
      edgecycle = edgecycle & currentpath;
    }
    edgecycle = edgecycle & cycle;

    patch[] toreturn = quadpatches(edgecycle, patchcorners, f, grad,
				     corners[0][0][0], corners[1][1][1]);
    if (alias(toreturn, null)) return subdividecube();
    return toreturn;
  }

  // Extracts the specified cube as a gridwithzeros object with
  // nx = ny = nz = 1.
  gridwithzeros getcube(int i, int j, int k) {
    gridwithzeros cube = new gridwithzeros;
    cube.grad = grad;
    cube.f = f;
    cube.nx = 1;
    cube.ny = 1;
    cube.nz = 1;
    cube.maxdepth = maxdepth;
    cube.corners = slice(corners,i,i+2,j,j+2,k,k+2);
    cube.xdirzeros = slice(xdirzeros,i,i+1,j,j+2,k,k+2);
    cube.ydirzeros = slice(ydirzeros,i,i+2,j,j+1,k,k+2);
    cube.zdirzeros = slice(zdirzeros,i,i+2,j,j+2,k,k+1);
    return cube;
  }

  // Returns an array of patches representing the surface.
  // The parameter reportactive should be an array of
  // length 6. Setting an entry to true indicates that the surface abuts the
  // corresponding face of the cube that bounds the entire grid.
  //
  // If reportactive == null, it is assumed that this is a top-level call;
  // a dot is printed to stdout for each cube drawn as a very rough
  // progress indicator.
  //
  // If reportactive != null, then it is assumed that the caller had a strong
  // reason to believe that this grid contains a part of the surface; the
  // grid will subdivide all the way to maxdepth if necessary to find points
  // on the surface.
  draw = new patch[](bool[] reportactive = null) {
    // A list of all the patches not already drawn but known
    // to contain part of the surface. This "queue" is
    // actually implemented as stack for simplicity, since
    // it does not make any difference. In a multi-threaded
    // version of the algorithm, a queue (shared across all threads)
    // would make more sense than a stack.
    triple[] queue = new triple[0];
    bool[][][] enqueued = new bool[nx][ny][nz];
    for (int i = 0; i < enqueued.length; ++i) {
      for (int j = 0; j < enqueued[i].length; ++j) {
	for (int k = 0; k < enqueued[i][j].length; ++k) {
	  enqueued[i][j][k] = false;
	}
      }
    }

    void enqueue(int i, int j, int k) {
      if (i >= 0 && i < nx
	  && j >= 0 && j < ny
	  && k >= 0 && k < nz
	  && !enqueued[i][j][k]) {
	queue.push((i,j,k));
	enqueued[i][j][k] = true;
      }
      if (!alias(reportactive, null)) {
	if (i < 0) reportactive[XLOW] = true;
	if (i >= nx) reportactive[XHIGH] = true;
	if (j < 0) reportactive[YLOW] = true;
	if (j >= ny) reportactive[YHIGH] = true;
	if (k < 0) reportactive[ZLOW] = true;
	if (k >= nz) reportactive[ZHIGH] = true;
      }
    }
    
    for (int i = 0; i < nx+1; ++i) {
      for (int j = 0; j < ny+1; ++j) {
	for (int k = 0; k < nz+1; ++k) {
	  if (i < nx && xdirzeros[i][j][k] != null) {
	    for (int jj = j-1; jj <= j; ++jj)
	      for (int kk = k-1; kk <= k; ++kk)
		enqueue(i, jj, kk);
	  }
	  if (j < ny && ydirzeros[i][j][k] != null) {
	    for (int ii = i-1; ii <= i; ++ii)
	      for (int kk = k-1; kk <= k; ++kk)
		enqueue(ii, j, kk);
	  }
	  if (k < nz && zdirzeros[i][j][k] != null) {
	    for (int ii = i-1; ii <= i; ++ii)
	      for (int jj = j-1; jj <= j; ++jj)
		enqueue(ii, jj, k);
	  }
	}
      }
    }

    if (!alias(reportactive, null) && queue.length == 0) {
      if (subdivide()) return draw(reportactive);
    }

    patch[] surface = new patch[0];

    while (queue.length > 0) {
      triple coord = queue.pop();
      int i = floor(coord.x);
      int j = floor(coord.y);
      int k = floor(coord.z);
      bool[] reportface = array(6, false);
      patch[] toappend = getcube(i,j,k).drawcube(reportface);
      if (reportface[XLOW]) enqueue(i-1,j,k);
      if (reportface[XHIGH]) enqueue(i+1,j,k);
      if (reportface[YLOW]) enqueue(i,j-1,k);
      if (reportface[YHIGH]) enqueue(i,j+1,k);
      if (reportface[ZLOW]) enqueue(i,j,k-1);
      if (reportface[ZHIGH]) enqueue(i,j,k+1);
      surface.append(toappend);
      if (settings.verbose > 1 && alias(reportactive, null)) write(stdout, '.');
    }
    if (settings.verbose > 1 && alias(reportactive, null)) write(stdout, '\n');
    return surface;
  };
}

// The external interface of this whole module. Accepts exactly one 
// function (throws an error if two or zero functions are specified).
// The function should be differentiable. (Whatever you do, do not
// pass in an indicator function!) Ideally, the zero locus of the
// function should be smooth; singularities will significantly slow
// down the algorithm and potentially give bad results.
//
// Returns a plot of the zero locus of the function within the
// rectangular solid with opposite corners at a and b.
//
// Additional parameters:
// n - the number of initial segments in each of the x, y, z directions.
// overlapedges - if true, the patches of the surface are slightly enlarged
//     to compensate for an artifact in which the viewer can see through the
//     boundary between patches. (Some of this may actually be a result of
//     edges not lining up perfectly, but I'm fairly sure a lot of it arises
//     purely as a rendering artifact.)
// nx - override n in the x direction
// ny - override n in the y direction
// nz - override n in the z direction
// maxdepth - the maximum depth to which the algorithm will subdivide in
//     an effort to find patches that closely approximate the true surface.
surface implicitsurface(real f(triple) = null, real ff(real,real,real) = null,
			 triple a, triple b,
			 int n = nmesh,
			 bool keyword overlapedges = false,
			 int keyword nx=n, int keyword ny=n,
			 int keyword nz=n,
			 int keyword maxdepth = 8) {
  if (f == null && ff == null)
    abort("implicitsurface called without specifying a function.");
  if (f != null && ff != null)
    abort("Only specify one function when calling implicitsurface.");
  if (f == null) f = new real(triple w) { return ff(w.x, w.y, w.z); };
  gridwithzeros grid = gridwithzeros(nx, ny, nz, f, a, b, maxdepth=maxdepth);
  patch[] patches = grid.draw();
  if (overlapedges) {
    for (int i = 0; i < patches.length; ++i) {
      triple center = patches[i].point(1/2,1/2);
      patches[i] = shift(center) * scale3(1.01) * shift(-center) * patches[i];
    }
  }
  return surface(...patches);
}