Asymptote: Drawing contours of spline surface

Question

I have a surface specified by samples of a function on a square grid (in particular, there is no formula for the function), and I want to draw a surface plot along with some contour lines.

As you can see in the image below, just calculating contour lines with the contour package doesn't work because the contour calculation doesn't account for the surface plot smoothing. (Actually, it doesn't seem to work even without smoothing, i.e., dropping Spline and using -- as the join operator for the contour lines).

import graph3;
import contour; 
import palette;
import x11colors; 
size3(200,IgnoreAspect);
settings.outformat="pdf";
settings.render=8;
currentprojection=perspective(20,8,2); 

// A function specified by its values on a 10 x 10 grid
real[][] z =
  {{0,0,0,0,0,0,0,0,0,0},
   {0,0.459792676157358,0.136148250882153,-0.163615908155485,-0.0259463996973455,-0.0840769425808095,0.62126438834251,0.518640203087028,1.03569619721556,0},
   {0,-0.604162941996047,-0.424941391270887,0.259335710132218,0.650164067086578,0.127505396459818,0.302247026415486,0.254979306532084,0.252221600145993,0},
   {0,0.0304849331645697,-0.524980730553991,0.854096681304495,0.480467910001354,0.258129875403682,0.692548824372049,-0.252410071666853,-0.167082932245632,0},
   {0,0.987911902753579,0.382866021557655,0.864790291964691,0.607114025159436,0.55793941674487,0.441528479985084,0.767320637646196,0.467364706957709,0},
   {0,0.115108721902316,0.516912837384094,1.01635622832271,0.660704188001264,0.382332982489028,-0.307456859893877,-0.382304124975139,1.06361557239264,0},
   {0,0.288185323850592,-0.117203025712584,0.803252469880009,0.714506160677705,0.151669905008431,-0.342341661772393,-0.0261977699813443,0.473208218671247,0},
   {0,0.726871540197522,0.0683933588314665,0.481762795614408,0.969192479689294,-0.40613837131098,0.142722121638684,0.635132993615454,0.237704962254851,0},
   {0,-0.769077284482047,-0.755640443336741,0.892723307000915,0.683470361768723,0.470367719595602,0.967570149005728,0.249827942333624,0.478612319838185,0},
   {0,0,0,0,0,0,0,0,0,0}};
// Dimensions of the array
int m = 10, n = 10; 

triple f(pair t) {
  int i=round(t.x);
  int j=round(t.y);
  return (i,j,z[i][j]);
}

// color and draw surface
surface s=surface(f,(0,0),(m-1,n-1),m-1,n-1,Spline);
s.colors(palette(s.map(zpart),Rainbow()));
draw(s); 

// draw a contour line with height 0.25 
real t = 0.25; 
guide[][] contours = contour(z,(0,0),(9,9),new real[]{t},join=operator ..);
guide3[][] liftedcontours = lift(new real (real x, real y) {return t;}, contours);
draw(liftedcontours[0][1],linewidth(2.0));

Answer 1

z is only discrete and there is no reason that the (smooth) contours lines of z are the contours lines of the smooth surface. You have to construct the function corresponding to the surface, it is possible through s.point and then you can use contour. Please find a possible code

import graph3;
import contour; 
import palette;
import x11colors; 
size3(200,IgnoreAspect);
//settings.outformat="pdf";
settings.render=8;
currentprojection=perspective(20,8,2); 

// A function specified by its values on a 10 x 10 grid
real[][] z =
  {{0,0,0,0,0,0,0,0,0,0},
   {0,0.459792676157358,0.136148250882153,-0.163615908155485,-0.0259463996973455,-0.0840769425808095,0.62126438834251,0.518640203087028,1.03569619721556,0},
   {0,-0.604162941996047,-0.424941391270887,0.259335710132218,0.650164067086578,0.127505396459818,0.302247026415486,0.254979306532084,0.252221600145993,0},
   {0,0.0304849331645697,-0.524980730553991,0.854096681304495,0.480467910001354,0.258129875403682,0.692548824372049,-0.252410071666853,-0.167082932245632,0},
   {0,0.987911902753579,0.382866021557655,0.864790291964691,0.607114025159436,0.55793941674487,0.441528479985084,0.767320637646196,0.467364706957709,0},
   {0,0.115108721902316,0.516912837384094,1.01635622832271,0.660704188001264,0.382332982489028,-0.307456859893877,-0.382304124975139,1.06361557239264,0},
   {0,0.288185323850592,-0.117203025712584,0.803252469880009,0.714506160677705,0.151669905008431,-0.342341661772393,-0.0261977699813443,0.473208218671247,0},
   {0,0.726871540197522,0.0683933588314665,0.481762795614408,0.969192479689294,-0.40613837131098,0.142722121638684,0.635132993615454,0.237704962254851,0},
   {0,-0.769077284482047,-0.755640443336741,0.892723307000915,0.683470361768723,0.470367719595602,0.967570149005728,0.249827942333624,0.478612319838185,0},
   {0,0,0,0,0,0,0,0,0,0}};
// Dimensions of the array
int m = 10, n = 10; 

triple f(pair t) {
  int i=round(t.x);
  int j=round(t.y);
  return (i,j,z[i][j]);
}

// color and draw surface
surface s=surface(f,(0,0),(m-1,n-1),m-1,n-1,Spline);
s.colors(palette(s.map(zpart),Rainbow()));
draw(s);
real g(real a, real b)
{
return s.point(a,b).z;
}
// draw a contour line with height 0.25 
real t = 0.25; 
guide[][] contours = contour(g,(0,0),(8.9,8.9),new real[]{t},50, join=operator ..);
//guide[][] contours = contour(z,(0,0),(9,9),new real[]{t},join=operator ..);
guide3[][] liftedcontours = lift(new real (real x, real y) {return t;}, contours);
draw(liftedcontours[0][1],linewidth(2.0));

And the picture

Notice that in contour(g,(0,0),(8.9,8.9),new real[]{t},50, join=operator ..); if 8.9 is replaced by 9 then an index error occurs. 50 (nx=50) is to reduce the accuracy, it is sufficient. A lower value produces also an index error...

O.G.