Cube C = new Cube(50,45,30,35,15,20,10,5,200);
//Cube C = new Cube(20,20,20,20,0,0,0,0,200);
boolean pick=false;
boolean picked=false; 
float threshold=0;
float m=0.25;
boolean showRefined=false;

class Cube {
   float s=1; // size
   int pi=0, pj=0, pk=0;  // id of picked corner
   float[][][] r = new float [2][2][2];  // values at corners
   pt V[][][] = new pt [2][2][2];    // corner points

  Cube(float ps) {s=ps; makeBalls();}
 
  Cube(float a, float b, float c, float d, float e, float f, float g, float h,  float ps) {s=ps; makeBalls();
    r[0][0][0]=a; r[0][0][1]=b; r[0][1][0]=c; r[0][1][1]=d; r[1][0][0]=e; r[1][0][1]=f; r[1][1][0]=g; r[1][1][1]=h;  }
 
  Cube(float a, float b, float c, float d, float e, float f, float g, float h) {makeBalls();
    r[0][0][0]=a; r[0][0][1]=b; r[0][1][0]=c; r[0][1][1]=d; r[1][0][0]=e; r[1][0][1]=f; r[1][1][0]=g; r[1][1][1]=h;  }

  void makeBalls() { for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  for (int k=0; k<2; k++) {V[i][j][k]=new pt(s*(2*i-1),s*(2*j-1),s*(2*k-1)); r[i][j][k]=s/25;}; 
            threshold=average();  }

  void changeRadius() {r[pi][pj][pk]+=(mouseX-pmouseX)/10.; if(r[pi][pj][pk]<0) r[pi][pj][pk]=0;  }
 
  float average() {float a=0; for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  for (int k=0; k<2; k++) a+=r[i][j][k]; return a/8;}
 
  void showBalls() { noStroke();
    for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  for (int k=0; k<2; k++) {
      if(r[i][j][k]<threshold) fill(blue); else fill(green); 
 //     if((i==pi)&&(j==pj)&&(k==pk)) fill(orange); 
    V[i][j][k].show(4*r[i][j][k]*m);};    }

  boolean pickBall(pt M) {boolean picked = false; float dm=10000.0; 
    for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  for (int k=0; k<2; k++) {
     float dd = M.disTo(V[i][j][k]); 
     if((dd<r[i][j][k]*10*m)&&(dd<dm)) {picked=true; dm=M.disTo(V[i][j][k]); pi=i; pj=j; pk=k;};  }; return picked;}
 
  void showEdges() {fill(green); noStroke(); 
    for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  cone(V[i][j][0], r[i][j][0]*m, V[i][j][1], r[i][j][1]*m, 12 );
    for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  cone(V[i][0][j], r[i][0][j]*m, V[i][1][j], r[i][1][j]*m, 12 );
    for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  cone(V[0][i][j], r[0][i][j]*m, V[1][i][j], r[1][i][j]*m, 12 );
    }
 
  void showThreshold() {fill(blue); noStroke(); 
    for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  cone(V[i][j][0], threshold*m, V[i][j][1], threshold*m, 12 );
    for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  cone(V[i][0][j], threshold*m, V[i][1][j], threshold*m, 12 );
    for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  cone(V[0][i][j], threshold*m, V[1][i][j], threshold*m, 12 );
    }

  void showVertices() {fill(red); noStroke(); 
    for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  interpolate(V[i][j][0], r[i][j][0], V[i][j][1], r[i][j][1]);
    for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  interpolate(V[i][0][j], r[i][0][j], V[i][1][j], r[i][1][j]);
    for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  interpolate(V[0][i][j], r[0][i][j], V[1][i][j], r[1][i][j]);
    }
 
  void makeIsoVertices(Mesh M) {fill(red); noStroke(); 
    for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  addIsoVertex(V[i][j][0], r[i][j][0], V[i][j][1], r[i][j][1], M);
    for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  addIsoVertex(V[i][0][j], r[i][0][j], V[i][1][j], r[i][1][j], M);
    for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  addIsoVertex(V[0][i][j], r[0][i][j], V[1][i][j], r[1][i][j], M);
  }

  void makeMesh(Mesh M) { 
     M.reset();
     addCornerVertices(M);
     addCubeTriangles(M);
     M.init();
     };
   
  void addCornerVertices(Mesh M) {for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  for (int k=0; k<2; k++) M.addVertex(V[i][j][k],r[i][j][k]); }
  void updateCornerVertices(Mesh M) {int v=0; for (int i=0; i<2; i++)  for (int j=0; j<2; j++)  for (int k=0; k<2; k++) M.rad[v++]=r[i][j][k]; }

  void addCubeTriangles(Mesh M) {
    M.addTriangle(0,1,3); M.addTriangle(0,3,2); 
    M.addTriangle(0,2,4); M.addTriangle(4,2,6); 
    M.addTriangle(2,3,6); M.addTriangle(3,7,6); 
    M.addTriangle(0,4,1); M.addTriangle(1,4,5); 
    M.addTriangle(1,5,3); M.addTriangle(3,5,7); 
    M.addTriangle(4,7,5); M.addTriangle(4,6,7); 
    };


  float value (pt P) {float x=(P.x/s+1)/2.; float y=(P.y/s+1)/2.; float z=(P.z/s+1)/2.;  
    return lin(lin(lin(r[0][0][0],x,r[1][0][0]),y,lin(r[0][1][0],x,r[1][1][0])),
             z,lin(lin(r[0][0][1],x,r[1][0][1]),y,lin(r[0][1][1],x,r[1][1][1]))) ; }
  
  float lin (float a, float x, float b) { return a+x*(b-a);}
 
  }  // end Cube class

void cone(pt A, float a, pt B, float b, int s) {
  vec T = A.vecTo(B); float d=T.norm(); T.makeUnit();
  vec U = new vec (1,0,0); vec N = cross(U,T);
  if (N.norm()<0.1) {U.setTo(0,1,0);  N.setTo(cross(U,T)); };
  N.makeUnit(); U.setTo(cross(N,T));
  float dd = A.disTo(B);
  float x=(a-b)/dd;
  float y=sqrt((dd-a+b)*(dd+a-b))/dd;
  pt PP = A.make();   pt QQ = B.make();
  PP.addScaledVec(a*x,T);   QQ.addScaledVec(b*x,T);
  pt P = new pt(); pt Q = new pt();
  beginShape(QUAD_STRIP);
  for (float t=0; t<TWO_PI; t+=TWO_PI/s) {
    P.setTo(PP);   Q.setTo(QQ);
    P.addScaledVec(a*y*cos(t),N);   Q.addScaledVec(b*y*cos(t),N); 
    P.addScaledVec(a*y*sin(t),U);   Q.addScaledVec(b*y*sin(t),U); 
   // P.showLineTo(Q);
   P.vert(); Q.vert();
    };
  endShape();
  }

void interpolate(pt A, float a, pt B, float b) {
  if ((min(a,b)<threshold)&&(threshold<max(a,b))) {
  float s = (threshold-a)/(b-a);
  pt P = A.make();  P.moveTowards(s,B);  P.show(threshold*4*m); 
  }; 
}

void addIsoVertex(pt A, float a, pt B, float b, Mesh M) {
  if ((min(a,b)<threshold)&&(threshold<max(a,b))) {
    float s = (threshold-a)/(b-a);
    pt P = A.make();  P.moveTowards(s,B);  M.addVertex(P,r);  
    };
  }