Asymptote: Is there an equivalent to silhouette for surfaces?

Question

Update

If there is no equivalent to silhouette could you comment on my following ideas: I can get the projection of my surface using the planeproject transformation. This is still a three dimension surface but maybe it is easier to extract the outline path or paths from this simpler surface? If this isn't possible neither maybe the following algorithm could work:

1. planeproject the surface
2. extract a dense mesh of points lying on this transformed surface (these points should ly within the same plane); or maybe it is even possible to extract the border points of the flat surface somehow (maybe via their slope being different from the other points on flat part of the surface)
3. use some concave hull algorithm (that would have to be implemented first, of course) to find the concave hull of the aforementioned points.
4. use the path along the concave hull as outline

Here I'm particularly uncertain if there is an easy way to do step 2. And I will have to think about what to do when the concave hull consists of several paths, e.g. like for a torus where you get the outer outline and the hole in the middle. Do you think this would be possible?

Original Question

The solids package defines revolution objects whose outlines can be accessed via silhouette. Example:

size(200);
import solids;
settings.render=0;
settings.prc=false;

currentprojection=perspective(4,4,3);
revolution hyperboloid=revolution(new real(real x) {return sqrt(1+x*x);},
                              -2,2,20,operator..,X);
draw(hyperboloid.silhouette(64),blue);

Is there an equivalent way to get the outline path of a surface? For example the outline of such a surface (I know that the object in this example is spherically symmetric and so I could define a revolution that describes it but I would like to know if there is a general way to get the outline of a surface):

import graph3;
import palette;
size(200);
currentprojection=orthographic(6,8,2);
viewportmargin=(1cm,0);

real c0=0.1;

real f(real r) {return r*(1-r/6)*exp(-r/3);}

triple f(pair t) {
  real r=t.x;
  real phi=t.y;
  real f=f(r);
  real s=max(min(c0/f,1),-1);
  real R=r*sqrt(1-s^2);
  return (R*cos(phi),R*sin(phi),r*s);
}

bool cond(pair t) {return f(t.x) != 0;}

real R=abs((20,20,20));
surface s=surface(f,(0,0),(R,2pi),100,8,Spline,cond);

s.colors(palette(s.map(abs),Gradient(palegreen,heavyblue)));

render render=render(compression=Low,merge=true);
draw(s,render);
draw(zscale3(-1)*s);

Answer 1

Okay, I think I've more or less got a solution.

The silhouette function is defined by the following code:

import graph3;
import contour;

// A bunch of auxiliary functions.

real fuzz = .001;

real umin(surface s) { return 0; }
real vmin(surface s) { return 0; }
pair uvmin(surface s) { return (umin(s), vmin(s)); }
real umax(surface s, real fuzz=fuzz) {
  if (s.ucyclic()) return s.index.length;
  else return s.index.length - fuzz;
}
real vmax(surface s, real fuzz=fuzz) {
  if (s.vcyclic()) return s.index[0].length;
  return s.index[0].length - fuzz;
}
pair uvmax(surface s, real fuzz=fuzz) { return (umax(s,fuzz), vmax(s,fuzz)); }

typedef real function(real, real);

function normalDot(surface s, triple eyedir) {
  real toreturn(real u, real v) {
    return dot(s.normal(u, v), eyedir);
  }
  return toreturn;
}

guide[] normalpathuv(surface s, triple eyedir, int n = ngraph) {
  return contour(normalDot(s, eyedir), uvmin(s), uvmax(s), new real[] {0}, nx=n)[0];
}

path3 onSurface(surface s, path p) {
  triple f(real t) {
    pair point = point(p,t);
    return s.point(point.x, point.y);
  }
  if (cyclic(p)) {
    guide3 toreturn = f(0);
    for (int i = 1; i < size(p); ++i)
      toreturn = toreturn -- f(i);
    toreturn = toreturn -- cycle;
    return toreturn;
  }
  return graph(f, 0, length(p));
}

/*
 * This method returns an array of paths that trace out all the
 * points on s at which s is parallel to eyedir.
 */
path3[] silhouetteNoEdges(surface s, triple eyedir, int n = ngraph) {
  guide[] uvpaths = normalpathuv(s, eyedir, n);
  path3[] toreturn = new path3[uvpaths.length];
  for (int i = 0; i < uvpaths.length; ++i) {
    toreturn[i] = onSurface(s, uvpaths[i]);
  }
  return toreturn;
}

/*
 * Now, add in the edges (if there are any).
 */
path3[] silhouette(surface s, triple eyedir, int n = ngraph) {
  path3[] toreturn = silhouetteNoEdges(s, eyedir, n);
  if (!s.ucyclic()) {
    toreturn.push(s.uequals(umin(s)));
    toreturn.push(s.uequals(umax(s)));
  }
  if (!s.vcyclic()) {
    toreturn.push(s.vequals(vmin(s)));
    toreturn.push(s.vequals(vmax(s)));
  }
  return toreturn;
}

After saving the code above in a file called silhouette.asy, here's how you draw a silhouette of your surface:

settings.outformat="pdf";
int resolutionfactor = 4;
settings.render=2.resolutionfactor;
settings.prc=false;

import silhouette;

size(200);
triple eye = (6,8,2);
currentprojection=orthographic(eye);
viewportmargin=(1cm,0);

real c0=0.1;

real f(real r) {return r*(1-r/6)*exp(-r/3);}

/* Note that the function below has been modified so as not to throw a divide by zero error
 * when f(r) == 0.
 */
triple f(pair t) {
  real r=t.x;
  real phi=t.y;
  real f=f(r);
  real s;
  //This assumes c0 > 0
  if (0 > f && f > -c0) s = -1;
  else if (0 <= f && f < c0) s = 1;
  else s = c0/f;
  //real s=max(min(c0/f,1),-1);
  real R=r*sqrt(1-s^2);
  return (R*cos(phi),R*sin(phi),r*s);
}

bool cond(pair t) {return f(t.x) != 0;}

real R=abs((20,20,20));
surface s=surface(f,(0,0),(R,2pi),nu=100,nv=8,Spline);

draw(silhouette(s,eye,n=200));
draw(s,surfacepen=emissive(white));

surface s2 = zscale3(-1)*s;
draw(silhouette(s2, eye, n=200));
draw(s2, surfacepen=emissive(white));

shipout(scale(resolutionfactor)*currentpicture.fit());

Here's the result: