# Asymptote: draw surface from data points in 3D

In the figure below, you can see 25 curves. Each of these curves was plotted from a list of points in a file mypoints_i.dat. Each of these files contains the exact same number of points.

How would you draw the surface described by the curves? A solution was proposed in this answer with a triangulation technique, which is not required here since all the curves are drawn from the same number of points. Additionally,

• the color of the surface should shade from blue to red (as the curves below);
• the surface should be as smooth as possible.

EDIT I manage to plot the surface but I am not able to use the smoothening operator operator.., so the surface is very rough, see below. This is because I plot the plane quadrangles one by one. Of course I could refine the points but this would result in a heavy file, and resource-consuming 3D manipulations.

There must be a way to plot the surface smoothly, right?

Note that the surface actually corresponds to (half of) the curves above, even though it is not obvious.

The code for the surface is the following:

//////// SURFACE PART I ////////

int increment=1;
for (int i=1; i<51-increment; i=i+increment)
{
string filename1 = filebasename + string(i) + "PartI.dat";
string filename2 = filebasename + string(i+increment) + "PartI.dat";
file in1=input(filename1).line().csv();
file in2=input(filename2).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
real[][] a2=in2.dimension(0,0);
a1=transpose(a1);
a2=transpose(a2);

real[] x=a1[0];
real[] y=a1[1];
real[] z=a1[3];

real[] x2=a2[0];
real[] y2=a2[1];
real[] z2=a2[3];

pen orbitpen=.7bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50));
pen projpen=.3bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50))+opacity(0.3);

int step=1;
for (int j=0; j<21-step; j=j+step)
{
path3 orbit1=graph(x,y,z,operator--);
path3 orbit2=reverse(graph(x2,y2,z2));

triple pointA1=(x[j],y[j],z[j]);
triple pointA2=(x[j+step],y[j+step],z[j+step]);
triple pointB1=(x2[j],y2[j],z2[j]);
triple pointB2=(x2[j+step],y2[j+step],z2[j+step]);

draw(surface(pointA1..pointA2--pointB2..pointB1--cycle),red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50))+opacity(1),currentlight);
}
}

EDIT2 Complete standalone example, including data files, available here.

• Could you sort all your points from farthest to closest to the POV? Jan 29, 2015 at 23:16
• @JohnKormylo: I would probably be able to do so for one POV, but I want to be able to manipulate my 3D object, in 3D within a pdf file (with prc). What do you have in mind? (BTW, I'm working on it, I add some elements to my post in a few hours.) Jan 29, 2015 at 23:30
• Never mind. Smoothing is the issue now. Jan 30, 2015 at 2:51
• See tex.stackexchange.com/a/164509/484 for an Asymptote solution to a similar problem. Jan 31, 2015 at 20:38
• @anderstood: If s was a surface created by graphing a function or a parametric function, then s.point gives a parametrization of the surface. For instance, s.point((0,0)) should give one of the four corners; and if the patches form (say) an 8 x 8 grid, then the opposite corner is given (roughly) by s.point((7.9999, 7.9999)). However, if the surface is constructed by some means other than graphing function(s), then s.point returns nonsense. Feb 3, 2015 at 3:59

I thought about it some more and managed a considerable simplification. Both this solution and the previous operate on the same principle suggested by the user Symbol 1: If you look under the hood at the operation of the smooth surface grapher, it actually plots only those points on a rectangular grid in uv-space. So, even if you have only those points (say, given in a matrix of triples), you can "trick" the grapher into believing you have a complete function--it will never know the difference, because it will only ever evaluate the function at the grid points. I believe this is also what O.G.'s "short solution" does.

Much of my previous solution was designed to figure out how to combine "rectangular" patches given in an arbitrary order. But since the points here are given in a regular order, that was overkill. (I mostly wrote it that way since I had already written the code, and the OP had already arranged the grid in rectangular patches in his/her own code.)

A note on the color palette: The draw method, when passed an array of pens, will arrange the "rectangular" patches in a regular (predictable) order and draw the first patch with the first pen, the second patch with the second pen, etc. So, technically, what is being drawn in this picture is not a smooth gradient; but in this particular example, it's close enough. By contrast, the surface.colors() method used in the answers of Symbol 1 and of O.G. colors the vertices rather than the patches (allowing for a smooth gradient, although not in prc) and stores this information in the surface object itself.

// Global Asymptote definitions can be put here.
import three;
import grid3;

usepackage("mathptmx");

// One can globally override the default toolbar settings here:
// settings.toolbar=true;

import graph3;

real xmin=-2, xmax=2;
real ymin=-2, ymax=1.5;
real zmin=-2.5, zmax=2.5;

limits((xmin,ymin,zmin),(xmax,ymax,zmax));

currentprojection=perspective(camera=(1.5,2,2.5));

unitsize(3cm,3cm,2cm);

real linewidth=1.1;
real linewidthprojections=.15;

string filebasename="./data/ForASY_mnl2ippGridBranche3T1_";

bool renderPRC = false;

if(renderPRC) {
// PRC TRUE
settings.prc=true;
settings.embed=true;
}
else {
// RASTERIZE
settings.outformat="png";
settings.prc=false;
settings.render=3;
}

/////// ORBITS IN 3D, SMOOTH ////////
for (int i=1; i<=50; i=i+2)
{
string filename1 = filebasename + string(i) + "PartI.dat";
file in1=input(filename1).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
a1=transpose(a1);

real[] x1=a1[0];
real[] y1=a1[1];
real[] z1=a1[3];

pen orbitpen=.7bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50));
pen projpen=.3bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50))+opacity(0.3);

path3 thepath1 = graph(x1, y1, z1, operator..);
draw(thepath1,orbitpen,currentlight);
}

//////// SURFACE PART I ////////

surface smoothSurface(triple[][] points) {
int nu = points.length - 1;
if (nu <= 0) abort("Grid must have at least two rows to produce a surface.");
int nv = points[0].length - 1;
if (nv <= 0) abort("Grid must have at least two columns to produce a surface.");
// Create a parametric function that is designed specifically for integer values.
triple f(pair uv) { return points[floor(uv.x)][floor(uv.y)]; }
// Now graph that parametric function:
return surface(f, (0,0), (nu, nv),
nu=nu, nv=nv,
usplinetype=Spline, vsplinetype=Spline);
}

int increment=1;

material[] surfacecolors = new material[0];

triple[][] points = new triple[0][0];

for (int i=1; i<51; i=i+increment)
{
string filename1 = filebasename + string(i) + "PartI.dat";
file in1=input(filename1).line().csv();
real[][] a1=in1.dimension(0,0); // 0 means up to the end of the file
a1=transpose(a1);

real[] x=a1[0];
real[] y=a1[1];
real[] z=a1[3];

triple[] currentrow = new triple[0];
int step = 1;
for (int j = 0; j < 21; j += step) {
currentrow.push((x[j],y[j],z[j]));

// Remember the color in which this patch should be drawn.
if (i < 51-increment && j < 21-step)
surfacecolors.push(
red*(1-sqrt(1-i/50)) + rgb(0,62/255,91/255)*(sqrt(1-i/50)) + opacity(1));
}
points.push(currentrow);
}

surface smoothsurface = smoothSurface(points);

draw(smoothsurface, surfacecolors);

////// CONTOUR OPENING //////

string fileImpactData = "./data/impactData.dat";
file in1=input(fileImpactData).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
a1=transpose(a1);

real[] x1=a1[0];
real[] y1=a1[1];
real[] z1=a1[3];

pen contourpen=green+1.5bp;

draw(graph(x1,y1,z1,operator--),contourpen,currentlight);

////// PLANES ///////
pen bg=gray(0.9)+opacity(0.2);
draw(surface((xmax,ymin,zmin)--(xmax,ymin,zmax)--(xmin,ymin,zmax)--(xmin,ymin,zmin)--cycle),bg);
draw(surface((xmin,ymax,zmin)--(xmin,ymax,zmax)--(xmin,ymin,zmax)--(xmin,ymin,zmin)--cycle),bg);
draw(surface((xmax,ymax,zmin)--(xmax,ymin,zmin)--(xmin,ymin,zmin)--(xmin,ymax,zmin)--cycle),bg);

////// GRID LINES ///////
pen gridpen=.2bp+gray(0.7);

grid3(XYgrid,Step=1,gridpen);
grid3(YXgrid,Step=.5,gridpen);
grid3(XZgrid,Step=1,gridpen);
grid3(ZXgrid,Step=2,gridpen);
grid3(YZgrid,Step=.5,gridpen);
grid3(ZYgrid,Step=2,gridpen);

// No-go zone
draw((xmax,1,zmin)--(xmin,1,zmin)--(xmin,1,zmax),black+1bp);

xaxis3(Label("$x_1$",MidPoint,align=Y-Z),Bounds(Both,Min),InTicks(Step=1),p=black);
yaxis3(Label("$x_2$",MidPoint,align=X-Z),Bounds(Both,Min),InTicks(Step=.5),p=black);
zaxis3(Label("$\dot x_2$",MidPoint,align=X-Y),Bounds(Both,Min),InTicks(Step=2),p=black);

Here's the result:

• I know thanks should be avoided, but still I wanted to warmly thank you. You helped me a lot. Feb 3, 2015 at 6:05
• Hello Charles. The very first one solution which combines the "rectangular" patches is very interesting and allows to deal with more general non orderer data (for rectangular finite elements for example even if the mesh is stored in an appropriate format). Now both solutions are equivalent, juste replace the data by a function.
– O.G.
Feb 4, 2015 at 8:13

Since the data are very well organized it is possible to construct a parametrized Bézier patches surface in the same spirit as the surface parametrized asymptote routine. In fact the surface described by a real function over a box(a,b) (only on a grid, f defined from the box to R) is an adaptation of a Scilab function.

For a surface described by a parametric function f over box(a,b) the process is the following (f takes value in 3D): if f depends on (t,s) variable, Dt a uniform subdivision in t with nu points, Ds a uniform subdivision in s with nv points. F the value of f on the grid, Fx, Fy, Fz the x, y, z values. Using the previous routine we construct 3 Bézier parametrized surfaces; (Fx,{1,...,nu},{1,...,nv}), (Fy,{1,...,nu},{1,...,nv}), (Fz,{1,...,nu},{1,...,nv}). At last we obtain the Bézier parametrized surface (Fx,Fy,Fz). It is a very naive method and in 2D it is well known that you can have some inappropriate behavior (see De Boor). If someone has a better parametrized choice (instead of uniform) please let me know.

// Global Asymptote definitions can be put here.
import three;
import grid3;

usepackage("mathptmx");

// One can globally override the default toolbar settings here:
// settings.toolbar=true;

import graph3;

real xmin=-2, xmax=2;
real ymin=-2, ymax=1.5;
real zmin=-2.5, zmax=2.5;

limits((xmin,ymin,zmin),(xmax,ymax,zmax));

currentprojection=perspective(camera=(1.5,2,2.5));

unitsize(3cm,3cm,2cm);

real linewidth=1.1;
real linewidthprojections=.15;

string filebasename="./data/ForASY_mnl2ippGridBranche3T1_";

bool renderPRC = false;

if(renderPRC) {
// PRC TRUE
settings.prc=true;
settings.embed=true;
}
else {
// RASTERIZE
settings.outformat="png";
settings.prc=false;
settings.render=3;
}

/////// ORBITS IN 3D, SMOOTH ////////
for (int i=1; i<=50; i=i+2)
{
string filename1 = filebasename + string(i) + "PartI.dat";
file in1=input(filename1).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
a1=transpose(a1);

real[] x1=a1[0];
real[] y1=a1[1];
real[] z1=a1[3];

pen orbitpen=.7bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50));
pen projpen=.3bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50))+opacity(0.3);

path3 thepath1 = graph(x1, y1, z1, operator..);
draw(thepath1,orbitpen,currentlight);
}

//////// SURFACE PART I ////////

triple[][] v=new triple[50][21];

int increment=1;
for (int i=1; i<52-increment; i=i+increment)
{
string filename1 = filebasename + string(i) + "PartI.dat";
file in1=input(filename1).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
a1=transpose(a1);
real[] x=a1[0];
real[] y=a1[1];
real[] z=a1[3];
int step=1;
for (int j=0; j<22-step; j=j+step)
{
v[i-1][j]=(x[j],y[j],z[j]);

}
}

int nu=50-1;
int nv=21-1;
real[] ipt=sequence(nu+1);
real[] jpt=sequence(nv+1);
real[][] fx=new real[nu+1][nv+1];
real[][] fy=new real[nu+1][nv+1];
real[][] fz=new real[nu+1][nv+1];

splinetype[] usplinetype=Spline;
splinetype[] vsplinetype=Spline;

for(int i=0; i <nu+1; ++i) {
real ui=i;
real[] fxi=fx[i];
real[] fyi=fy[i];
real[] fzi=fz[i];
for(int j=0; j < nv+1; ++j) {
pair z=(ui,j);
fxi[j]=v[i][j].x;
fyi[j]=v[i][j].y;
fzi[j]=v[i][j].z;
}
}
real[][][] sx=bispline(fx,ipt,jpt);//,usplinetype[0],vsplinetype[0]);
real[][][] sy=bispline(fy,ipt,jpt);//,usplinetype[1],vsplinetype[1]);
real[][][] sz=bispline(fz,ipt,jpt);//,usplinetype[2],vsplinetype[2]);

surface s=surface(sx.length);
s.index=new int[nu][nv];
int k=-1;
for(int i=0; i < nu; ++i) {
int[] indexi=s.index[i];
for(int j=0; j < nv; ++j)
indexi[j]=++k;
}
for(int k=0; k < sx.length; ++k) {
triple[][] Q=new triple[4][];
real[][] Px=sx[k];
real[][] Py=sy[k];
real[][] Pz=sz[k];
for(int i=0; i < 4 ; ++i) {
real[] Pxi=Px[i];
real[] Pyi=Py[i];
real[] Pzi=Pz[i];
Q[i]=new triple[] {(Pxi[0],Pyi[0],Pzi[0]),
(Pxi[1],Pyi[1],Pzi[1]),
(Pxi[2],Pyi[2],Pzi[2]),
(Pxi[3],Pyi[3],Pzi[3])};
}
s.s[k]=patch(Q);
}

draw(s,red);
////// CONTOUR OPENING //////

string fileImpactData = "./data/impactData.dat";
file in1=input(fileImpactData).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
a1=transpose(a1);

real[] x1=a1[0];
real[] y1=a1[1];
real[] z1=a1[3];

pen contourpen=green+1.5bp;

draw(graph(x1,y1,z1,operator--),contourpen,currentlight);

////// PLANES ///////
pen bg=gray(0.9)+opacity(0.2);
draw(surface((xmax,ymin,zmin)--(xmax,ymin,zmax)--(xmin,ymin,zmax)--(xmin,ymin,zmin)--cycle),bg);
draw(surface((xmin,ymax,zmin)--(xmin,ymax,zmax)--(xmin,ymin,zmax)--(xmin,ymin,zmin)--cycle),bg);
draw(surface((xmax,ymax,zmin)--(xmax,ymin,zmin)--(xmin,ymin,zmin)--(xmin,ymax,zmin)--cycle),bg);

////// GRID LINES ///////
pen gridpen=.2bp+gray(0.7);

grid3(XYgrid,Step=1,gridpen);
grid3(YXgrid,Step=.5,gridpen);
grid3(XZgrid,Step=1,gridpen);
grid3(ZXgrid,Step=2,gridpen);
grid3(YZgrid,Step=.5,gridpen);
grid3(ZYgrid,Step=2,gridpen);

// No-go zone
draw((xmax,1,zmin)--(xmin,1,zmin)--(xmin,1,zmax),black+1bp);

xaxis3(Label("$x_1$",MidPoint,align=Y-Z),Bounds(Both,Min),InTicks(Step=1),p=black);
yaxis3(Label("$x_2$",MidPoint,align=X-Z),Bounds(Both,Min),InTicks(Step=.5),p=black);
zaxis3(Label("$\dot x_2$",MidPoint,align=X-Y),Bounds(Both,Min),InTicks(Step=2),p=black);

The code is not optimized. You obtain the picture. I do not manage the color palette.

It is possible to have a very short solution (but not very readable). The "Bézier parametrized" process (50 lines) is replaced by 5 lines.

// Global Asymptote definitions can be put here.
import three;
import grid3;
import palette;

usepackage("mathptmx");

// One can globally override the default toolbar settings here:
// settings.toolbar=true;

import graph3;

real xmin=-2, xmax=2;
real ymin=-2, ymax=1.5;
real zmin=-2.5, zmax=2.5;

limits((xmin,ymin,zmin),(xmax,ymax,zmax));

currentprojection=perspective(camera=(1.5,2,2.5));

unitsize(3cm,3cm,2cm);

real linewidth=1.1;
real linewidthprojections=.15;

string filebasename="./data/ForASY_mnl2ippGridBranche3T1_";

bool renderPRC = false;

if(renderPRC) {
// PRC TRUE
settings.prc=true;
settings.embed=true;
}
else {
// RASTERIZE
settings.outformat="png";
settings.prc=false;
settings.render=3;
}

/////// ORBITS IN 3D, SMOOTH ////////
for (int i=1; i<=50; i=i+2)
{
string filename1 = filebasename + string(i) + "PartI.dat";
file in1=input(filename1).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
a1=transpose(a1);

real[] x1=a1[0];
real[] y1=a1[1];
real[] z1=a1[3];

pen orbitpen=.7bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50));
pen projpen=.3bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50))+opacity(0.3);

path3 thepath1 = graph(x1, y1, z1, operator..);
draw(thepath1,orbitpen,currentlight);
}

//////// SURFACE PART I ////////

triple[][] v=new triple[50][21];

int increment=1;
for (int i=1; i<52-increment; i=i+increment)
{
string filename1 = filebasename + string(i) + "PartI.dat";
file in1=input(filename1).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
a1=transpose(a1);
real[] x=a1[0];
real[] y=a1[1];
real[] z=a1[3];
int step=1;
for (int j=0; j<22-step; j=j+step)
{
v[i-1][j]=(x[j],y[j],z[j]);

}
}
triple f (pair t) {
return(v[round(t.x)][round(t.y)]);
}

surface ns=surface(f,(0,0),(49,20),49,20,Spline);
ns.colors(palette(ns.map(zpart),Rainbow()));
draw(ns,render(merge=true));

////// CONTOUR OPENING //////

string fileImpactData = "./data/impactData.dat";
file in1=input(fileImpactData).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
a1=transpose(a1);

real[] x1=a1[0];
real[] y1=a1[1];
real[] z1=a1[3];

pen contourpen=green+1.5bp;

draw(graph(x1,y1,z1,operator--),contourpen,currentlight);

////// PLANES ///////
pen bg=gray(0.9)+opacity(0.2);
draw(surface((xmax,ymin,zmin)--(xmax,ymin,zmax)--(xmin,ymin,zmax)--(xmin,ymin,zmin)--cycle),bg);
draw(surface((xmin,ymax,zmin)--(xmin,ymax,zmax)--(xmin,ymin,zmax)--(xmin,ymin,zmin)--cycle),bg);
draw(surface((xmax,ymax,zmin)--(xmax,ymin,zmin)--(xmin,ymin,zmin)--(xmin,ymax,zmin)--cycle),bg);

////// GRID LINES ///////
pen gridpen=.2bp+gray(0.7);

grid3(XYgrid,Step=1,gridpen);
grid3(YXgrid,Step=.5,gridpen);
grid3(XZgrid,Step=1,gridpen);
grid3(ZXgrid,Step=2,gridpen);
grid3(YZgrid,Step=.5,gridpen);
grid3(ZYgrid,Step=2,gridpen);

// No-go zone
draw((xmax,1,zmin)--(xmin,1,zmin)--(xmin,1,zmax),black+1bp);

xaxis3(Label("$x_1$",MidPoint,align=Y-Z),Bounds(Both,Min),InTicks(Step=1),p=black);
yaxis3(Label("$x_2$",MidPoint,align=X-Z),Bounds(Both,Min),InTicks(Step=.5),p=black);
zaxis3(Label("$\dot x_2$",MidPoint,align=X-Y),Bounds(Both,Min),InTicks(Step=2),p=black);

A similar picture with a color palette.

• Sorry for my late comment and thank you for your very efficient solution. I have never used palette, but I am not sure if it can easily be parametrized as a function of the file number (rather than z as in ns.map(zpart)). Anyway, I used Charles Staat's code as it came first. The computation time to activate navigation in pdf (with prc) is the same, as well as the size of the file. Again, thank you! Feb 8, 2015 at 3:34
• You're welcome. In fact my goal was only to have the smooth surface and to give some explanations. The palette option comes from the asymptote gallery and I missed it from your message.
– O.G.
Feb 8, 2015 at 14:13
• Could you provide a more specific reference for the bad behavior of this algorithm? De Boor did a whole lot of stuff. Mar 11, 2015 at 5:28
• You're right. "A pratical guide to spline", revised edition, Springer : chapter 16, Taut Splines, Periodic splines, cardinal splines, pages 277-279 example 19 : choice of parametrization is important (yes it is !). I made some tests a few years ago. It is in relation with the first solution I gave in tex.stackexchange.com/questions/231515/draw-a-smooth-surface and the spline option : a very bad result.
– O.G.
Mar 12, 2015 at 12:27

Consider this example. Does this result illustrate what you need?

As you can see in the following code, the 2-array P stores the coordinates and s is a surface constructed by P. Then s is colored by an array of pens.

import three;
import palette;

size(12cm);
currentprojection=orthographic(1,1,1.5);
currentlight=(1,0,1);

triple P00=-X-Y+0.5*Z, P03=-X+Y, P33=X+Y, P30=X-Y;

triple[][] P={
{P00,P00+(-0.5,0.5,-1),P03+(0,-0.5,1),P03},
{P00+(0.5,-0.5,-1),(-0.5,-0.5,0.5),(-0.5,0.5,-1.5),P03+(0.5,0,1)},
{P30+(-0.5,0,1),(0.5,-0.5,-1.5),(0.5,0.5,1),P33+(-0.5,0,1)},
{P30,P30+(0,0.5,1),P33+(0,-0.5,1),P33}
};

surface s=surface(patch(P));
// s.colors(palette(s.map(zpart),Rainbow()));

draw(s);
draw(sequence(new path3(int i){
return s.s[i].external();},s.s.length), bp+orange);

if(!is3D())
shipout(bbox(Fill(lightgrey)));

It looks good. However in another example it warns

Here we determine the colors of vertexes (vertex shading). Since the PRC output format does not support vertex shading of Bezier surfaces, PRC patches are shaded with the mean of the four vertex colors.

So there is nothing to do with asymptote as long as you want PRC.

Warning: on my computer, OpenGL has a hard time to handle thousands of bezier surfaces. So please test a lot before you present it.

• In the given example, P contains not only the points (4 triples), but three more rows to define the tangents (3*4 triples). The given result is for one quadrangle; if I had to calculate the 12 spline parameters so that the surface is smooth, that'd be extremely tedious... I'd like Asymptote to calculate the parameters of the smooth surface for me. I know it is not an easy problem (with several possible solutions), but I'm surprised I spent so many hours on this and still see no end to it--not blaming asymptote but my skills. I did not manage to make Charles code work for open surfaces yet. Feb 1, 2015 at 18:59
• @anderstood Sorry I did not study the case carefully. But I saw, in the second link I provide, an construction surface s=surface(f,(0,0),(pi,2pi),100,Spline); Doesn't the keyword Spline realize your dream? Feb 2, 2015 at 1:46
• no worries, that was not obvious, it took me quite some time to understand what P was. The construction you indicate is similar to that used by Charles Staat (see comments below OP). I tried with that too, but in vain. I think the main reason why I fail is that I don't have an explicit function f, but a list of points, so the calculation of the derivatives (to make a smooth surface) is not straightforward. Feb 2, 2015 at 1:51
• @anderstood I had a small experiment that if you send 100 to assign the size of samples, the function f is called 101*101 times. That means you can patch your data to make it looks like a function, right? Feb 2, 2015 at 2:04
• Also in asymptote manual there is a declaration surface surface(triple[][] f, bool[][] cond={});. Which is followed by // return the surface described by a matrix f. Feb 2, 2015 at 2:14