import java.util.Vector;
import java.util.Iterator;

public class Contours {

  public static Vector contour(Vector triangles, double level) {
    Vector contours = new Vector();
    Iterator iter = triangles.iterator();
    while(iter.hasNext()) {
      Triangle t = (Triangle)iter.next();
      Edge e = getContour(t,level);
      if (e != null) {
        contours.add(e); 
      }
    }
    return contours;
  }

  /*-------------------------------------------------------------------------
    Create a contour slice through a triangle t
    The contour "level" is a horizontal plane perpendicular to the z axis,
    ie: The equation of the contour plane Ax + By + Cz + D = 0
    has A = 0, B = 0, C = 1, D = -level
    Return
      null if the contour plane doesn't cut the facet
      Edge of line points if it does cut the facet
      null for an unexpected occurrence 
    If a vertex touches the contour plane nothing needs to be drawn!?
    adapted from code by Paul Bourke - http://astronomy.swin.edu.au/~pbourke/projection/contour/index.html
   */
  public static Edge getContour(Triangle t, double level) {

    /*
      Evaluate the equation of the plane for each vertex
     sidea = A * t.p1.x + B * t.p1.y + C * t.p1.z + D;
     sideb = A * t.p2.x + B * t.p2.y + C * t.p2.z + D;
     sidec = A * t.p3.x + B * t.p3.y + C * t.p3.z + D;
     */
    double sidea = t.p1.z - level;
    double sideb = t.p2.z - level;
    double sidec = t.p3.z - level;

    /* Are all the vertices on one side */
    if (SIGN(sidea) == SIGN(sideb) && SIGN(sidea) == SIGN(sidec)) {
      return null;
    }

    /* Is p0 the only point on a side by itself */
    if ((SIGN(sidea) != SIGN(sideb)) && (SIGN(sidea) != SIGN(sidec))) {
      Edge e = new Edge();
      e.p1 = new Point3d( t.p1.x - sidea * (t.p3.x - t.p1.x) / (sidec - sidea), t.p1.y - sidea * (t.p3.y - t.p1.y) / (sidec - sidea), t.p1.z - sidea * (t.p3.z - t.p1.z) / (sidec - sidea) );
      e.p2 = new Point3d( t.p1.x - sidea * (t.p2.x - t.p1.x) / (sideb - sidea), t.p1.y - sidea * (t.p2.y - t.p1.y) / (sideb - sidea), t.p1.z - sidea * (t.p2.z - t.p1.z) / (sideb - sidea) );
      return e;
    }

    /* Is p1 the only point on a side by itself */
    if ((SIGN(sideb) != SIGN(sidea)) && (SIGN(sideb) != SIGN(sidec))) {
      Edge e = new Edge();
      e.p1 = new Point3d( t.p2.x - sideb * (t.p3.x - t.p2.x) / (sidec - sideb), t.p2.y - sideb * (t.p3.y - t.p2.y) / (sidec - sideb), t.p2.z - sideb * (t.p3.z - t.p2.z) / (sidec - sideb) );
      e.p2 = new Point3d( t.p2.x - sideb * (t.p1.x - t.p2.x) / (sidea - sideb), t.p2.y - sideb * (t.p1.y - t.p2.y) / (sidea - sideb), t.p2.z - sideb * (t.p1.z - t.p2.z) / (sidea - sideb) );
      return e;
    }

    /* Is p2 the only point on a side by itself */
    if ((SIGN(sidec) != SIGN(sidea)) && (SIGN(sidec) != SIGN(sideb))) {
      Edge e = new Edge();
      e.p1 = new Point3d( t.p3.x - sidec * (t.p1.x - t.p3.x) / (sidea - sidec), t.p3.y - sidec * (t.p1.y - t.p3.y) / (sidea - sidec), t.p3.z - sidec * (t.p1.z - t.p3.z) / (sidea - sidec) );
      e.p2 = new Point3d( t.p3.x - sidec * (t.p2.x - t.p3.x) / (sideb - sidec), t.p3.y - sidec * (t.p2.y - t.p3.y) / (sideb - sidec), t.p3.z - sidec * (t.p2.z - t.p3.z) / (sideb - sidec) );
      return e;
    }

    /* Shouldn't get here */
    System.err.println("strangeness in contouring");
    return null;
  }

  public static int SIGN(double a) {
    return (a > 0) ? 1 : -1;
  }

}