// smoothing
// Written by Jarek Rossignac August 2006
float M=10;                    // values range from 0 to M
int N=30;                      // number of samples
float V[] = new float[N];     // values at samples
boolean C[] = new boolean[N];     // if constrained

float V1[] = new float[N];     // values at samples
float L1[] = new float[N];     // laplace adjustments

float V2[] = new float[N];     // values at samples
float L2[] = new float[N];     // laplace adjustments
float P2[] = new float[N]; 

float V3[] = new float[N];     // values at samples
float L3[] = new float[N];     // laplace adjustments
float R3[] = new float[N]; 
float P3[] = new float[N]; 
float Q3[] = new float[N]; 

int n1=0, n2=0, n3=0;          // Iteration counters
float t=0, s=0;               //  Error measures for C3

boolean showP=false,  showQ=false, showR=false, showL=false; 
char stp=' '; // why C1 stopped

void setup() { 
  size(800, 800, P3D);  
  PFont font = loadFont("Courier-48.vlw"); textFont(font, 14);  // font for writing labels on screen
  for (int i=0; i<N; i++) V[i]=V1[i]=random(0,M/3)+M/2;
  for (int i=0; i<N; i++) C[i]=false; 
  V[0]=M*0.85; V[1]=M*0.85;   V[N-2]=M*0.85;  V[N-1]=M*0.85; 
  int m=int(0.7*N); V[m]=M*0.30; C[m]=true; 
  reset();          // reserts V1, V2, V3 to V and reserts counters
  for (int i=0; i<2*N; i++) fairC3(); setFrom3();
  } 
  
void reset() {
    for (int i=0; i<N; i++) { V1[i]=V[i]; V2[i]=V[i]; V3[i]=V[i];  };
    C[0]=C[1]=true;  C[N-1]=C[N-2]=true; 
    lorInit();         // init Lorenzo smoothing
    initFairC3();
    n1=0; n2=0; n3=0;   }



void setFrom3() {for (int i=0; i<N; i++) V[i]=V3[i]; reset();}
 
void draw() { 
   background(255);   
   stroke(0,0,0); showCurve(); 
   stroke(0,0,0); showConstraints();
   String ss="n1="+n1+""+stp+", n2="+n2+", n3="+n3+""+stp+", s="+s+", t="+t; fill(0); text(ss,10,15);   
   stroke(0,250,0); showCurve2();
   stroke(0,0,250); showCurve3();
   stroke(250,0,0); showCurve1();
   if (showP) showP3(); if (showQ) showQ3(); if (showR) showR3(); if (showL) showL3();
   }

void showCurve() { beginShape(LINE_STRIP);  for (int i=0; i<N; i++) vertex(i*width/N,V[i]/M*height); endShape();   }
void showCurve1() { beginShape(LINE_STRIP);  for (int i=0; i<N; i++) vertex(i*width/N,V1[i]/M*height); endShape();   }
void showCurve2() { beginShape(LINE_STRIP);  for (int i=0; i<N; i++) vertex(i*width/N,V2[i]/M*height); endShape();   }
void showCurve3() { beginShape(LINE_STRIP);  for (int i=0; i<N; i++) vertex(i*width/N,V3[i]/M*height); endShape();   }
void showConstraints() {for (int i=0; i<N; i++) if (C[i]) ellipse(i*width/N,V[i]/M*height,5,5);   }

void keyPressed() {  
   if (key=='1') {fairC1();};  if (key=='!') {for (int i=0; i<2*N; i++) fairC1(); }; 
   if (key=='2') {fairC2();};  if (key=='@') {for (int i=0; i<N; i++) fairC2(); }; 
   if (key=='3') {fairC3();};  if (key=='#') {for (int i=0; i<2*N; i++) fairC3(); };  
   if (key=='p') showP=!showP;   if (key=='q') showQ=!showQ;   if (key=='r') showR=!showR; if (key=='l') showL=!showL;
   if (key=='i') setFrom3();  
   };

void mousePressed() {float x = mouseX; x=x/width*N; float y = float(mouseY)/height*M; 
int i=int(x+0.5); V[i]=y; C[i]=!C[i];  reset(); }

void fairC1 () {if(lorStep()) n1++; }   // Lorenzo

void fairC1b () {aim1(); tuck1(0.5); aim1(); tuck1(-0.5); snap1(); n1++;}   // back&forth laplace V+=sL(V); V-=sL(V);
void aim1 () {for (int i=1; i<N-1; i++) L1[i]=(V1[i-1]+V1[i+1])/2-V1[i]; }
void tuck1 (float s) {for (int i=0; i<N-1; i++) V1[i]+=s*L1[i]; }
void snap1 () {for (int i=0; i<N; i++) if (C[i]) V1[i]=V[i]; }
  
void fairC2() {aim2(); fix2(); tuck2(0.5);  snap2(); n2++; }              // biaplace V+=ssL(L(V));
void aim2 () {for (int i=1; i<N-1; i++) L2[i]=(V2[i-1]+V2[i+1])/2-V2[i]; }  // L(V)
void fix2 () {for (int i=1; i<N-1; i++) P2[i]=(L2[i-1]+L2[i+1])/2-L2[i]; }  // P=L(L(V))
void tuck2 (float s) {for (int i=0; i<N-1; i++) V2[i]-=s*s*P2[i]; }         // V+=P
void snap2 () {for (int i=0; i<N; i++) if (C[i]) V2[i]=V[i]; }

void initFairC3() {
    for (int i=2; i<N-2; i++) {P3[i]=0; Q3[i]=0;};
    for (int i=2; i<N-2; i++) L3[i]=(V3[i-1]+V3[i+1])/2-V3[i]; 
                 // for (int i=2; i<N-2; i++) R3[i]=- ((L3[i-1]+L3[i+1])/2-L3[i]); 
    for (int i=2; i<N-2; i++) R3[i]= - ( V3[i-2]/2 - 2*V3[i-1] + 3*V3[i] - 2*V3[i+1] + V3[i+2]/2); 
  }
float CP(int i) {float v = P3[i]; if (C[i]) v=0; return(v);}

void fairC3 () {if(jrStep()) n3++; }   // Lorenzo

boolean jrStep() {
    s=0; for (int i=2; i<N-2; i++) if (!C[i]) s+=sq(R3[i]); 
    if (s<0.000001) {println("s="+s+" too small"); stp='s'; return(false);   };        // early break (residues are zero)
    for (int i=2; i<N-2; i++) P3[i]+=R3[i]/s; 
    for (int i=2; i<N-2; i++) if (C[i]) Q3[i]=0; else Q3[i]=3*CP(i) -2*CP(i-1) -2*CP(i+1) +CP(i-2)/2 +CP(i+2)/2; 
    t = 0; for (int i=2; i<N-2; i++) if (!C[i]) t+=Q3[i]*P3[i];  
    if (s*t <= 6/1000) {println("s*t too small");  stp='t'; return(false);   };
    for (int i=2; i<N-2; i++) if (!C[i]) R3[i]-=Q3[i]/t;
    for (int i=2; i<N-2; i++) if (!C[i]) V3[i]+=P3[i]/t;
    stp=' '; return(true);
    }

void showL3() {stroke(100,100,100); for (int i=2; i<N-2; i++)  if (!C[i]) line(i*width/N , V3[i]/M*height , i*width/N , (V3[i]+L3[i])/M*height);  }
void showP3() {stroke(0,200,200); for (int i=2; i<N-2; i++)  if (!C[i])  line(i*width/N , V3[i]/M*height , i*width/N , (V3[i]+P3[i]/t)/M*height);  }
void showQ3() {stroke(200,0,200);for (int i=2; i<N-2; i++)   if (!C[i]) line(i*width/N , V3[i]/M*height , i*width/N , (V3[i]+Q3[i])/M*height);  }
void showR3() {stroke(200,200,0); for (int i=2; i<N-2; i++)  if (!C[i])  line(i*width/N , V3[i]/M*height , i*width/N , (V3[i]+R3[i])/M*height);  }