// Expand small image into large Mesh2 height field 
// using a pigment function to evaluate pixels values
// and splines to fill the mesh2 inbetween pixels values.
// Saves Smooth Mesh2 
// mesh2 centered at <0,0,0> and faces like an heightfeild 
// Auth: Leroy Whetstone 
// Email whetstoneleroy@gmail.com
// screen 600*600

#version 3.6;
 global_settings{assumed_gamma 1.0 max_trace_level 50} 

#include "colors.inc"

//============= Changes ==============================================================================

 //------- Display changes---------
 background{Blue*.5}

 #declare ShowOld=.2;// 0= no show; >0 radius; show old pixel locations with spheres

 #declare ShowPig =0;// int 0,1,2 
                     // 0=use gradient y pigment 
                     // 1=use gradient y pigment scaled with Height
                     // 2=use the Image pigment 
//----------- Change for Mesh2 -----------------------------------------------------------------------

 #declare Image=pigment{image_map { bmp "Small4.bmp"} } //hf map
          

 #declare Spline_Type=1;//spline types 0-linear 1-natural 
 #declare Cuts=8;  // int >=1 divisions between pixels; each cut expands mesh by 1 pov unit
 #declare Height=6; // float !=0 scale height of Mesh2

 #declare Normalize=1;// y/n using splines may produce height values less than Zero
                      // or greater than One before scaling with 'Height' 
                      // 'yes' here will bound Heights between 0-1 before scaling


 //--------Save to file ------------------
 // comment out to not save
 // #declare SavMesh= "LandMesh4.inc";  // save the mesh2 object


//============ end of changes =============================================================================

 #declare Isize=max_extent(Image);#debug concat("\n\nImage size= <",vstr(2,Isize,",",0,0),">\n")

 #declare Pfnc = function {pigment {Image scale<Isize.x,Isize.y,1>} }
 #declare Qfnc = function {Pfnc(x,y,z).gray}

#declare A=(Isize.x-1)*Cuts; #declare aL=A+1;
#declare B=(Isize.y-1)*Cuts; #declare bL=B+1;

#declare Jtab=array[aL][bL];

 #debug concat("\nMesh Size= <"str(A,0,0)","str(Height,0,0)","str(B,0,0)">\n\n")
 #debug concat("Mesh  Points  ",str(aL*bL,0,0),"\n") 
 #debug concat("Mesh Trangles ",str(2*A*B,0,0),"\n\n") 
 

#macro Hloc(J) // find height 
 #local TG=<J.x+.5,-J.y-.5,0>; 
 Qfnc(TG.x,TG.y,0);
#end

//---------- First Fill: with pixels values----------
#declare Cnt = 0;
#while (Cnt < aL)
 #declare Cnt2 = 0;
 #while (Cnt2 < bL)
  #declare Tp=<Cnt/Cuts,Cnt2/Cuts,0>;
  #declare Jtab[Cnt][Cnt2]=Hloc(Tp);// get height from Image & fill array
 #declare Cnt2 = Cnt2+Cuts; 
 #end
#declare Cnt = Cnt+Cuts; 
#end

//the rest uses splines to fill array

#if (Cuts>1)
 //----------- line fill 1x -------------
 #declare Cnt2=0;
 // make spline
 #while (Cnt2<bL)
 #declare SplA=spline{#if (Spline_Type=0) linear_spline 
                      #else natural_spline
                      #end
   #declare Cnt = 0;
   #while (Cnt < aL)
      #ifdef(Jtab[Cnt][Cnt2]) Cnt/A,<0,Jtab[Cnt][Cnt2]>  #end
   #declare Cnt = Cnt+1; 
   #end
  }
// use spline
  #declare Cnt = 1;
  #while (Cnt < aL)
   #ifndef(Jtab[Cnt][Cnt2]) 
    #declare Vt=SplA(Cnt/A); 
    #declare Jtab[Cnt][Cnt2]=Vt.y; 
   #end
  #declare Cnt = Cnt+1; 
  #end
 #declare Cnt2 = Cnt2+Cuts; 
 #end

 //----------- line fill 2y -------------
 #declare Cnt=0;
// make spline
 #while (Cnt<aL)
  #declare SplA=spline{#if (Spline_Type=0) linear_spline 
                       #else natural_spline
                       #end
   #declare Cnt2 = 0;
   #while (Cnt2 < bL)
       #ifdef(Jtab[Cnt][Cnt2]) 
       Cnt2/B,<0,Jtab[Cnt][Cnt2]> 
       #end
   #declare Cnt2 = Cnt2+1; 
   #end
  }
//todo: insert average deform

//use spline
  #declare Cnt2 = 1;
  #while (Cnt2 < bL)
   #ifndef(Jtab[Cnt][Cnt2]) 
     #declare Vt=SplA(Cnt2/B); 
     #declare Jtab[Cnt][Cnt2]=Vt.y; 
   #end
  #declare Cnt2 = Cnt2+1; 
  #end
 #declare Cnt = Cnt+1; 
 #end

#end

//-------------- Normalize ---------------------------------------------
#if(Normalize)
 #declare Lowp=999;
 #declare Hip=-999;
 #declare Cnt = 0;
 #while (Cnt < aL)
  #declare Cnt2 = 0;
  #while (Cnt2 < bL)
   #declare Ip=Jtab[Cnt][Cnt2];        //read array
   #if(Ip<Lowp) #declare Lowp=Ip; #end //get low point
   #if(Ip>Hip) #declare Hip=Ip; #end   //get high point
  #declare Cnt2 = Cnt2+1; 
  #end
 #declare Cnt = Cnt+1; 
 #end
  #declare Dc=1/(Hip-Lowp);
   #if(Dc=1)
    #debug"  Normalize Not needed\n"
    #debug"Low Point =0 High point =1\n\n"
   #else
    #debug"Normalize ON\n"
    #debug concat(" Low point = "str(Lowp,3,3)"\n")
    #debug concat(" Hi  point  = "str(Hip,3,3)"\n")
    #debug concat(" adjustment = "str(Dc,3,3)"\n\n")

  #declare Cnt = 0;
  #while (Cnt < aL)
   #declare Cnt2 = 0;
   #while (Cnt2 < bL)
    #declare Jtab[Cnt][Cnt2]=(Jtab[Cnt][Cnt2]-Lowp)*Dc;
   #declare Cnt2 = Cnt2+1; 
   #end
  #declare Cnt = Cnt+1; 
  #end
 #end//no Norm
#end//Norm

//---------------------------Mesh2-------------------------------------------------

#declare Pnts=array[aL*bL];
#declare Facs=array[2*A*B];
#declare Norm=array[aL*bL];

//points
   #declare N = 0;
   #declare Cnt2 = 0;
   #while (Cnt2 < bL)
    #declare Cnt = 0;
     #while (Cnt < aL)
      #declare Pnts[N]=<Cnt-(aL/2),Jtab[Cnt][Cnt2]*Height,(bL/2)-Cnt2>;
     #declare N=N+1;
     #declare Cnt = Cnt+1; 
     #end
    #declare Cnt2 = Cnt2+1; 
    #end
//faces
   #declare N = 0;
    #declare Z = 0;
     #while (Z < B)
    #declare X = 0;
      #while (X < A)
       #declare Facs[N]=<Z*(aL)+X,Z*(aL)+X+1,(Z+1)*(aL)+X+1>;
       #declare N=N+1; 
       #declare Facs[N]=<Z*(aL)+X,(Z+1)*(aL)+X+1,(Z+1)*(aL)+X>;
       #declare N=N+1; 
     #declare X = X+1; 
     #end
   #declare Z = Z+1; 
   #end
//norms
 #declare Cnt = 0;
 #while (Cnt < aL*bL)
         #declare Norm[Cnt]=<0,0,0>;
 #declare Cnt = Cnt+1; 
 #end

 #declare Cnt = 0;
 #while (Cnt < 2*A*B)
     #declare Nbase=Facs[Cnt];
     //#debug concat(" Tri < " vstr(3,Nbase,",",3,3)," >\n")
    #declare Face=vcross(Pnts[Nbase.x]-Pnts[Nbase.z],Pnts[Nbase.x]-Pnts[Nbase.y]);
    //#debug concat("   Face= < " vstr(3,Face,",",3,3)," >\n")
    #declare Norm[Nbase.x]=Norm[Nbase.x]+Face;
    #declare Norm[Nbase.y]=Norm[Nbase.y]+Face;
    #declare Norm[Nbase.z]=Norm[Nbase.z]+Face;

 #declare Cnt = Cnt+1; 
 #end
 #while (Cnt <aL*bL)
    #if((Norm[Cnt].x+Norm[Cnt].y+Norm[Cnt].z)!=0)  
      #declare Norm[Cnt]=vnormalize(Norm[Cnt]);
    #else #debug concat("Normal vector ZERO at point " str(Cnt,3,3),"\n") 
    #end
 //#debug concat("   Cnt= < " str(Cnt,3,3),"\n")
 #declare Cnt = Cnt+1; 
 #end


//--------- mesh2--------------
mesh2{
  vertex_vectors{aL*bL 
    #declare Cnt=0;
    #while(Cnt<aL*bL)
     Pnts[Cnt] 
     #declare Cnt=Cnt+1;
    #end
    }
  normal_vectors {aL*bL, 
    #declare Cnt=0;
    #while(Cnt<aL*bL)
     Norm[Cnt] 
     #declare Cnt=Cnt+1;
    #end
    } 

  face_indices{2*A*B
    #declare Cnt=0;
    #while(Cnt<2*A*B)
     Facs[Cnt] 
     #declare Cnt=Cnt+1;
    #end
}
 
// pigments

 #if(ShowPig=2) pigment{Image translate -.5 scale <A+5,B+5,1> rotate x*90}
 #else
  pigment{gradient y 
  #if(ShowPig=1) scale Height #end
  }
 #end
 }//mesh2

//show old pixel loc
#if(ShowOld)
    #declare Cnt=0;
    #while(Cnt<aL)
     #declare Cnt2=0;
     #while(Cnt2<bL)
     sphere{<Cnt-(aL/2),Jtab[Cnt][Cnt2]*Height-.2,(bL/2)-Cnt2>,ShowOld 
            pigment{Red}}
     #declare Cnt2=Cnt2+Cuts;
    #end
    #declare Cnt=Cnt+Cuts;
    #end
#end    


#ifdef(SavMesh)
 #fopen MyFile SavMesh write
 #write (MyFile "mesh2{\n")
 #write (MyFile "vertex_vectors {\n",aL*bL,",")
    #declare Cnt=0;
    #declare Cnt2=0;
    #while(Cnt<aL*bL)
     #if (Cnt!=aL*bL-1)
     #write (MyFile "<",vstr(3,Pnts[Cnt] ,",",0,9),">,")            
     #else
     #write (MyFile "<",vstr(3,Pnts[Cnt],",",0,9),">")            
     #end
     #if (Cnt2=3)
      #declare Cnt2=0;
       #write (MyFile "\n")            
    #end
     #declare Cnt2=Cnt2+1;
     #declare Cnt=Cnt+1;
  #end
 #write (MyFile "}\n")


 #write (MyFile "normal_vectors {\n",aL*bL,",")            
        #declare Cnt=0;
        #declare Cnt2=1;
    #declare Cnt=0;
    #while(Cnt<aL*bL)
  #while(Cnt<aL*bL)
     //PntAr[Cnt] 
   #if (Cnt!=aL*bL-1)
   #write (MyFile "<",vstr(3,Norm[Cnt],",",0,9),">,")            
   #else
   #write (MyFile "<",vstr(3,Norm[Cnt],",",0,9),">")            
   #end
     #if (Cnt2=5)
      #declare Cnt2=0;
      #write (MyFile "\n")            
     #end
   #declare Cnt2=Cnt2+1;
   #declare Cnt=Cnt+1;
  #end
 #write (MyFile "}\n")
#end

 #write (MyFile "face_indices {\n",2*A*B,",")            
        #declare Cnt=0;
        #declare Cnt2=1;
  #while(Cnt<2*A*B)
     //PntAr[Cnt] 
   #if (Cnt!=2*A*B-1)
   #write (MyFile "<",vstr(3,Facs[Cnt],",",0,0),">,")            
   #else
   #write (MyFile "<",vstr(3,Facs[Cnt],",",0,0),">")            
   #end
     #if (Cnt2=5)
      #declare Cnt2=0;
       #write (MyFile "\n")            
     #end
   #declare Cnt2=Cnt2+1;
   #declare Cnt=Cnt+1;
  #end
 #write (MyFile "}\n")
 #write (MyFile "}\n")
#end

//----------------------show and tell----------------------------------------------------------

#declare Cat=(bL/2)*tan(radians(67));//#debug concat("Cat =",str(Cat,3,3),"\n\n\n")

camera{ location <0,Cat,0>
        right x*image_width/image_height
        look_at <0,0,0>
        sky z
        }

light_source{ <Cat,Cat,-Cat> color White}
light_source{ <0,Cat,0> color White}