import three; import math;
size(8cm,0);
currentprojection=obliqueX;

real h=2;
triple A =(0,0,h), B=(h,0,0), C=(0,h,0), D=(0,0,0);
triple Ip=midpoint(A--C), J=midpoint(A--B);
triple K=shift((0,0,-0.25*h))*A;

triple M=interp(K,J,intersect(K,J,normal(new triple[]{B,C,D}),D));
triple Np=interp(K,Ip,intersect(K,Ip,normal(new triple[]{B,C,D}),D));

dot("$A$", A, align=N);  dot("$B$", B, align=S);
dot("$C$", C, align=S);  dot("$D$", D, align=W);
dot("$I$", Ip, align=N); dot("$J$", J, align=W);
dot("$K$", K, align=NE); dot("$M$", M, align=SE);
dot("$N$", Np, align=S);

draw(A--B--C--cycle3^^B--M^^C--Np^^J--M^^Ip--Np);
draw(A--D--C^^D--B^^J--K^^K--Ip, dashed);

import three;
size(4cm,0);

currentprojection=perspective((45,45,30));
path3 YZ=plane((0,4,0),(0,0,4));

draw("$x$",O--X,Arrow);
draw("$y$",O--Y,Arrow);
draw("$z$",O--Z,Arrow);
draw(YZ);

label(scale(5)*project("A",Y,Z,(0,2,2)));

//Author Jens Schwaiger

// For reading (simple) obj-files to polyhedrons
// `simple' means that 1) all vertex statements should come before the face statements,
//                        that 2) face informations with respect to texture and/or normal vectors are ignored and
//                        that 3) face statements only contain positive numbers (no relative positions).
//  The reading process only takes into account lines starting with "v" or "f"
import three;
import light;


string truncate(string s, string mark="/"){ // cuts string s after the first appearence of the string mark
  int pos=find(s,mark);
  if(pos==-1) return s;
  return erase(s,pos,-1);
}

string[] truncate(string[] s, string mark="/"){ // extension to arrays of strings
  if(s.length==0) abort("array has length 0");
  string[] tmp;
  for(int i=0;i<s.length;++i) tmp[i]=truncate(s[i],mark);
  return tmp;
}

typedef guide3[] polyhedron;

real sf(polyhedron p){
  if(p.length==0) abort("empty polyhedron");
  real maxmin(int i){
    return max(abs(max(p[i])),abs(min(p[i])));
  }
  return 10*max(sequence(maxmin,p.length));
}

triple cofgr(guide3 g){
  int l=length(g);
  if(l==0) return point(g,0);
  triple cgr=(0,0,0);
  for(int i=0;i<l;++i) cgr=cgr+point(g,i);
  return(1/l)*cgr;
}


polyhedron operator *(transform3 T, polyhedron p){
  polyhedron pt;
  for(int i=0;i<p.length;++i) pt[i]=T*p[i];
  return pt;
}

void draw(picture pic=currentpicture, polyhedron[] pol, pen[] pn={currentpen}, projection P=currentprojection){
  pen[] pntmp=pn;
  pntmp.cyclic(true);
  for(int i=0;i<pol.length;++i) for(int j=0;j<(pol[i]).length;++j) draw(pic,project(pol[i][j],P),pntmp[i]);
}
void filldraw(picture pic=currentpicture,
              polyhedron[] pol,
              pen[] fcol={ red,blue,green,orange, gray},
              real op=0.75,
              pen[] dcol={black+1bp},
              projection P=currentprojection,
              light lt=currentlight)
{
  write("Sorting the faces with respect to their distance from the camera.");
  fcol.cyclic(true);
  dcol.cyclic(true);
  polyhedron join;
  int[] nofpoly;
  int[] constarray(int i, int n){int[] cr; for(int j=0;j<n;++j) cr[j]=i; return cr;}
  for(int i=0;i<pol.length;++i){join.append(pol[i]);nofpoly.append(constarray(i,pol[i].length));}
  triple cam;
  cam=P.camera;
  if(P.infinity) cam=sf(join)*P.camera;
  real[][] dist=new real[join.length][3];
  real[][] distsort=new real[join.length][3];
  for(int i=0;i<join.length;++i){dist[i][0]=abs(cofgr(join[i])-cam);dist[i][1]=i;dist[i][2]=nofpoly[i];}
  distsort=sort(dist);
  write("Sorting done.");
  write("Filling the faces.");
  for(int i=distsort.length-1;i>-1;--i){
    int d1=(int)distsort[i][1];
    int d2=(int)distsort[i][2];
    triple u,v,w;
    u=point(join[d1],0)-cofgr(join[d1]);
    v=point(join[d1],1)-cofgr(join[d1]);
    w=unit(cross(u,v));
    if(dot(w,cam-cofgr(join[d1]))<0) w=-w;
    real contrast=0.5*(dot(lt.source,w)+1);
    filldraw(pic,project(join[d1],P),contrast*fcol[d2]+opacity(op),dcol[d2]);
  }
  write("Filling done.");
}

polyhedron objtopolyhedron(string datafile){
  file in=word(line(input(datafile,check=true)));

  triple[] vert;
  polyhedron pol;
  void Vertex(real x,real y ,real z){vert.push((x,y,z));}
  void Face(int[] vertnr){
    guide3 gh;
    for(int i=0;i<vertnr.length;++i) gh=gh--vert[vertnr[i]-1];
    gh=gh--cycle;
    pol.push(gh);
  }
  write("Reading data from "+datafile+".");
  while(true) {
    string[] str=in;
    if(str.length==0) break;
    str=truncate(str);
    if(str[0]=="v") Vertex((real) str[1],(real) str[2],(real) str[3]);
    if(str[0]=="f") {
      int[] vertnr;
      for(int i=1;i<str.length;++i) vertnr[i-1]=(int) str[i];
      Face(vertnr);
    }

    if(eof(in)) break;
  }
  close(in);
  write("Reading done.");
  return pol;
}


size(12cm,true);
currentprojection=orthographic(5,4,1);
currentlight=(-5,5,7.5);

// More data may be found at http://www.cs.technion.ac.il/~irit/data/Viewpoint/

polyhedron[] parr={rotate(90,Y)*rotate(90,Z)*objtopolyhedron("galleon.obj")};
pen[] fcol1={rgb(0.5,0.2,0.2)};
pen[] dcol1={darkgrey};
filldraw(parr,fcol=fcol1,dcol=dcol1, op=1.0);
shipout(bbox(3mm,darkblue+3bp+miterjoin,FillDraw(paleyellow)));

//Author Jens Schwaiger
// For reading (simple) obj-files with groups to polyhedrons
// `simple' means that 1) all vertex statements should come before the face statements,
//                        that 2) face informations with respect to texture and/or normal vectors are ignored and
//                         that 3) face statements only contain positive numbers (no relative positions).
//  The reading process only takes into account lines starting with "v" or "f" or "g" (group)
import three;
import light;

string truncate(string s, string mark="/"){ // cuts string s after the first appearence of the string mark
  int pos=find(s,mark);
  if(pos==-1) return s;
  return erase(s,pos,-1);
}

string[] truncate(string[] s, string mark="/"){ // extension to arrays of strings
  if(s.length==0) abort("array has length 0");
  string[] tmp;
  for(int i=0;i<s.length;++i) tmp[i]=truncate(s[i],mark);
  return tmp;
}

int find(string[] S, string s){
  bool[] test;
  for(int i=0;i<S.length;++i) test[i]=(S[i]==s);
  return find(test);
}

typedef guide3[] polyhedron;

real sf(polyhedron p){
  real maxmin(int i){return max(abs(max(p[i])),abs(min(p[i])));};
  if(p.length==0) abort("empty polyhedron");
  return 100*max(sequence(maxmin,p.length));
};

triple cofgr(guide3 g){
  int l=length(g);
  if(l==0) return point(g,0);
  triple cgr=(0,0,0);
  for(int i=0;i<l;++i) cgr=cgr+point(g,i);
  return(1/l)*cgr;
};

polyhedron operator *(transform3 T, polyhedron p){polyhedron pt;
  for(int i=0;i<p.length;++i) pt[i]=T*p[i];
  return pt;
}

polyhedron[] operator *(transform3 T, polyhedron[] p){polyhedron[] pt;
  for(int i=0;i<p.length;++i) pt[i]=T*p[i];
  return pt;
}
void draw(picture pic=currentpicture, polyhedron[] pol, pen[] pn={currentpen}, projection P=currentprojection){
  pen[] pntmp=pn;
  pntmp.cyclic(true);
  for(int i=0;i<pol.length;++i) for(int j=0;j<(pol[i]).length;++j) draw(pic,project(pol[i][j],P),pntmp[i]);
};
void filldraw(picture pic=currentpicture,
              polyhedron[] pol,
              pen[] fcol={ red,blue,green,orange, gray},
              real op=0.75,
              pen[] dcol={black+1bp},
              projection P=currentprojection,
              light lt=currentlight)
{
  write("Sorting the faces with respect to their distance from the camera.");
  fcol.cyclic(true);
  dcol.cyclic(true);
  polyhedron join;
  int[] nofpoly;
  int[] constarray(int i, int n){int[] cr; for(int j=0;j<n;++j) cr[j]=i; return cr;};
  for(int i=0;i<pol.length;++i){join.append(pol[i]);nofpoly.append(constarray(i,pol[i].length));};
  triple cam;
  cam=P.camera;
  if(P.infinity) cam=sf(join)*P.camera;
  real[][] dist=new real[join.length][3];
  real[][] distsort=new real[join.length][3];
  for(int i=0;i<join.length;++i){dist[i][0]=abs(cofgr(join[i])-cam);dist[i][1]=i;dist[i][2]=nofpoly[i];};
  //write("Sorting the faces with respect to their distance from the camera.");
  distsort=sort(dist);
  write("Sorting done.");
  write("Filling the faces.");
  for(int i=distsort.length-1;i>-1;--i){
    int d1=(int)distsort[i][1];
    int d2=(int)distsort[i][2];
    triple u,v,w;
    u=point(join[d1],0)-cofgr(join[d1]);
    v=point(join[d1],1)-cofgr(join[d1]);
    w=unit(cross(u,v));
    if(dot(w,cam-cofgr(join[d1]))<0) w=-w;
    real contrast=0.5*(dot(lt.source,w)+1);
    filldraw(pic,project(join[d1],P),contrast*fcol[d2]+opacity(op),dcol[d2]);
  };
  write("Filling done.");
};

polyhedron[] objtopolyhedrongroup(string datafile){
  file in=word(line(input(datafile,check=true)));
  triple[] vert;
  polyhedron[] polh=new polyhedron[];
  polh[0]=new guide3[];
  string[] G;
  void Vertex(real x,real y ,real z){vert.push((x,y,z));};
  void Face(int[] vertnr, int groupnr){
    guide3 gh;
    for(int i=0;i<vertnr.length;++i) gh=gh--vert[vertnr[i]-1];
    gh=gh--cycle;
    polh[groupnr].push(gh);
  };
  write( "Reading data from "+datafile+". This may take a while."       );
  int groupnr;
  while(true) {
    string[] str=in;
    if(str.length==0) break;
    str=truncate(str);
    if(str[0]=="g" && str.length>1) {
      int tst=find(G,str[1]);
      if(tst==-1)        {G.push(str[1]);groupnr=G.length-1; polh[groupnr]=new guide3[];};
      if(tst>-1) groupnr=tst;
    };
    if(str[0]=="v") Vertex((real) str[1],(real) str[2],(real) str[3]);
    if(str[0]=="f") {
      int[] vertnr;
      for(int i=1;i<str.length;++i) vertnr[i-1]=(int) str[i];
      Face(vertnr,groupnr);
    };
    if(eof(in)) break; };
  write("Number of groups: ", (string) G.length);
  write("Groups and their names");
  write(G);
  close(in);
  write("Reading done.");
  write("Number of faces contained in the groups:");
  for(int j=0;j<polh.length;++j) write(G[j],": ",(string) polh[j].length);
  return polh;
};

size(15cm,true);
currentprojection=orthographic(5,4,1);
//currentprojection=perspective(10,24,-1);
//currentprojection=perspective(10,24,3);
//currentprojection=perspective(10,24,-20);
//currentprojection=perspective(50,40,10);
currentlight=(1,1,1);
//currentlight=(0.1,0,0);
//currentlight=(0.5,0.25,1);
//currentlight=(0.3,0.1,0.1);

/*More interesting data, e.g., triceratops.obj, may be found at HERE*/
polyhedron[] parr=rotate(65,Z)*rotate(90,Y)*rotate(90,Z)*objtopolyhedrongroup("galleon.obj");

pen[] fcol1={
  darkred,brown,darkred+orange,heavyred,heavyred,darkred+orange,palegreen+blue+lightgrey,
  darkred,darkred,yellow,darkred,white,white,white,white,white,white};
// Information for defining fcol1 is taken from a first run and reading the number of groups and their names
pen[] dcol1={darkgrey};

filldraw(parr,fcol=fcol1,dcol=dcol1, op=1.0);
shipout(bbox(3mm,darkblue+3bp+miterjoin,FillDraw(fillpen=paleblue)));

// From the Asymptote forum
import three;
import light;

typedef path3[] shape;

shape operator *(transform3 T, shape p){
  shape os;
  for(path3 g:p) os.push(T*g);
  return os;
}


path3 path(triple[] T){
  path3 P;
  for(triple i:T) P=P--i;
  return P;
}

void addshapes(face[] F, shape[] shp, pen drawpen=currentpen, pen fillpen=white)
{
  for(int i=0; i < shp.length; ++i)
    for(int j=0; j < shp[i].length; ++j) {
      path3 g=shp[i][j];
      picture pic=F.push(g);
      filldraw(pic,g,currentlight.intensity(F[F.length-1].point)*fillpen, drawpen);
    }
}

shape cylinder(real R=1, real H=1, int n=18){
  shape Cyl;
  triple[] CTop;
  triple[] CBot;
  for(int i=0; i <= n; ++i)
    CBot.push((R*cos(2pi*i/n), R*sin(2pi*i/n),0));
  CTop = CBot+(0,0,H);
  for(int i=0; i < n; ++i)
    Cyl.push(CBot[i]--CBot[i+1]--CTop[i+1]--CTop[i]--cycle3);

  path3 P=path(CBot)--cycle3;
  Cyl.push(P);
  Cyl.push(shift(H*Z)*P);

  return Cyl;
}

shape rightslab(real x=1, real y=1, real z=1){
  shape slab;
  slab[0] = (0,0,0)--(1,0,0)--(1,1,0)--(0,1,0)--cycle3;
  slab[1] = (0,0,0)--(1,0,0)--(1,0,1)--(0,0,1)--cycle3;
  slab[2] = (1,0,0)--(1,1,0)--(1,1,1)--(1,0,1)--cycle3;
  slab[3] = (1,1,0)--(0,1,0)--(0,1,1)--(1,1,1)--cycle3;
  slab[4] = (0,1,0)--(0,0,0)--(0,0,1)--(0,1,1)--cycle3;
  slab[5] = (0,0,1)--(1,0,1)--(1,1,1)--(0,1,1)--cycle3;
  return scale(x,y,z)*slab;
}

size(10cm,0);
triple cam=(1600,200,150);
//currentprojection=orthographic(1600,800,400);
currentprojection=perspective(cam); //Far away!
currentlight=rotate(-45,Z)*(cam+(0,0,1000));

real Blen = 180;
real Bwdt = 30;
real Bhgt = 3;
real Clen = 130;
real Cwdt = 50;
real Chgt = 50;
real cylr = 7.5;
real cylh = 37.0 ;

shape slab1 = shift(-Bwdt/2*Y-Bhgt/2*Z+Clen/2*X)*rightslab(Blen,Bwdt,Bhgt);
shape slab2 = shift(-Cwdt/2*Y-Chgt/2*Z-Clen/2*X)*rightslab(Clen,Cwdt,Chgt);
shape cyl1 = shift((Blen+Clen/2-2*cylr)*X-(cylh/2)*Z)*cylinder(R=cylr, H=cylh);

shape[] group1={slab1};
shape[] group2={slab2};
shape[] group3={cyl1};

face[] faces;
addshapes(faces, group1, drawpen=linewidth(0.25bp), fillpen=opacity(0.35)+red);
addshapes(faces, group2, drawpen=linewidth(0.25bp), fillpen=opacity(0.5)+lightblue);
addshapes(faces, group3, drawpen=linewidth(0.25bp), fillpen=opacity(0.5)+lightyellow);
add(faces);

shipout(format="pdf", bbox(3mm,white));