boolean showTerrain=false;
int n=64;                                   // size of grid. Must be >2!
p3[][] G = new p3 [n][n];                   // array of vertices
v3[][] L = new v3 [n][n];                 // laplace vectors for smoothing
float w,h,ww,hh;                             // width, height of cell in absolute and normalized units
//float ax=342, az=1557;                          // camera angles 
float ax=800, az=0;                          // camera angles 
boolean showImage=false;
float zmax=400, zmin=0;

void initTerrain() {
  ww=1.0/(n-1); hh=1.0/(n-1);                                      // set intial width and height of a cell
  w=height*ww; h=height*hh;                                         // set intial width and height of a cell in normalized [0,1]x[0,1]
  resetVertices();
  }


void resetVertices() { for (int i=0; i<n; i++) for (int j=0; j<n; j++) G[i][j]=new p3(i*w,j*h,0);  } 

void showTerrain() {
  lights();
  translate(width/2,height/2,-150);
  rotateX(PI/2.-ax*PI/height/2.);
  rotateZ(-az*PI/width/4.);
  translate(-width/2,-height/2,0);  
  if  (mousePressed) if (!keyPressed)  {ax+=mouseY-pmouseY;  az+=mouseX-pmouseX; println("ax="+ax+", az="+az);}  ;
  if (showImage) scale(1,1,0); 
  paintTerrain();
  }

void paintTerrain() {         // renders the terrain with color as a function of height
   noStroke(); fill(200,200,200);        
   for (int i=0; i<n-1; i++) {
     beginShape(QUAD_STRIP); 
     for (int j=0; j<n; j++) { 
           float zi=(G[i][j].z-zmin)/(zmax-zmin), zip=(G[i+1][j].z-zmin)/(zmax-zmin); 
  //         fill(zi*256,zi*256,zi*256);
           G[i][j].vert();  G[i+1][j].vert(); 
         };
     endShape();
     };
   }

void distanceTransform() {
  zmax=0;
  for (int i=0; i<n; i++) for (int j=0; j<n; j++) {
    pt S= new pt(i*h,j*h);
    if (C.contains(S)) {
       closest=C.project(S); 
       G[i][j].z=S.disTo(closest); 
       zmax=max(zmax,G[i][j].z);
       };
    };
  }

void ballTransform() {
  zmax=400;
  for (int i=0; i<n; i++) for (int j=0; j<n; j++) {
    pt S= new pt(i*h,j*h);
    if (C.contains(S)) {
       closest=C.project(S); 
       scanSphere(S,S.disTo(closest));
      };
    };
  }
  
void stabilityTransform() {
  zmax=400;
  for (int i=0; i<n; i++) for (int j=0; j<n; j++) {
    pt S= new pt(i*h,j*h);
    if (C.contains(S)) {
       closest=C.project(S); 
       scanDisk(S,S.disTo(closest));
      };
    };
  }

void scanSphere(pt O, float r) {
   float cx=O.x/h, cy=O.y/h, cr=r/h;
   for (int i=max(0,int(cx-cr)); i<min(n,cx+cr); i++) for (int j=max(0,int(cy-cr)); j<min(n,cy+cr); j++) 
   if(sq(cx-i)+sq(cy-j)<sq(cr)) 
  {
    float z=sqrt(sq(cr)-sq(cx-i)-sq(cy-j))*h;    
    if(z>G[i][j].z) G[i][j].z=z;
    zmax=max(zmax,G[i][j].z);
   };
  }

void scanDisk(pt O, float r) {
   float cx=O.x/h, cy=O.y/h, cr=r/h;
   for (int i=max(0,int(cx-cr)); i<min(n,cx+cr); i++) for (int j=max(0,int(cy-cr)); j<min(n,cy+cr); j++) 
   if(sq(cx-i)+sq(cy-j)<sq(cr)) 
  {
    float z=cr*h;    
    if(z>G[i][j].z) G[i][j].z=z;
    zmax=max(zmax,G[i][j].z);
   };
  }


void XYcubicFilter() {
   for(int k=0; k<20; k++) {
    for (int i=2; i<n-2; i++) for (int j=2; j<n-2; j++) L[i][j]=average(cubicFitX(i,j),cubicFitY(i,j));
    for (int i=2; i<n-2; i++) for (int j=2; j<n-2; j++) {pt S= new pt(i*w,j*h); if (C.contains(S)) G[i][j].addScaledVec(0.5,L[i][j]); };
   };
  }

v3 cubicFitX(int i, int j) {return(cubic(G[i-2][j], G[i-1][j],G[i][j],G[i+1][j],G[i+2][j] ));}
v3 cubicFitY(int i, int j) {return(cubic(G[i][j-2], G[i][j-1],G[i][j],G[i][j+1],G[i][j+2] ));}
v3 cubic (p3 A, p3 B, p3 C, p3 D, p3 E) 
    {return(new v3 ( (-A.x+4*B.x-6*C.x+4*D.x-E.x)/6, (-A.y+4*B.y-6*C.y+4*D.y-E.y)/6, (-A.z+4*B.z-6*C.z+4*D.z-E.z)/6  ));}