private import math;
import three;

real epsilon=10*realEpsilon;

// Routines for hidden surface removal (via binary space partition):
// Structure face is derived from picture.
struct face {
  picture pic;
  transform t;
  frame fit;
  triple normal,point;
  triple min,max;
  void operator init(path3 p) {
    this.normal=normal(p);
    if(this.normal == O) abort("path is linear");
    this.point=point(p,0);
    min=min(p);
    max=max(p);
  }
  face copy() {
    face f=new face;
    f.pic=pic.copy();
    f.t=t;
    f.normal=normal;
    f.point=point;
    f.min=min;
    f.max=max;
    add(f.fit,fit);
    return f;
  }
}

picture operator cast(face f) {return f.pic;}
face operator cast(path3 p) {return face(p);}

struct line {
  triple point;
  triple dir;
}

private line intersection(face a, face b)
{
  line L;
  L.point=intersectionpoint(a.normal,a.point,b.normal,b.point);
  L.dir=unit(cross(a.normal,b.normal));
  return L;
}

struct half {
  pair[] left,right;

  // Sort the points in the pair array z according to whether they lie on the
  // left or right side of the line L in the direction dir passing through P.
  // Points exactly on L are considered to be on the right side.
  // Also push any points of intersection of L with the path operator --(... z)
  // onto each of the arrays left and right.
  void operator init(pair dir, pair P ... pair[] z) {
    pair lastz;
    pair invdir=dir != 0 ? 1/dir : 0;
    bool left,last;
    for(int i=0; i < z.length; ++i) {
      left=(invdir*z[i]).y > (invdir*P).y;
      if(i > 0 && last != left) {
        pair w=extension(P,P+dir,lastz,z[i]);
        this.left.push(w);
        this.right.push(w);
      }
      if(left) this.left.push(z[i]);
      else this.right.push(z[i]);
      last=left;
      lastz=z[i];
    }
  }
}

struct splitface {
  face back,front;
}

// Return the pieces obtained by splitting face a by face cut.
splitface split(face a, face cut, projection P)
{
  splitface S;

  void nointersection() {
    if(abs(dot(a.point-P.camera,a.normal)) >=
       abs(dot(cut.point-P.camera,cut.normal))) {
      S.back=a;
      S.front=null;
    } else {
      S.back=null;
      S.front=a;
    }
  }

  if(P.infinity) {
    P=P.copy();
    static real factor=1/sqrtEpsilon;
    P.camera *= factor*max(abs(a.min),abs(a.max),
                           abs(cut.min),abs(cut.max));
  }

  if((abs(a.normal-cut.normal) < epsilon ||
      abs(a.normal+cut.normal) < epsilon)) {
    nointersection();
    return S;
  }

  line L=intersection(a,cut);

  if(dot(P.camera-L.point,P.camera-P.target) < 0) {
    nointersection();
    return S;
  }

  pair point=a.t*project(L.point,P);
  pair dir=a.t*project(L.point+L.dir,P)-point;
  pair invdir=dir != 0 ? 1/dir : 0;
  triple apoint=L.point+cross(L.dir,a.normal);
  bool left=(invdir*(a.t*project(apoint,P))).y >= (invdir*point).y;

  real t=intersect(apoint,P.camera,cut.normal,cut.point);
  bool rightfront=left ^ (t <= 0 || t >= 1);

  face back=a, front=a.copy();
  pair max=max(a.fit);
  pair min=min(a.fit);
  half h=half(dir,point,max,(min.x,max.y),min,(max.x,min.y),max);
  if(h.right.length == 0) {
    if(rightfront) front=null;
    else back=null;
  } else if(h.left.length == 0) {
    if(rightfront) back=null;
    else front=null;
  }
  if(front != null)
    clip(front.fit,operator --(... rightfront ? h.right : h.left)--cycle,
         zerowinding);
  if(back != null)
    clip(back.fit,operator --(... rightfront ? h.left : h.right)--cycle,
         zerowinding);
  S.back=back;
  S.front=front;
  return S;
}

// A binary space partition
struct bsp
{
  bsp back;
  bsp front;
  face node;

  // Construct the bsp.
  void operator init(face[] faces, projection P) {
    if(faces.length != 0) {
      this.node=faces.pop();
      face[] front,back;
      for(int i=0; i < faces.length; ++i) {
        splitface split=split(faces[i],this.node,P);
        if(split.front != null) front.push(split.front);
        if(split.back != null) back.push(split.back);
      }
      this.front=bsp(front,P);
      this.back=bsp(back,P);
    }
  }

  // Draw from back to front.
  void add(frame f) {
    if(back != null) back.add(f);
    add(f,node.fit,group=true);
    if(labels(node.fit)) layer(f); // Draw over any existing TeX layers.
    if(front != null) front.add(f);
  }
}

void add(picture pic=currentpicture, face[] faces,
         projection P=currentprojection)
{
  int n=faces.length;
  face[] Faces=new face[n];
  for(int i=0; i < n; ++i)
    Faces[i]=faces[i].copy();

  pic.add(new void (frame f, transform t, transform T,
                    pair m, pair M) {
            // Fit all of the pictures so we know their exact sizes.
            face[] faces=new face[n];
            for(int i=0; i < n; ++i) {
              faces[i]=Faces[i].copy();
              face F=faces[i];
              F.t=t*T*F.pic.T;
              F.fit=F.pic.fit(t,T*F.pic.T,m,M);
            }

            bsp bsp=bsp(faces,P);
            if(bsp != null) bsp.add(f);
          });

  for(int i=0; i < n; ++i) {
    picture F=Faces[i].pic;
    pic.userBox3(F.userMin3(), F.userMax3());
    pic.bounds.append(F.T, F.bounds);
    // The above 2 lines should be replaced with a routine in picture which
    // copies only sizing data from another picture.
  }
}