Return to povray.rr CVS log

Up to [local] / OpenXM / src / asir-contrib / packages / src

Experimental version of sm1 interface of asir-contrib.
sm1_xyz() ===>  sm1.xyz().

Warning: Porting asir-contrib for the new "modular" asir is a work in progress.
asir-contrib will not work perfectly if you cvsup this change.
Checkout asir-contrib with the tag RELEASE_1_2_2 if you need to USE
asir-contrib.

/* $OpenXM: OpenXM/src/asir-contrib/packages/src/povray.rr,v 1.8 2003/05/18 02:20:22 takayama Exp $ */

/* This is an experimental asir program to model
   z=f(x,y) for povray.
*/

/* This has not yet been ported to be used with the new sm1 for module. 
   2003.05.18 */

/* 
References:
  povray implicit algebraic surface
  Polyray system   (this is an excellent cite).
  http://wims.unice.fr/~wims/en_tool~geometry~polyray.en.html
  http://pages.infinit.net/gollum/polyray/poly1.htm

  All constrains is expressed by F(x1,x2,...,xn)>=0
  http://www.hyperfun.org
     (http://wwwcis.k.hosei.ac.jp/~F-rep/HF_proj.html)

  http://www.povray.org/links/
    regular 3D polygon generator.
    (3D polygon) http://www.swin.edu.au/astronomy/pbourke/geometry/polyhedra/
    
  A tutorial:
    http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/surface/basic.html

  Introduction in Japanese.
    http://www.jks.is.tsukuba.ac.jp/~akaho/cg/3d/povray.html

  Check out also surf.
*/

/* A note for Debian
    apt-get install povray
    apt-get install povray-misc
    apt-get install povray-manual

*/

/* File IO functions. These should be implemented as an asir
   native functions.
*/
File_file_name = 0$
def open_file_for_write(FileName) {
  extern File_file_name, Sm1_proc;
  File_file_name = FileName;
  S = " ("+FileName+") (w) file /asir_file_tmp set ";
  if (Sm1_proc < 0) sm1_start();
  print(S);
  sm1(Sm1_proc,S);
  return(FileName);
}

def put_bytes(Fp,S) {
  extern Sm1_proc;
  if (type(S) == 0 || type(S) == 1) {
    S = rtostr(S);
    S = " asir_file_tmp "+S+" (string) dc "+" writestring ";
    sm1(Sm1_proc,S);
  }else if (type(S) == 7) {
    sm1(Sm1_proc," asir_file_tmp ");
    ox_push_cmo(Sm1_proc,S);
    sm1(Sm1_proc," writestring ");
  }else error("cannot write this object.");
}

def close_file_for_write(Fp) {
  extern Sm1_proc;
  sm1(Sm1_proc," asir_file_tmp closefile 0");
  ox_pop_cmo(Sm1_proc); /* Wait for closing the file. */
}


/* Generating a shade file */
def povray_prolog(FileName,Camera, Light) {
  /* Introduction to povray:
     See for example, http://www.jks.is.tsukuba.ac.jp/~akaho/cg/3d   
     The next prolog is taken from this cite.
  */
  NL = asciitostr([10]);
  Fp = open_file_for_write(FileName);
  S = "#version 3.0 "+NL+
      "global_settings { assumed_gamma 2.2 } "+NL+
      "#include \"colors.inc\" "+NL+
      "background { color Blue } "+NL+
      "camera { "+NL+
/*      "location <0, 5, -30> "+NL+ */
      "location <10, 5, -20> "+NL+ 
      "right <4/3, 0, 0> "+NL+
      "up <0, 1, 0> "+NL+
      "sky <0, 1, 0> "+NL+
      "direction <0, 0, 1.21> "+NL+
      "look_at <0, 0, 0> "+NL+
      "} "+NL+
      "light_source { <5, 20, -10> colour White } "+NL+
      "#declare RED = texture {"+NL+
      "   pigment {color rgb<0.8, 0.2,0.2>}"+NL+
      "   finish { ambient 0.2 diffuse 0.5}"+NL+
      "}"+NL+
      "#declare GREEN = texture {"+NL+
      "   pigment {color rgb<0.8, 0.8, 0.2>}"+NL+
      "   finish { ambient 0.2 diffuse 0.5}"+NL+
      "}"+NL$
 

  put_bytes(Fp,S);
  return(Fp);
}

def povray_epilog(Fp) {
  close_file_for_write(Fp);
}

def povray_spheres(Fp,Slist) {
  NL = asciitostr([10]);
  S = " union { "+NL;
  put_bytes(Fp,S);
  S = " ";
  Param = 
   "pigment { "+NL+
   "color Silver "+NL+
   "filter 0.0 "+NL+
   "} "+NL+
   "finish { "+NL+
   "ambient 0.2 "+NL+
   "diffuse 0.5 "+NL+
   "phong 0.6 "+NL+
   "phong_size 7 "+NL+
   "reflection 0.8 "+NL+
   "} "+NL$

  N = length(Slist);
  for (I=0; I<N; I++) {
    Sphere = Slist[I];
    S = " sphere { <"+rtostr(Sphere[0])+","+rtostr(Sphere[1])+","+
                      rtostr(Sphere[2])+">,"+rtostr(Sphere[3])+NL;
    put_bytes(Fp,S);
    put_bytes(Fp,Param);
    put_bytes(Fp," } "+NL+NL);
    put_bytes(Fp,10);
  }
  put_bytes(Fp,"} "+NL);
}

/* Example input

    povray_prolog("t.pov",0,0);
    povray_spheres(0,[ [0,3,0,3],[1,1,0,2] ]);
    povray_epilog(0);

*/

def povray_plane(Fp,Dist) {
NL=asciitostr([10]);
S=
"plane {"+NL+
"  <0,1,0>,"+rtostr(Dist)+NL+
"  pigment {"+NL+
"    checker"+NL+
"      color Yellow"+NL+
"      color Blue"+NL+
"    scale 5"+NL+
"  }"+NL+
"  finish {"+NL+
"    ambient 0.2"+NL+
"    diffuse 0.8"+NL+
"    phong 0.6"+NL+
"    phong_size 7"+NL+
"    reflection 0.1"+NL+
"  }"+NL+
"}"+NL$
  put_bytes(0,S);
}

/*  Generate mesh objects */

def povray_mesh(Fp,Slist) {
  /* Slist is a union of triangles.
     Example:
     [ [[0,0,0],[1,0,0],[0,1,0]],  
       [[1,0,0],[0,1,0],[1,1,0.5]]
     ]
  */
  NL = asciitostr([10]);
  S = "mesh { "+NL;
  put_bytes(Fp,S);
  S = " ";
  Param = 
   " pigment { "+NL+
   "   color Silver "+NL+
   "   filter 0.0 "+NL+
   " } "+NL+
   " finish { "+NL+
   "   ambient 0.2 "+NL+
   "   diffuse 0.5 "+NL+
   "   phong 0.6 "+NL+
   "   phong_size 7 "+NL+
   "   reflection 0.8 "+NL+
   " } "+NL$

  N = length(Slist);
  for (I=0; I<N; I++) {
    Mesh = Slist[I];
    S = " triangle { <"+rtostr(Mesh[0][0])+","+rtostr(Mesh[0][1])+","+
                        rtostr(Mesh[0][2])+">,"+
                    "<"+rtostr(Mesh[1][0])+","+rtostr(Mesh[1][1])+","+
                        rtostr(Mesh[1][2])+">,"+
                    "<"+rtostr(Mesh[2][0])+","+rtostr(Mesh[2][1])+","+
                        rtostr(Mesh[2][2])+"> "+NL$
    if (P==0) {
      S = S+ "  texture { RED } "+NL$
      P = 1;
    }else{
      S = S+ "  texture { GREEN } "+NL$
      P = 0;
    }
    S = S+ " } "+NL;
    put_bytes(Fp,S);
  }
  put_bytes(Fp,Param);
  put_bytes(Fp,"} "+NL);
}

def test_mesh() {
    povray_prolog("t.pov",0,0);
    povray_mesh(0,[ [[0,0,0],[1,0,0],[0,1,0]],    /* triangle 1 */
                    [[1,0,0],[0,1,0],[1,1,0.5]]   /* triangle 2 */
                  ]);
    povray_epilog(0);
}

/* test mesh */
def test_f(X,Y) {
  return(deval((X^2-(Y-3)^2)/4));
}
def test_triangle1(X,Y) {
  DX = 0.5$ DY=0.5$  /* Use camera at <10,5,-20> */
  return([
          [X,test_f(X,Y),Y],
          [X+DX,test_f(X+DX,Y),Y],
          [X,test_f(X,Y+DY),Y+DY]]);
}
def test_triangle2(X,Y) {
  DX = 0.5$ DY=0.5$  /* Use camera at <10,5,-20> */
  return([
          [X,test_f(X,Y),Y],
          [X-DX,test_f(X-DX,Y),Y],
          [X,test_f(X,Y-DY),Y-DY]]);
}
def test_mesh2() {
    DX = 0.5$ DY=0.5$  /* Use camera at <10,5,-20> */
    Mesh=[];  
    for (X=-5; X<5; X += DX) {
      for (Y=0; Y<5; Y += DY) {
         Mesh = append(Mesh,[test_triangle1(X,Y),test_triangle2(X+DX,Y+DY)]);
      }
    }
    povray_prolog("t.pov",0,0);
    povray_plane(0,0);
    povray_mesh(0,Mesh);
    povray_epilog(0);
    print("See the picture by  ")$
    print("   povray +it.pov +D0 ")$
    print("or ")$
    print("   povray +it.pov +L/usr/local/lib/povray3/include +D0 +P +H200 +W200")$
    povray_render("t.pov");
}

def povray_render(F)  {
   /* If FreeBSD, */
   shell("   povray +i"+F+" +L/usr/local/lib/povray3/include +D0 +P +H200 +W200&")$
}


/* A general function to generate meshes for z=f(x,y). */
/*  povray_plot3d_mesh(F|xrange=[xmin,xmax],yrange=[ymin,ymax],
                         z_bounding_box=[zmin,zmax],camera=[x,y,z],
                         floor=bool);
   Ex. povray_plot3d_mesh(sin(2*x)*y^2);
*/
#define XMAX 5
#define XMIN (-5)
#define YMAX 5
#define YMIN (-1)
#define ZMAX 5
#define ZMIN (-5)
Povray_bound_zmin = 2*ZMIN $  /* Need space before $. 
                                 Otherwise, some preprocessors
                                 cannot replace ZMIN by the value. */
Povray_bound_zmax = 2*ZMAX $
def povray_plot3d_mesh(F) {

    Workfile = "t.pov";
    Xmin = XMIN; Xmax = XMAX;
    Ymin = YMIN; Ymax = YMAX;
    P = getopt(floor);  /* type(getopt(floor)) does not work. */
    if (type(P) == -1) {
      Need_floor = 1;
    }else{
      if (type(P)==1 && P != 0)  {
         Need_floor = 1;
      }else{
         Need_floor = 0;
      }
    }

    Z_bound = povray_get_z_bound(F,Xmin,Xmax,Ymin,Ymax);
    print("Z_bound=",0); print(Z_bound);
    Zmin = Z_bound[0];
    Zmax = Z_bound[1];
    
    F = subst(F,x,(XMIN-XMAX)*(x-Xmax)/(Xmin-Xmax)+XMAX,
                y,(YMIN-YMAX)*(y-Ymax)/(Ymin-Ymax)+YMAX);
    F = (ZMIN-ZMAX)*(F-Zmax)/(Zmin-Zmax)+ZMAX;

    NewOrigin = [Xmax-XMAX*(Xmin-Xmax)/(XMIN-XMAX),
                 Ymax-YMAX*(Ymin-Ymax)/(YMIN-YMAX),
                 (ZMIN-ZMAX)*(-Zmax)/(Zmin-Zmax)+ZMAX];
    NewOrigin = map(deval,NewOrigin);
    Distance = NewOrigin[2];

    /* Generate mesh */
    DX = 0.5$ DY=0.5$  /* Use camera at <10,5,-20> */
    Mesh=[];  
    Range =[map(subst,
             [(XMIN-XMAX)*(x-Xmax)/(Xmin-Xmax)+XMAX,
              (YMIN-YMAX)*(y-Ymax)/(Ymin-Ymax)+YMAX],x,XMIN,y,YMIN),
            map(subst,
             [(XMIN-XMAX)*(x-Xmax)/(Xmin-Xmax)+XMAX,
              (YMIN-YMAX)*(y-Ymax)/(Ymin-Ymax)+YMAX],x,XMAX,y,YMAX)];
    print(Range); 

    for (X=XMIN; X<XMAX; X += DX) {
      for (Y=YMIN; Y<YMAX; Y += DY) {
         Mesh = append(Mesh,
                 [povray_triangle1(F,X,Y,DX,DY),
                  povray_triangle2(F,X+DX,Y+DY,DX,DY)
                 ]);
      }
    }
    povray_prolog(Workfile,0,0);
    if (Need_floor)  povray_plane(0,Distance);
    povray_mesh(0,Mesh);
    povray_spheres(0,[[NewOrigin[0],NewOrigin[2],NewOrigin[1],0.3]]); 
     /* Small sphere to show the place of the origin */
    povray_epilog(0);

    print("See the picture by  ")$
    print("   povray +i"+Workfile+" +D0 ")$
    print("or ")$
    print("   povray +i"+Workfile+" +L/usr/local/lib/povray3/include +D0 +P +H200 +W200")$
    povray_render(Workfile);
    return(Mesh);
}

def povray_get_z_bound(F,Xmin,Xmax,Ymin,Ymax) {
   Zmin = Zmax = deval(subst(F,x,(Xmax-Xmin)/2,y,(Ymax-Ymin)/2));
   DX = deval((Xmax-Xmin)/5);
   DY = deval((Ymax-Ymin)/5);
   for (X=Xmin; X<Xmax; X += DX) {
     for (Y=Ymin; Y<Ymax; Y += DY) {
       Z = deval(subst(F,x,X,y,Y));
       if (Z < Zmin) { 
         Zmin=Z;
       }else if (Z > Zmax) {
         Zmax =Z;
       }
     }
   }
   return([Zmin,Zmax]);
}

def povray_function(F,X,Y) {
   extern Povray_bound_zmin, Povray_bound_zmax;
   Z = deval(subst(F,x,X,y,Y));
   if (Z > Povray_bound_zmax) return(Povray_bound_zmax);
   if (Z < Povray_bound_zmin) return(Povray_bound_zmin);
   return(Z);
}

def povray_triangle1(F,X,Y,DX,DY) {
  return([
          [X,povray_function(F,X,Y),Y],
          [X+DX,povray_function(F,X+DX,Y),Y],
          [X,povray_function(F,X,Y+DY),Y+DY]]);
}
def povray_triangle2(F,X,Y,DX,DY) {
  return([
          [X,povray_function(F,X,Y),Y],
          [X-DX,povray_function(F,X-DX,Y),Y],
          [X,povray_function(F,X,Y-DY),Y-DY]]);
}


/* The next function patches f(x,y,z)=0 by small squares
   by using phc pack or pari(root,*) and normal vectors.
*/

/* Not yet written */



end$