Asymptote: Is there an equivalent to silhouette for surfaces?

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);

render render=render(compression=Low,merge=true);
draw(s,render);
draw(zscale3(-1)*s); • I'm pretty sure there's no easy way to do this, otherwise the solids package would not need to define a special silhouette function for objects of type revolution. In fact, even the function in solids is not bug-free; you might check and see if it can handle your surface. The dimple may prove especially challenging. – Charles Staats Sep 26 '13 at 21:23
• Philipp, I think you would lose too much information by projecting your surface onto a plane. For instance, if you tried this with the hyperbola example, then it would no longer be possible to tell which of the two circles was supposed to be visible. What you really want to do is to put together a path3[] tracing out all the points at which the projection vector is tangent to the surface, together with all the "edge paths" of the surface (excluding those at which two edges are meeting to form a smooth surface). Then, add in the surface drawn with e.g. surfacepen=emissive(white),... – Charles Staats Sep 27 '13 at 19:49
• ... which draws the surface in pure white. (Replace white with whatever your background color is.) The purpose of this is to hide from view the 3d paths that are behind the surface. The most difficult part of this will probably be to form paths tracing out all the points at which the surface is tangent to the projection vector (or equivalently, the normal to the surface is normal to the projection vector). Good luck! – Charles Staats Sep 27 '13 at 19:53
• Philip, I'm sure it's possible, but I would not know how to do it offhand. A surface is really a one-dimensional array of Bezier patches, so you want to look in three_surface.asy for the definitions of struct patch and struct surface to see what information is actually available to you. (Just google "three_surface.asy" to find the file. Use a search tool to find all three occurrences of struct in the file.) – Charles Staats Sep 27 '13 at 23:54
• I've managed to provide a partial implementation of my suggestions that at least works for a torus; see tex.stackexchange.com/a/135438/484 – Charles Staats Sep 28 '13 at 5:20

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.length;
return s.index.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);
}
}

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

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 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]);
}
}

/*
* 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)));
}
}

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: 