diff options
Diffstat (limited to 'Master/texmf-dist/asymptote/smoothcontour3.asy')
-rw-r--r-- | Master/texmf-dist/asymptote/smoothcontour3.asy | 748 |
1 files changed, 383 insertions, 365 deletions
diff --git a/Master/texmf-dist/asymptote/smoothcontour3.asy b/Master/texmf-dist/asymptote/smoothcontour3.asy index 216f3269d92..4fe26654283 100644 --- a/Master/texmf-dist/asymptote/smoothcontour3.asy +++ b/Master/texmf-dist/asymptote/smoothcontour3.asy @@ -19,69 +19,7 @@ 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; -} +import math; /***********************************************/ /******** CREATING BEZIER PATCHES **************/ @@ -126,7 +64,7 @@ function[] bernstein = new function[] {B03, B13, B23, B33}; // the specified boundary path. However, the entries in the array // remain intact. patch patchwithnormals(path3 external, triple[] u0normals, triple[] u1normals, - triple[] v0normals, triple[] v1normals) + triple[] v0normals, triple[] v1normals) { assert(cyclic(external)); assert(length(external) == 4); @@ -275,8 +213,8 @@ patch patchwithnormals(path3 external, triple[] u0normals, triple[] u1normals, 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."); + write("Warning: unable to solve matrix for specifying edge normals " + + "on bezier patch. Using coons patch."); return patch(external); } @@ -284,22 +222,95 @@ patch patchwithnormals(path3 external, triple[] u0normals, triple[] u1normals, 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]); + solution[position+1], + solution[position+2]); } } return patch(controlpoints); } -// A wrapper for the previous function when the normal direction +// This function attempts to produce a Bezier triangle +// with the specified boundary path and normal directions at the +// edge midpoints. The bezier triangle should be normal to +// n1 at point(external, 0.5), +// normal to n2 at point(external, 1.5), and +// normal to n3 at point(external, 2.5). +// 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). +patch trianglewithnormals(path3 external, triple n1, + triple n2, triple n3) { + assert(cyclic(external)); + assert(length(external) == 3); + // Use the formal symbols a3, a2b, abc, etc. to denote the control points, + // following the Wikipedia article on Bezier triangles. + triple a3 = point(external, 0), a2b = postcontrol(external, 0), + ab2 = precontrol(external, 1), b3 = point(external, 1), + b2c = postcontrol(external, 1), bc2 = precontrol(external, 2), + c3 = point(external, 2), ac2 = postcontrol(external, 2), + a2c = precontrol(external, 0); + + // Use orthogonal projection to ensure that the normal vectors are + // actually normal to the boundary path. + triple tangent = dir(external, 0.5); + n1 -= dot(n1,tangent)*tangent; + n1 = unit(n1); + + tangent = dir(external, 1.5); + n2 -= dot(n2,tangent)*tangent; + n2 = unit(n2); + + tangent = dir(external, 2.5); + n3 -= dot(n3,tangent)*tangent; + n3 = unit(n3); + + real wild = 2 * wildnessweight; + real[][] matrix = { {n1.x, n1.y, n1.z}, + {n2.x, n2.y, n2.z}, + {n3.x, n3.y, n3.z}, + { wild, 0, 0}, + { 0, wild, 0}, + { 0, 0, wild} }; + real[] rightvector = + { dot(n1, (a3 + 3a2b + 3ab2 + b3 - 2a2c - 2b2c)) / 4, + dot(n2, (b3 + 3b2c + 3bc2 + c3 - 2ab2 - 2ac2)) / 4, + dot(n3, (c3 + 3ac2 + 3a2c + a3 - 2bc2 - 2a2b)) / 4 }; + + // The inner control point that minimizes the sum of squares of + // the mixed partials on the corners. + triple tameinnercontrol = + ((a2b + a2c - a3) + (ab2 + b2c - b3) + (ac2 + bc2 - c3)) / 3; + rightvector.append(wild * new real[] + {tameinnercontrol.x, tameinnercontrol.y, tameinnercontrol.z}); + real[] solution = leastsquares(matrix, rightvector, warn=false); + if (solution.length == 0) { // if the matrix was singular + write("Warning: unable to solve matrix for specifying edge normals " + + "on bezier triangle. Using coons triangle."); + return patch(external); + } + triple innercontrol = (solution[0], solution[1], solution[2]); + return patch(external, innercontrol); +} + +// A wrapper for the previous functions 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. +// are three or four segments. patch patchwithnormals(path3 external, triple normalat(triple)) { assert(cyclic(external)); assert(1 <= length(external) && length(external) <= 4); + if (length(external) == 3) { + triple n1 = normalat(point(external, 0.5)); + triple n2 = normalat(point(external, 1.5)); + triple n3 = normalat(point(external, 2.5)); + return trianglewithnormals(external, n1, n2, n3); + } while (length(external) < 4) external = external -- cycle; triple[] u0normals = new triple[3]; triple[] u1normals = new triple[3]; @@ -315,82 +326,6 @@ patch patchwithnormals(path3 external, triple normalat(triple)) { } /***********************************************/ -/*********** 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 *************/ /***********************************************/ @@ -466,8 +401,8 @@ struct int_to_intset { void add(int key, int value) { for (int i = 0; i < keys.length; ++i) { if (keys[i] == key) { - values[i].add(value); - return; + values[i].add(value); + return; } } keys.push(key); @@ -561,9 +496,9 @@ int[] makecircle(edge[] edges) { if (adjacentvertices.length != 2) return null; for (int v : adjacentvertices) { if (v != lastvertex) { - lastvertex = currentvertex; - currentvertex = v; - break; + lastvertex = currentvertex; + currentvertex = v; + break; } } } while (currentvertex != startvertex); @@ -597,9 +532,10 @@ string operator cast(positionedvector vv) { // The angle, in degrees, between two vectors. real angledegrees(triple a, triple b) { - real lengthprod = abs(a) * abs(b); + real dotprod = dot(a,b); + real lengthprod = max(abs(a) * abs(b), abs(dotprod)); if (lengthprod == 0) return 0; - return aCos(dot(a,b) / lengthprod); + return aCos(dotprod / lengthprod); } // A path (single curved segment) between two points. At each point @@ -622,10 +558,10 @@ path3 pathbetween(positionedvector v1, positionedvector v2) { // 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) { + bool warn=false) { real[][] matrix = {{v1.x, v2.x}, - {v1.y, v2.y}, - {v1.z, v2.z}}; + {v1.y, v2.y}, + {v1.z, v2.z}}; real[] desiredanswer = {toproject.x, toproject.y, toproject.z}; return leastsquares(matrix, desiredanswer, warn=warn); } @@ -635,7 +571,7 @@ real[] projecttospan_findcoeffs(triple toproject, triple v1, triple v2, // 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 mincoeff = 0.05) { real[] coeffs = projecttospan_findcoeffs(toproject, v1, v2, warn=false); real a, b; if (coeffs.length == 0) { @@ -701,25 +637,25 @@ path3 bisector(path3 edgecycle, int[] savevertices) { int opposite = i + 3; triple vec = unit(point(edgecycle, opposite) - point(edgecycle, i)); real[] coeffsbegin = projecttospan_findcoeffs(vec, forwarddirections[i], - backwarddirections[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]); + 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]; + * coeffsbegin[1] * coeffsend[1]; if (goodness > maxgoodness) { maxgoodness = goodness; directionout = coeffsbegin[0] * forwarddirections[i] + - coeffsbegin[1] * backwarddirections[i]; + coeffsbegin[1] * backwarddirections[i]; directionin = -(coeffsend[0] * forwarddirections[opposite] + - coeffsend[1] * backwarddirections[opposite]); + coeffsend[1] * backwarddirections[opposite]); chosenindex = i; } } @@ -738,7 +674,7 @@ path3 bisector(path3 edgecycle, int[] savevertices) { // 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 facenorm, triple edge1normout, triple edge2normout) { triple dir1 = cross(v1.direction, facenorm); real dotprod = dot(dir1, edge1normout); @@ -769,15 +705,16 @@ triple normalout(int face) { // 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) { + int face, int edge1face, int edge2face) { return pathinface(v1, v2, normalout(face), normalout(edge1face), - normalout(edge2face)); + normalout(edge2face)); } /***********************************************/ /******** DRAWING IMPLICIT SURFACES ************/ /***********************************************/ +// DEPRECATED // Quadrilateralization: // Produce a surface (array of *nondegenerate* Bezier patches) with a // specified three-segment boundary. The surface should approximate the @@ -787,11 +724,12 @@ path3 pathinface(positionedvector v1, positionedvector v2, // 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 +// the usual trend in mathematics. However, *before the introduction of bezier +// triangles,* the pathwithnormals algorithm +// did a poor job of choosing a good surface when the boundary path did // not consist of four positive-length segments. patch[] triangletoquads(path3 external, real f(triple), triple grad(triple), - triple a, triple b) { + triple a, triple b) { static real epsilon = 1e-3; assert(length(external) == 3); assert(cyclic(external)); @@ -867,6 +805,65 @@ patch[] triangletoquads(path3 external, real f(triple), triple grad(triple), patchwithnormals(quad2, grad)}; } +// Attempts to fill the path external (which should by a cyclic path consisting of +// three segments) with bezier triangle(s). Returns an empty array if it fails. +// +// In more detail: A single bezier triangle is computed using trianglewithnormals. The normals of +// the resulting triangle at the midpoint of each edge are computed. If any of these normals +// is in the negative f direction, the external triangle is subdivided into four external triangles +// and the same procedure is applied to each. If one or more of them has an incorrectly oriented +// edge normal, the function gives up and returns an empty array. +// +// Thus, the returned array consists of 0, 1, or 4 bezier triangles; no other array lengths +// are possible. +// +// This function assumes that the path orientation is consistent with f (and its gradient) +// -- i.e., that +// at a corner, (tangent in) x (tangent out) is in the positive f direction. +patch[] maketriangle(path3 external, real f(triple), + triple grad(triple), bool allowsubdivide = true) { + assert(cyclic(external)); + assert(length(external) == 3); + triple m1 = point(external, 0.5); + triple n1 = unit(grad(m1)); + triple m2 = point(external, 1.5); + triple n2 = unit(grad(m2)); + triple m3 = point(external, 2.5); + triple n3 = unit(grad(m3)); + patch beziertriangle = trianglewithnormals(external, n1, n2, n3); + if (dot(n1, beziertriangle.normal(0.5, 0)) >= 0 && + dot(n2, beziertriangle.normal(0.5, 0.5)) >= 0 && + dot(n3, beziertriangle.normal(0, 0.5)) >= 0) + return new patch[] {beziertriangle}; + + if (!allowsubdivide) return new patch[0]; + + positionedvector m1 = positionedvector(m1, n1); + positionedvector m2 = positionedvector(m2, n2); + positionedvector m3 = positionedvector(m3, n3); + path3 p12 = pathbetween(m1, m2); + path3 p23 = pathbetween(m2, m3); + path3 p31 = pathbetween(m3, m1); + patch[] triangles = maketriangle(p12 & p23 & p31 & cycle, f, grad=grad, + allowsubdivide=false); + if (triangles.length < 1) return new patch[0]; + + triangles.append(maketriangle(subpath(external, -0.5, 0.5) & reverse(p31) & cycle, + f, grad=grad, allowsubdivide=false)); + if (triangles.length < 2) return new patch[0]; + + triangles.append(maketriangle(subpath(external, 0.5, 1.5) & reverse(p12) & cycle, + f, grad=grad, allowsubdivide=false)); + if (triangles.length < 3) return new patch[0]; + + triangles.append(maketriangle(subpath(external, 1.5, 2.5) & reverse(p23) & cycle, + f, grad=grad, allowsubdivide=false)); + if (triangles.length < 4) return new patch[0]; + + return triangles; +} + + // 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). @@ -891,14 +888,14 @@ bool checkptincube(triple pt, triple a, triple b) { 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)); + (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) { + triple a, triple b) { return checkptincube(testpt, a, b) && check_fpt_zero(testpt, f, grad); } @@ -910,8 +907,8 @@ bool checkpt(triple testpt, real f(triple), triple grad(triple), // 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) { + real f(triple), triple grad(triple), + triple a, triple b, bool usetriangles) { assert(corners.cyclic); // The tolerance for considering two points "essentially identical." @@ -924,8 +921,8 @@ patch[] quadpatches(path3 edgecycle, positionedvector[] corners, if (corners.length == 2) return new patch[0]; corners.delete(i); edgecycle = subpath(edgecycle, 0, i) - & subpath(edgecycle, i+1, length(edgecycle)) - & cycle; + & subpath(edgecycle, i+1, length(edgecycle)) + & cycle; --i; assert(length(edgecycle) == corners.length); } @@ -937,9 +934,9 @@ patch[] quadpatches(path3 edgecycle, positionedvector[] corners, 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))); + 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))); + dir(edgecycle, 1, sign=1, normalize=false))); if (area0 < areatolerance && area1 < areatolerance) return new patch[0]; else return null; } @@ -947,13 +944,14 @@ patch[] quadpatches(path3 edgecycle, positionedvector[] corners, for (int i = 0; i < length(edgecycle); ++i) { if (angledegrees(dir(edgecycle,i,sign=1), - dir(edgecycle,i+1,sign=-1)) > 80) { + dir(edgecycle,i+1,sign=-1)) > 80) { return null; } } if (length(edgecycle) == 3) { - patch[] toreturn = triangletoquads(edgecycle, f, grad, a, b); + patch[] toreturn = usetriangles ? maketriangle(edgecycle, f, grad) + : triangletoquads(edgecycle, f, grad, a, b); if (toreturn.length == 0) return null; else return toreturn; } @@ -974,8 +972,9 @@ patch[] quadpatches(path3 edgecycle, positionedvector[] corners, & 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); + & subpath(edgecycle, bisectorindices[1], 5+bisectorindices[0]) & cycle; + toreturn = usetriangles ? maketriangle(secondpatch, f, grad) + : triangletoquads(secondpatch, f, grad, a, b); if (toreturn.length == 0) return null; toreturn.push(patchwithnormals(firstpatch, grad)); } else { @@ -984,7 +983,7 @@ patch[] quadpatches(path3 edgecycle, positionedvector[] corners, & subpath(edgecycle, bisectorindices[1], 6+bisectorindices[0]) & cycle; toreturn = new patch[] {patchwithnormals(firstpatch, grad), - patchwithnormals(secondpatch, grad)}; + patchwithnormals(secondpatch, grad)}; } return toreturn; } @@ -995,8 +994,8 @@ 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) ); + (f(v + epsilon*Y) - f(v - epsilon*Y)) / (2 epsilon), + (f(v + epsilon*Z) - f(v - epsilon*Z)) / (2 epsilon) ); }; } @@ -1015,18 +1014,18 @@ 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) + 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); + 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); } } } @@ -1045,8 +1044,8 @@ T[][] slice(T[][] a, int start1, int end1, int start2, int end2) { return toreturn; } T[][][] slice(T[][][] a, int start1, int end1, - int start2, int end2, - int start3, int end3) { + 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); @@ -1062,8 +1061,8 @@ T[][] slice(T[][] a, int start1, int end1, int start2, int end2) { return toreturn; } T[][][] slice(T[][][] a, int start1, int end1, - int start2, int end2, - int start3, int end3) { + 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); @@ -1083,73 +1082,74 @@ struct gridwithzeros { triple grad(triple); real f(triple); int maxdepth; + bool usetriangles; // 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); - } - } + 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); - } - } + 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); - } - } + 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); + } + } } } } @@ -1159,14 +1159,15 @@ struct gridwithzeros { // 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) { + real f(triple), triple a, triple b, + int maxdepth = 6, bool usetriangles) { this.nx = nx; this.ny = ny; this.nz = nz; grad = nGrad(f); this.f = f; this.maxdepth = maxdepth; + this.usetriangles = usetriangles; corners = make3dgrid(a, b, nx, ny, nz, f); xdirzeros = new positionedvector[nx][ny+1][nz+1]; ydirzeros = new positionedvector[nx+1][ny][nz+1]; @@ -1174,11 +1175,11 @@ struct gridwithzeros { 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; - } + 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; + } } } @@ -1206,18 +1207,18 @@ struct gridwithzeros { 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); - } - } + 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); + } + } } } @@ -1225,23 +1226,23 @@ struct gridwithzeros { 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; - } - } - } + 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; + } + } + } } } @@ -1249,23 +1250,23 @@ struct gridwithzeros { 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; - } - } - } + 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; + } + } + } } } @@ -1273,23 +1274,23 @@ struct gridwithzeros { 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; - } - } - } + 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; + } + } + } } } @@ -1316,14 +1317,14 @@ struct gridwithzeros { void pushifnonnull(positionedvector v) { if (v != null) { - zeroedges.push(edge(currentface, nextface)); - zeros.push(v); + 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]; + if (zeroedges[i] == e) return zeros[i]; } return null; } @@ -1365,7 +1366,7 @@ struct gridwithzeros { patch[] subdividecube() { if (!subdivide()) { - return new patch[0]; + return new patch[0]; } return draw(reportactive); } @@ -1386,19 +1387,31 @@ struct gridwithzeros { 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]); + 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(); + return subdividecube(); } edgecycle = edgecycle & currentpath; } edgecycle = edgecycle & cycle; + + { // Ensure the outward normals are pointing in the same direction as the gradient. + triple tangentin = patchcorners[0].position - precontrol(edgecycle, 0); + triple tangentout = postcontrol(edgecycle, 0) - patchcorners[0].position; + triple normal = cross(tangentin, tangentout); + if (dot(normal, patchcorners[0].direction) < 0) { + edgecycle = reverse(edgecycle); + patchcorners = patchcorners[-sequence(patchcorners.length)]; + patchcorners.cyclic = true; + } + } + patch[] toreturn = quadpatches(edgecycle, patchcorners, f, grad, - corners[0][0][0], corners[1][1][1]); + corners[0][0][0], corners[1][1][1], usetriangles); if (alias(toreturn, null)) return subdividecube(); return toreturn; } @@ -1413,6 +1426,7 @@ struct gridwithzeros { cube.ny = 1; cube.nz = 1; cube.maxdepth = maxdepth; + cube.usetriangles = usetriangles; 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); @@ -1434,6 +1448,7 @@ struct gridwithzeros { // grid will subdivide all the way to maxdepth if necessary to find points // on the surface. draw = new patch[](bool[] reportactive = null) { + if (alias(reportactive, null)) progress(true); // 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 @@ -1444,49 +1459,49 @@ struct gridwithzeros { 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; - } + 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; + && 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; + 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); - } - } + 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); + } + } } } @@ -1510,9 +1525,9 @@ struct gridwithzeros { 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 (alias(reportactive, null)) progress(); } - if (settings.verbose > 1 && alias(reportactive, null)) write(stdout, '\n'); + if (alias(reportactive, null)) progress(false); return surface; }; } @@ -1540,22 +1555,25 @@ struct gridwithzeros { // 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) { + 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, + bool keyword usetriangles=true) { 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); + gridwithzeros grid = gridwithzeros(nx, ny, nz, f, a, b, maxdepth=maxdepth, + usetriangles=usetriangles); patch[] patches = grid.draw(); if (overlapedges) { for (int i = 0; i < patches.length; ++i) { - triple center = patches[i].point(1/2,1/2); + triple center = (patches[i].triangular ? + patches[i].point(1/3, 1/3) : patches[i].point(1/2,1/2)); patches[i] = shift(center) * scale3(1.01) * shift(-center) * patches[i]; } } |