/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor                            Date:    10.12.99 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Internal class for executing boolean operations           *
 *           on Polyhedra                                              *
 *                                                                     *
 ***********************************************************************/

#include <CLHEP/Geometry/Plane3D.h>

using namespace HepGeom;

#define INITIAL_SIZE 200
#define CRAZY_POINT  Point3D<double>(-10.e+6, -10.e+6, -10.e+6)
#define GRANULARITY  10.e+5;

#define SWAP(A,B) w = A; A = B; B = w

#define OP_UNION         0    // Operations
#define OP_INTERSECTION  1
#define OP_SUBTRACTION   2

#define OUT_OF_PLANE     0    // Face vs face statuses
#define ON_PLANE         1
#define INTERSECTION     2
#define EDGE             3
#define NON_PLANAR_FACE  4

#define UNKNOWN_FACE     0    // Face statuses 
#define ORIGINAL_FACE   -1
#define NEW_FACE        -2
#define UNSUITABLE_FACE -3
#define DEFECTIVE_FACE  -4

// -------------------------------------------- Simplified STL vector ---
template<class T>
class vector {
 private:
  int cur_size, max_size;
  T * v;

 public:
  vector(): cur_size(0), max_size(INITIAL_SIZE), v(new T[INITIAL_SIZE]) {}
  ~vector() { delete [] v; }  

  void      clear()                 { cur_size = 0; }
  int       size()            const { return cur_size; }
  T &       operator[](int i)       { return v[i]; }
  const T & operator[](int i) const { return v[i]; }
  T &       front()                 { return v[0]; }
  const T & front()           const { return v[0]; }
  T &       back()                  { return v[cur_size-1]; }
  const T & back()            const { return v[cur_size-1]; }
  void      pop_back()              { cur_size--; }
  void      push_back(const T & a) {
    if (cur_size == max_size) {
      T * w     = v;
      max_size *= 2;
      v         = new T[max_size];
      for (int i=0; i<cur_size; i++) v[i] = w[i];
      v[cur_size++] = a;
      delete [] w;
    }else{
      v[cur_size++] = a;
    }
  }
};

// ---------------------------------------------------- Extended node ---
class ExtNode {
 public:
  Point3D<double> v;
  int        s;

 public: 
  ExtNode(Point3D<double> vertex=Point3D<double>(), int status=0)
    : v(vertex), s(status) {}
  ~ExtNode() {}

  ExtNode(const ExtNode & node) : v(node.v), s(node.s) {}

  ExtNode & operator=(const ExtNode & node) {
    v = node.v;
    s = node.s;
    return *this;
  }
};

// ---------------------------------------------------- Extended edge ---
class ExtEdge {
 public:
  int       i1, i2;           // end points
  int       iface1;           // native face
  int       iface2;           // neighbouring face
  int       ivis;             // visibility: +1 (visible), -1 (invisible)
  int       inext;            // index of next edge 

 public:
  ExtEdge(int k1=0, int k2=0, int kface1=0, int kface2=0, int kvis=0) :
    i1(k1), i2(k2), iface1(kface1), iface2(kface2), ivis(kvis), inext(0) {}

  ~ExtEdge() {};

  ExtEdge(const ExtEdge & edge) :
    i1(edge.i1), i2(edge.i2), iface1(edge.iface1), iface2(edge.iface2),
    ivis(edge.ivis), inext(edge.inext) {}

  ExtEdge & operator=(const ExtEdge & edge) {
    i1     = edge.i1;
    i2     = edge.i2;
    iface1 = edge.iface1;
    iface2 = edge.iface2;
    ivis   = edge.ivis;
    inext  = edge.inext;
    return *this;
  }

  void invert() {
    int w;
    SWAP(i1, i2);
  }
};

// ---------------------------------------------------- Extended face ---
class ExtFace {
 public:
  int        iedges[4];        // indices of original edges
  Plane3D<double> plane;            // face plane
  double     rmin[3], rmax[3]; // bounding box
  int        iold;             // head of the list of the original edges
  int        inew;             // head of the list of the new edges
  int        iprev;            // index of previous face 
  int        inext;            // index of next face 

 public:
  ExtFace(int iedge=0) : iold(iedge), inew(0), iprev(iprev), inext(0) {}
  ~ExtFace() {}

  ExtFace(const ExtFace & face) :
    plane(face.plane), iold(face.iold), inew(face.inew),
    iprev(face.iprev), inext(face.inext)
  { 
    int i;
    for (i=0; i<4; i++) { iedges[i] = face.iedges[i]; }
    for (i=0; i<3; i++) { rmin[i] = face.rmin[i]; rmax[i] = face.rmax[i]; }
  }

  ExtFace & operator=(const ExtFace & face) {
    int i;
    for (i=0; i<4; i++) { iedges[i] = face.iedges[i]; }
    plane  = face.plane;
    for (i=0; i<3; i++) { rmin[i] = face.rmin[i]; rmax[i] = face.rmax[i]; }
    iold   = face.iold;
    inew   = face.inew;
    iprev  = face.iprev;
    inext  = face.inext;
    return *this;
  }

  void invert();
};

// ---------------------------------------------------- Global arrays ---
static vector<ExtNode> nodes;        // vector of nodes
static vector<ExtEdge> edges;        // vector of edges
static vector<ExtFace> faces;        // vector of faces

// ---------------------------------------------------- List of faces ---
class FaceList {
 private:
  int ihead;
  int ilast;

 public:
  FaceList() : ihead(0), ilast(0) {}
  ~FaceList() {}

  void clean() { ihead = 0; ilast = 0; }  
  int front()  { return ihead; }

  void push_back(int i) {
    if (ilast == 0) { ihead = i; } else { faces[ilast].inext = i; } 
    faces[i].iprev = ilast;
    faces[i].inext = 0;
    ilast = i;
  }

  void remove(int i) {
    if (ihead == i) {
      ihead = faces[i].inext;
    }else{
      faces[faces[i].iprev].inext = faces[i].inext;
    }
    if (ilast == i) {
      ilast = faces[i].iprev; 
    }else{
      faces[faces[i].inext].iprev = faces[i].iprev;
    }
    faces[i].iprev = 0;
    faces[i].inext = 0;
  }
};

// --------------------- Polyhedron with extended access to
//                       its members from the BooleanProcessor class ---
class ExtPolyhedron : public HepPolyhedron {
  friend class BooleanProcessor;
  virtual HepPolyhedron& operator = (const HepPolyhedron& from) {
    return HepPolyhedron::operator = (from);
  }
};

// ----------------------------------------- Boolean processor class ---
class BooleanProcessor {
 private:
  int             processor_error;   // is set in case of error
  int             operation;  // 0 (union), 1 (intersection), 2 (subtraction)
  int             ifaces1, ifaces2;  // lists of faces
  int             iout1,   iout2;    // lists of faces with status "out"
  int             iunk1,   iunk2;    // lists of faces with status "unknown"
  double          rmin[3], rmax[3];  // intersection of bounding boxes
  double          del;               // precision (tolerance)   

  FaceList        result_faces;      // list of accepted faces
  FaceList        suitable_faces;    // list of suitable faces
  FaceList        unsuitable_faces;  // list of unsuitable faces
  FaceList        unknown_faces;     // list of unknown faces

  vector<int>     external_contours; // heads of external contours
  vector<int>     internal_contours; // heads of internal contours

 private:
  void   takePolyhedron(const HepPolyhedron & p, double, double, double);
  double findMinMax();
  void   selectOutsideFaces(int & ifaces, int & iout);
  int    testFaceVsPlane(ExtEdge & edge); 
  void   renumberNodes(int & i1, int & i2, int & i3, int & i4);
  int    testEdgeVsEdge(ExtEdge & edge1, ExtEdge & edge2); 
  void   removeJunkNodes() { while(nodes.back().s != 0) nodes.pop_back(); }
  void   divideEdge(int & i1, int & i2);
  void   insertEdge(const ExtEdge & edge);
  void   caseII(ExtEdge & edge1, ExtEdge & edge2);
  void   caseIE(ExtEdge & edge1, ExtEdge & edge2);
  void   caseEE(ExtEdge & edge1, ExtEdge & edge2);
  void   testFaceVsFace(int iface1, int iface2); 
  void   invertNewEdges(int iface);
  void   checkDoubleEdges(int iface);
  void   assembleFace(int what, int iface);
  void   assembleNewFaces(int what, int ihead);
  void   initiateLists();
  void   assemblePolyhedra();
  void   findABC(double x1, double y1, double x2, double y2,
                 double &a, double &b, double &c) const;
  int    checkDirection(double *x, double *y) const;
  int    checkIntersection(int ix, int iy, int i1, int i2) const;
  void   mergeContours(int ix, int iy, int kext, int kint);
  int    checkTriangle(int iedge1, int iedge2, int ix, int iy) const;
  void   triangulateContour(int ix, int iy, int ihead);
  void   modifyReference(int iface, int i1, int i2, int iref);
  void   triangulateFace(int iface);
  HepPolyhedron createPolyhedron();

 public:
  BooleanProcessor() {}
  ~BooleanProcessor() {}

  HepPolyhedron execute(int op,
                        const HepPolyhedron &a,
                        const HepPolyhedron &b);

  void draw();
  void draw_edge(int, int);
  void draw_contour(int, int, int);
  void draw_faces(int, int, int);
  void print_face(int);
  void print_edge(int);
  int get_processor_error() const {return processor_error;}
};

inline void ExtFace::invert()
/***********************************************************************
 *                                                                     *
 * Name: ExtFace::invert()                           Date:    28.02.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Invert face                                               *
 *                                                                     *
 ***********************************************************************/
{
  int iEprev, iEcur, iEnext;

  iEprev = 0; iEcur = iold;
  while (iEcur > 0) {
    edges[iEcur].invert();
    iEnext = edges[iEcur].inext;
    edges[iEcur].inext = iEprev;
    iEprev = iEcur;
    iEcur  = iEnext;
  }
  if (iold > 0) iold = iEprev;

  iEprev = 0; iEcur = inew;
  while (iEcur > 0) {
    edges[iEcur].invert();
    iEnext = edges[iEcur].inext;
    edges[iEcur].inext = iEprev;
    iEprev = iEcur;
    iEcur  = iEnext;
  }
  if (inew > 0) inew = iEprev;

  plane = Plane3D<double>(-plane.a(), -plane.b(), -plane.c(), -plane.d());
}

void BooleanProcessor::takePolyhedron(const HepPolyhedron & p,
                                      double dx, double dy, double dz)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::takePolyhedron            Date:    16.12.99 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Transfer Polyhedron to internal representation            *
 *                                                                     *
 ***********************************************************************/
{
  int i, k, nnode, iNodes[5], iVis[4], iFaces[4];
  int dnode = nodes.size() - 1;
  int dface = faces.size() - 1;

  //   S E T   N O D E S

  //  for (i=1; i <= p.GetNoVertices(); i++) {  
  //    nodes.push_back(ExtNode(p.GetVertex(i)));
  //  }

  Point3D<double> ppp;        
  for (i=1; i <= p.GetNoVertices(); i++) {  
    ppp = p.GetVertex(i);
    ppp.setX(ppp.x()+dx);
    ppp.setY(ppp.y()+dy);
    ppp.setZ(ppp.z()+dz);
    nodes.push_back(ExtNode(ppp));
  }

  //   S E T   F A C E S

  for (int iface=1; iface <= p.GetNoFacets(); iface++) {
    faces.push_back(ExtFace(edges.size()));

    //   S E T   F A C E   N O D E S

    p.GetFacet(iface, nnode, iNodes, iVis, iFaces);
    for (i=0; i<nnode; i++) {
      if (iNodes[i] < 1 || iNodes[i] > p.GetNoVertices()) processor_error = 1;
      if (iFaces[i] < 1 || iFaces[i] > p.GetNoFacets())   processor_error = 1;
      iNodes[i] += dnode;
      iFaces[i] += dface;
    }

    //   S E T   E D G E S

    iNodes[nnode] = iNodes[0];
    faces.back().iedges[3] = 0;
    for (i=0; i<nnode; i++) {
      faces.back().iedges[i] = edges.size();
      edges.push_back(ExtEdge(iNodes[i], iNodes[i+1],
                              iface+dface, iFaces[i], iVis[i]));
      edges.back().inext     = edges.size();
    }
    edges.back().inext = 0;

    //   S E T   F A C E   M I N - M A X

    for (i=0; i<3; i++) {
      faces.back().rmin[i] = nodes[iNodes[0]].v[i];
      faces.back().rmax[i] = nodes[iNodes[0]].v[i];
    }
    for (i=1; i<nnode; i++) {
      for (k=0; k<3; k++) {
        if (faces.back().rmin[k] > nodes[iNodes[i]].v[k])
            faces.back().rmin[k] = nodes[iNodes[i]].v[k];
        if (faces.back().rmax[k] < nodes[iNodes[i]].v[k])
            faces.back().rmax[k] = nodes[iNodes[i]].v[k];
      }
    }

    //   S E T   F A C E   P L A N E

    Normal3D<double> n = (nodes[iNodes[2]].v-nodes[iNodes[0]].v).cross
                    (nodes[iNodes[3]].v-nodes[iNodes[1]].v);
    Point3D<double>  p(0,0,0);
    
    for (i=0; i<nnode; i++) { p += nodes[iNodes[i]].v; }
    p *= (1./nnode);
    faces.back().plane = Plane3D<double>(n.unit(), p);

    //   S E T   R E F E R E N C E   T O   T H E   N E X T   F A C E

    faces.back().inext = faces.size(); 
  }
  faces.back().inext = 0; 
}

double BooleanProcessor::findMinMax()
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::findMinMax                Date:    16.12.99 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Find min-max (bounding) boxes for polyhedra               *
 *                                                                     *
 ***********************************************************************/
{
  if (ifaces1 == 0 || ifaces2 == 0) return 0;

  int    i, iface;
  double rmin1[3], rmax1[3];
  double rmin2[3], rmax2[3];

  //   F I N D   B O U N D I N G   B O X E S

  for (i=0; i<3; i++) {
    rmin1[i] = faces[ifaces1].rmin[i];
    rmax1[i] = faces[ifaces1].rmax[i];
    rmin2[i] = faces[ifaces2].rmin[i];
    rmax2[i] = faces[ifaces2].rmax[i];
  }

  iface = faces[ifaces1].inext;
  while(iface > 0) {
    for (i=0; i<3; i++) {
      if (rmin1[i] > faces[iface].rmin[i]) rmin1[i] = faces[iface].rmin[i];
      if (rmax1[i] < faces[iface].rmax[i]) rmax1[i] = faces[iface].rmax[i];
    }
    iface = faces[iface].inext;
  }

  iface = faces[ifaces2].inext;
  while(iface > 0) {
    for (i=0; i<3; i++) {
      if (rmin2[i] > faces[iface].rmin[i]) rmin2[i] = faces[iface].rmin[i];
      if (rmax2[i] < faces[iface].rmax[i]) rmax2[i] = faces[iface].rmax[i];
    }
    iface = faces[iface].inext;
  }

  //   F I N D   I N T E R S E C T I O N   O F   B O U N D I N G   B O X E S

  for (i=0; i<3; i++) {
    rmin[i] = (rmin1[i] > rmin2[i]) ? rmin1[i] : rmin2[i]; 
    rmax[i] = (rmax1[i] < rmax2[i]) ? rmax1[i] : rmax2[i]; 
  }

  //   F I N D   T O L E R A N C E

  double del1 = 0;
  double del2 = 0;
  for (i=0; i<3; i++) {
    if ((rmax1[i]-rmin1[i]) > del1) del1 = rmax1[i]-rmin1[i]; 
    if ((rmax2[i]-rmin2[i]) > del2) del2 = rmax2[i]-rmin2[i]; 
  }
  return ((del1 < del2) ? del1 : del2) / GRANULARITY;
}

void BooleanProcessor::selectOutsideFaces(int & ifaces, int & iout)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::selectOutsideFaces        Date:    10.01.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Preselection of outside faces                             *
 *                                                                     *
 ***********************************************************************/
{
  int i, outflag, iface = ifaces, *prev;
  Point3D<double> mmbox[8] = { Point3D<double>(rmin[0],rmin[1],rmin[2]),
                               Point3D<double>(rmax[0],rmin[1],rmin[2]),
                               Point3D<double>(rmin[0],rmax[1],rmin[2]),
                               Point3D<double>(rmax[0],rmax[1],rmin[2]),
                               Point3D<double>(rmin[0],rmin[1],rmax[2]),
                               Point3D<double>(rmax[0],rmin[1],rmax[2]),
                               Point3D<double>(rmin[0],rmax[1],rmax[2]),
                               Point3D<double>(rmax[0],rmax[1],rmax[2]) };
  prev = &ifaces;
  while (iface > 0) {

    //   B O U N D I N G   B O X   vs  B O U N D I N G   B O X 

    outflag = 0;
    for (i=0; i<3; i++) {
      if (faces[iface].rmin[i] > rmax[i] + del) { outflag = 1; break; }
      if (faces[iface].rmax[i] < rmin[i] - del) { outflag = 1; break; }
    }

    //   B O U N D I N G   B O X   vs  P L A N E

    if (outflag == 0) {
      int npos = 0, nneg = 0;
      double d;
      for (i=0; i<8; i++) {
        d = faces[iface].plane.distance(mmbox[i]);
        if (d > +del) npos++;
        if (d < -del) nneg++;
      }
      if (npos == 8 || nneg == 8) outflag = 1;
    }

    //   U P D A T E   L I S T S

    if (outflag == 1) {
      *prev = faces[iface].inext;
      faces[iface].inext = iout;
      iout = iface;
    }else{
      prev = &faces[iface].inext;
    }
    iface = *prev;
  }
}

int BooleanProcessor::testFaceVsPlane(ExtEdge & edge)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::testFaceVsPlane           Date:    19.01.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Find result of intersection of face by plane              *
 *                                                                     *
 ***********************************************************************/
{
  int        iface = edge.iface1;
  Plane3D<double> plane = faces[edge.iface2].plane;
  int        i, nnode, npos = 0, nneg = 0, nzer = 0;
  double     dd[5];

  //   F I N D   D I S T A N C E S

  nnode = (faces[iface].iedges[3] == 0) ? 3 : 4;
  for (i=0; i<nnode; i++) {
    dd[i] = plane.distance(nodes[edges[faces[iface].iedges[i]].i1].v);
    if (dd[i] > del) {
      npos++; 
    }else if (dd[i] < -del) {
      nneg++;
    }else{
      nzer++; dd[i] = 0;
    }  
  }

  //   S O M E   S I M P L E   C A S E S  ( N O   I N T E R S E C T I O N )

  if (npos == nnode || nneg == nnode)   return OUT_OF_PLANE;
  if (nzer == 1     && nneg == 0)       return OUT_OF_PLANE;
  if (nzer == 1     && npos == 0)       return OUT_OF_PLANE;
  if (nzer == nnode)                    return ON_PLANE;
  if (nzer == 3)                        return NON_PLANAR_FACE;

  //   F I N D   I N T E R S E C T I O N

  int       ie1 = 0, ie2 = 0, s1 = 0, s2 = 0, status, nint = 0; 
  enum      { PLUS_MINUS, MINUS_PLUS, ZERO_ZERO, ZERO_PLUS, ZERO_MINUS };

  dd[nnode] = dd[0];
  for (i=0; i<nnode; i++) {
    if (dd[i] > 0) {
      if (dd[i+1] >= 0) continue;
      status = PLUS_MINUS;
    }else if (dd[i] < 0) {
      if (dd[i+1] <= 0) continue;
      status = MINUS_PLUS;
    }else{   
      status = ZERO_ZERO;
      if (dd[i+1] > 0) status = ZERO_PLUS;
      if (dd[i+1] < 0) status = ZERO_MINUS;
    }
    switch (nint) {
    case 0:
      ie1 = i; s1 = status; nint++; break;
    case 1:
      ie2 = i; s2 = status; nint++; break;
    default:
      return NON_PLANAR_FACE;
    }
  }
  if (nint != 2)                        return NON_PLANAR_FACE;

  //   F O R M   I N T E R S E C T I O N   S E G M E N T

  if (s1 != ZERO_ZERO && s2 != ZERO_ZERO) {
    if (s1 == s2)                       return NON_PLANAR_FACE;
    int     iedge, i1 = 0, i2 = 0, ii[2];
    double  d1 = 0., d2 = 0., dd ;
    ii[0] = ie1; ii[1] = ie2;
    for (i=0; i<2; i++) {
      iedge = faces[iface].iedges[ii[i]];
      while (iedge > 0) {
        i1 = edges[iedge].i1;
        i2 = edges[iedge].i2;
        d1 = plane.distance(nodes[i1].v);
        d2 = plane.distance(nodes[i2].v);
        if (d1 > del) {
          if (d2 < -del) { ii[i] = nodes.size(); break; } // +-
        }else if (d1 < -del) {
          if (d2 >  del) { ii[i] = nodes.size(); break; } // -+
        }else{ 
          ii[i] = i1; break;                              // 0+ or 0-
        }
        iedge = edges[iedge].inext;
      }
      if (ii[i] == nodes.size()) {
        dd = d2-d1; d1 = d1/dd; d2 = d2/dd;
        nodes.push_back(ExtNode(d2*nodes[i1].v-d1*nodes[i2].v, iedge));
      }  
    }
    edge.inext = 0;
    if (s1 == MINUS_PLUS || s1 == ZERO_PLUS) {
      edge.i1 = ii[1];
      edge.i2 = ii[0];
    }else{
      edge.i1 = ii[0];
      edge.i2 = ii[1];
    }
    return INTERSECTION;
  }else{
    if (npos == nneg)                   return NON_PLANAR_FACE;
    edge.inext = (s1 == ZERO_ZERO) ? ie1+1 : ie2+1;
    if (s1 == ZERO_PLUS || s2 == ZERO_MINUS) {
      edge.i1 = edges[faces[iface].iedges[ie2]].i1;
      edge.i2 = edges[faces[iface].iedges[ie1]].i1;
    }else{
      edge.i1 = edges[faces[iface].iedges[ie1]].i1;
      edge.i2 = edges[faces[iface].iedges[ie2]].i1;
    }
    return EDGE;
  }
}

void BooleanProcessor::renumberNodes(int & i1, int & i2, int & i3, int & i4)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::renumberNodes             Date:    19.01.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Renumber nodes and remove last temporary node.            *
 *           Remark: In principal this routine can be replaced just    *
 *           with i1 = i2;                                             *
 *           Removal of temporary nodes provides additional control    *
 *           on number of nodes, that is very useful for debugging.    *
 *                                                                     *
 ***********************************************************************/
{
  if (i1 == i2) return;
  if (nodes[i1].s == 0 || nodes.back().s == 0) { i1 = i2; return; }

  int ilast = nodes.size()-1;
  if (i1 == ilast) { i1 = i2; nodes.pop_back(); return; }
  if (i2 == ilast) { i2 = i1; }
  if (i3 == ilast) { i3 = i1; }
  if (i4 == ilast) { i4 = i1; }
  nodes[i1] = nodes.back(); i1 = i2; nodes.pop_back();
}

int BooleanProcessor::testEdgeVsEdge(ExtEdge & edge1, ExtEdge & edge2)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::testEdgeVsEdge            Date:    19.01.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Find common part of two edges                             *
 *                                                                     *
 ***********************************************************************/
{
  int    i, ii = 0;
  double d, dd = 0.;

  for (i=0; i<3; i++) {
    d = nodes[edge1.i1].v[i]-nodes[edge1.i2].v[i];
    if (d < 0.) d = -d; 
    if (d > dd) { dd = d; ii = i; }
  }
  double t1 = nodes[edge1.i1].v[ii];
  double t2 = nodes[edge1.i2].v[ii];
  double t3 = nodes[edge2.i1].v[ii];
  double t4 = nodes[edge2.i2].v[ii];
  if (t2-t1 < 0.) { t1 = -t1; t2 = -t2; t3 = -t3; t4 = -t4; }
 
  if (t3 <= t1+del || t4 >= t2-del) return 0;
  if (t3 > t2+del) {
    renumberNodes(edge2.i1, edge1.i2, edge1.i1, edge2.i2);
  }else if (t3 < t2-del) {
    renumberNodes(edge1.i2, edge2.i1, edge1.i1, edge2.i2);
  }
  if (t4 < t1-del) {
    renumberNodes(edge2.i2, edge1.i1, edge1.i2, edge2.i1);
  }else if (t4 > t1+del) {
    renumberNodes(edge1.i1, edge2.i2, edge1.i2, edge2.i1);
  }
  return 1;
}

void BooleanProcessor::divideEdge(int & i1, int & i2)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::divideEdge                Date:    24.01.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Unify the nodes and divide edge on two parts by the node. *
 *                                                                     *
 ***********************************************************************/
{
  int iedges[2];
  iedges[0] = nodes[i1].s;
  iedges[1] = nodes[i2].s;

  //   U N I F Y   N O D E S
  
  if      (i1 < i2) { i2 = i1; }
  else if (i1 > i2) { i1 = i2; }
  else              { iedges[1] = 0; }  
  if (iedges[0] == iedges[1]) return;

  int ie1, ie2, inode = i1;
  nodes[inode].s = 0;
  for (int i=0; i<2; i++) {

    //   F I N D   C O R R E S P O N D I N G   E D G E

    if ((ie1 = iedges[i]) == 0) continue;
    ie2 = faces[edges[ie1].iface2].iedges[0];
    while (ie2 > 0) {
      if (edges[ie2].i1 == edges[ie1].i2 &&
          edges[ie2].i2 == edges[ie1].i1) break;
      ie2 = edges[ie2].inext;
    }

    //   D I V I D E   E D G E S
    
    edges.push_back(edges[ie1]);
    edges[ie1].inext = edges.size() - 1;
    edges[ie1].i2    = inode;
    edges.back().i1  = inode;
    
    edges.push_back(edges[ie2]);
    edges[ie2].inext = edges.size() - 1;
    edges[ie2].i2    = inode;
    edges.back().i1  = inode;
  } 
}

void BooleanProcessor::insertEdge(const ExtEdge & edge)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::insertEdge                Date:    24.01.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Insert edge to the list of new edges                      *
 *                                                                     *
 ***********************************************************************/
{
  int iface = edge.iface1;
  edges.push_back(edge);
  edges.back().inext = faces[iface].inew;
  faces[iface].inew  = edges.size() - 1;
}

void BooleanProcessor::caseII(ExtEdge & edge1, ExtEdge & edge2)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::caseII                    Date:    19.01.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Intersection/Intersection case                            *
 *                                                                     *
 ***********************************************************************/
{
  divideEdge(edge1.i1, edge2.i2);
  divideEdge(edge1.i2, edge2.i1);
  edge1.ivis = 1;
  edge2.ivis = 1;
  insertEdge(edge1);
  insertEdge(edge2);
}

void BooleanProcessor::caseIE(ExtEdge &, ExtEdge &)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::caseIE                    Date:    19.01.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Intersection/Edge-touch case                              *
 *                                                                     *
 ***********************************************************************/
{
  processor_error = 1;
  std::cout
    << "BooleanProcessor::caseIE : unimplemented case"
    << std::endl;
}

void BooleanProcessor::caseEE(ExtEdge &, ExtEdge &)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::caseEE                    Date:    19.01.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Edge-touch/Edge-touch case                                *
 *                                                                     *
 ***********************************************************************/
{
  processor_error = 1;
  std::cout
    << "BooleanProcessor::caseEE : unimplemented case"
    << std::endl;
}

void BooleanProcessor::testFaceVsFace(int iface1, int iface2)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::testFaceVsFace            Date:    11.01.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Find result (an edge) of intersection of two faces        *
 *                                                                     *
 ***********************************************************************/
{
  ExtEdge edge1, edge2;
  int     irep1, irep2;

  //   M I N - M A X

  for (int i=0; i<3; i++) {
    if (faces[iface1].rmin[i] > faces[iface2].rmax[i] + del) return;
    if (faces[iface1].rmax[i] < faces[iface2].rmin[i] - del) return;
  }

  //   F A C E - 1   vs   P L A N E - 2

  edge1.iface1 = iface1;
  edge1.iface2 = iface2;
  irep1        = testFaceVsPlane(edge1); 
  if (irep1 == OUT_OF_PLANE || irep1 == ON_PLANE) {
    removeJunkNodes();
    return;
  }

  //   F A C E - 2   vs   P L A N E - 1

  edge2.iface1 = iface2;
  edge2.iface2 = iface1;
  irep2        = testFaceVsPlane(edge2); 
  if (irep2 == OUT_OF_PLANE || irep2 == ON_PLANE) {
    removeJunkNodes();
    return;
  }

  //   C H E C K   F O R   N O N P L A N A R   F A C E 

  if (irep1 == NON_PLANAR_FACE || irep2 == NON_PLANAR_FACE) {
    removeJunkNodes();
    return;
  }

  //   F I N D   I N T E R S E C T I O N   P A R T

  if (testEdgeVsEdge(edge1, edge2) == 0) return;

  //   C O N S I D E R   D I F F E R E N T   C A S E S  

  if (irep1 == INTERSECTION && irep2 == INTERSECTION) caseII(edge1, edge2);
  if (irep1 == INTERSECTION && irep2 == EDGE)         caseIE(edge1, edge2);
  if (irep1 == EDGE         && irep2 == INTERSECTION) caseIE(edge2, edge1);
  if (irep1 == EDGE         && irep2 == EDGE)         caseEE(edge1, edge2);
  removeJunkNodes();

}

void BooleanProcessor::invertNewEdges(int iface)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::invertNewEdges            Date:    04.02.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Invert direction of new edges                             *
 *                                                                     *
 ***********************************************************************/
{
  int iedge = faces[iface].inew;
  while (iedge > 0) {
    edges[iedge].invert();
    iedge = edges[iedge].inext;
  }
}

void BooleanProcessor::checkDoubleEdges(int)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::checkDoubleEdges          Date:    04.02.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Eliminate duplication of edges                            *
 *                                                                     *
 ***********************************************************************/
{

}

void BooleanProcessor::assembleFace(int what, int iface)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::assembleFace              Date:    19.02.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Assemble face                                             *
 *                                                                     *
 ***********************************************************************/
{
  //   A S S E M B L E   N E W   F A C E 

  int ihead;      // head of the list of edges for new face
  int icur;       // current edge in the list - last edge inserted to the list
  int *ilink;     // pointer to the current link 
  int ifirst;     // first node of a contour
  int *i;         // pointer to the index of the current edge in a loop
  int ioldflag=0; // is set if an edge from iold has been taken

#define INSERT_EDGE_TO_THE_LIST(A) \
*ilink = A; ilink = &edges[A].inext; *ilink = 0

  ilink = &ihead;
  for(;;) {
    if (faces[iface].inew == 0) break;

    //   S T A R T   N E W   C O N T O U R

    icur   = faces[iface].inew;
    faces[iface].inew = edges[icur].inext;
    INSERT_EDGE_TO_THE_LIST(icur);
    ifirst = edges[icur].i1; 
    
    //   C O N S T R U C T   T H E   C O N T O U R

    for (;;) {
      i = &faces[iface].inew;
      while(*i > 0) {
        if (edges[*i].i1 == edges[icur].i2) break;
        i = &edges[*i].inext;
      }
      if (*i == 0) {
        i = &faces[iface].iold;
        while(*i > 0) {
          if (edges[*i].i1 == edges[icur].i2) ioldflag = 1;
          if (edges[*i].i1 == edges[icur].i2) break;
          i = &edges[*i].inext;
        }
      }
      if (*i > 0) {
        icur = *i;
        *i = edges[icur].inext;
        INSERT_EDGE_TO_THE_LIST(icur);
        if (edges[icur].i2 == ifirst) { break; } else { continue; }
      }else{
        processor_error = 1;
        std::cerr
          << "BooleanProcessor::assembleFace(" << iface << ") : "
          << "could not find next edge of the contour"
          << std::endl;
        faces[iface].inew = DEFECTIVE_FACE;
        return;
      }
    }
  }

  //   C H E C K   O R I G I N A L   C O N T O U R

  int iedge;
  iedge = faces[iface].iold;
  if (what == 0 && ioldflag == 0 && iedge > 0) {
    for (;;) {
      if (edges[iedge].inext > 0) {
        if (edges[iedge].i2 == edges[edges[iedge].inext].i1) {
          iedge = edges[iedge].inext;
        }else{
          break;
        }
      }else{
        if (edges[iedge].i2 == edges[faces[iface].iold].i1) {
          edges[iedge].inext = ihead;   // set new face
          return;
        }else{
          break;
        }
      }
    }
  }

  //   M A R K   U N S U I T A B L E   N E I G H B O U R I N G   F A C E S

  int iface2;
  iedge = faces[iface].iold;
  while(iedge > 0) {
    iface2 = edges[iedge].iface2;
    if (faces[iface2].inew == 0) faces[iface2].inew = UNSUITABLE_FACE;
    iedge = edges[iedge].inext;
  }
  faces[iface].iold = ihead;            // set new face
}

void BooleanProcessor::assembleNewFaces(int what, int ihead)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::assembleNewFaces          Date:    30.01.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Assemble internal or external parts of faces              *
 *                                                                     *
 ***********************************************************************/
{
  int iface = ihead;
  while(iface > 0) {
    if (faces[iface].inew > 0) {
      if (what != 0) invertNewEdges(iface);
      checkDoubleEdges(iface);
      assembleFace(what, iface);
      faces[iface].inew =
        (faces[iface].iold == 0) ? UNSUITABLE_FACE : NEW_FACE;
    }
    iface = faces[iface].inext;
  }
}

void BooleanProcessor::initiateLists()
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::initiateLists             Date:    28.02.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Initiate lists of faces.                                  *
 *                                                                     *
 ***********************************************************************/
{
  int i, iface;

  //   R E S E T   L I S T S   O F   F A C E S

  result_faces.clean();
  suitable_faces.clean();
  unsuitable_faces.clean();
  unknown_faces.clean();

  //   I N I T I A T E   T H E   L I S T S

  iface = iout1;
  while (iface > 0) {
    i     = iface;
    iface = faces[i].inext;
    if (operation == OP_INTERSECTION) {
      unsuitable_faces.push_back(i);
      faces[i].inew = UNSUITABLE_FACE;
    }else{
      suitable_faces.push_back(i);
      faces[i].inew = ORIGINAL_FACE;
    }
  }
  iface = iout2;
  while (iface > 0) {
    i     = iface;
    iface = faces[i].inext;
    if (operation == OP_UNION) {
      suitable_faces.push_back(i);
      faces[i].inew = ORIGINAL_FACE;
    }else{
      unsuitable_faces.push_back(i);
      faces[i].inew = UNSUITABLE_FACE;
    }
  }

  iface = iunk1;
  while (iface > 0) {
    i     = iface;
    iface = faces[i].inext;
    unknown_faces.push_back(i);
  }
  iface = iunk2;
  while (iface > 0) {
    i     = iface;
    iface = faces[i].inext;
    if (operation == OP_SUBTRACTION) faces[i].invert();
    unknown_faces.push_back(i);
  }

  iface = ifaces1;
  while (iface > 0) {
    i     = iface;
    iface = faces[i].inext;
    switch(faces[i].inew) {
    case UNKNOWN_FACE:
      unknown_faces.push_back(i);
      break;
    case ORIGINAL_FACE: case NEW_FACE:
      suitable_faces.push_back(i);
      break;
    case UNSUITABLE_FACE:
      unsuitable_faces.push_back(i);
      break;
    default:
      faces[i].iprev = 0;
      faces[i].inext = 0;
      break;
    }
  }
  iface = ifaces2;
  while (iface > 0) {
    i     = iface;
    iface = faces[i].inext;
    if (operation == OP_SUBTRACTION) faces[i].invert();
    switch(faces[i].inew) {
    case UNKNOWN_FACE:
      unknown_faces.push_back(i);
      break;
    case ORIGINAL_FACE: case NEW_FACE:
      suitable_faces.push_back(i);
      break;
    case UNSUITABLE_FACE:
      unsuitable_faces.push_back(i);
      break;
    default:
      faces[i].iprev = 0;
      faces[i].inext = 0;
      break;
    }
  }
  ifaces1 = ifaces2 = iout1 = iout2 = iunk1 = iunk2 = 0;
}

void BooleanProcessor::assemblePolyhedra()
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::assemblePolyhedra()       Date:    10.12.99 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Collect suitable faces and remove unsuitable ones.        *
 *                                                                     *
 ***********************************************************************/
{
  int i, iedge, iface;

  //   L O O P   A L O N G   S U I T A B L E   F A C E S

  iface = suitable_faces.front();
  while(iface > 0) {
    i = iface;
    iedge = faces[i].iold;
    while(iedge > 0) {
      iface = edges[iedge].iface2;
      if (faces[iface].inew == UNKNOWN_FACE) {
        unknown_faces.remove(iface);
        suitable_faces.push_back(iface);
        faces[iface].inew = ORIGINAL_FACE;
      }
      iedge = edges[iedge].inext;
    }
    iface = faces[i].inext;
    suitable_faces.remove(i);
    result_faces.push_back(i);
  }
  if (unknown_faces.front() == 0) return;

  //   L O O P   A L O N G   U N S U I T A B L E   F A C E S

  iface = unsuitable_faces.front();
  while(iface > 0) {
    i = iface;
    iedge = faces[i].iold;
    while(iedge > 0) {
      iface = edges[iedge].iface2;
      if (faces[iface].inew == UNKNOWN_FACE) {
        unknown_faces.remove(iface);
        unsuitable_faces.push_back(iface);
        faces[iface].inew = UNSUITABLE_FACE;
      }
      iedge = edges[iedge].inext;
    }
    iface = faces[i].inext;
    unsuitable_faces.remove(i);
  }
}

inline void
BooleanProcessor::findABC(double x1, double y1, double x2, double y2,
                          double &a, double &b, double &c) const
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::findABC                   Date:    07.03.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Find line equation Ax+By+C=0                              *
 *                                                                     *
 ***********************************************************************/
{
  double w;
  a  = y1 - y2;
  b  = x2 - x1;
  w  = std::abs(a)+std::abs(b);
  a /= w;
  b /= w;
  c  = -(a*x2 + b*y2);
}

int BooleanProcessor::checkDirection(double *x, double *y) const
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::checkDirection            Date:    06.03.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Check direction of line 1-4                               *
 *                                                                     *
 ***********************************************************************/
{
  double a1, b1, c1, a2, b2, c2, d1, d2;

  //   T E S T   L I N E   1 - 4   V S   E X T E R N A L   C O N T O U R

  findABC(x[0], y[0], x[1], y[1], a1, b1, c1);
  findABC(x[1], y[1], x[2], y[2], a2, b2, c2);
  d1 = a1*x[4] + b1*y[4] + c1;
  d2 = a2*x[4] + b2*y[4] + c2;
  if (d1 <= del && d2 <= del)            return 1;
  if (! (d1 > del && d2 > del)) {
    if ( a1*x[2] + b1*y[2] + c1 >= -del) return 1;
  }

  //   T E S T   L I N E   1 - 4   V S   I N T E R N A L   C O N T O U R
  
  findABC(x[3], y[3], x[4], y[4], a1, b1, c1);
  findABC(x[4], y[4], x[5], y[5], a2, b2, c2);
  d1 = a1*x[1] + b1*y[1] + c1;
  d2 = a2*x[1] + b2*y[1] + c2;
  if (d1 <= del && d2 <= del)            return 1;
  if (!(d1 > del && d2 > del)) {
    if ( a1*x[5] + b1*y[5] + c1 >= -del) return 1;
  }
  return 0;
}

int BooleanProcessor::checkIntersection(int ix, int iy, int i1, int i2) const
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::checkDirection            Date:    06.03.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Check line i1-i2 on intersection with contours            *
 *                                                                     *
 ***********************************************************************/
{
  //  F I N D   L I N E   E Q U A T I O N

  double x1, y1, x2, y2, a1, b1, c1;
  x1 = nodes[i1].v[ix];
  y1 = nodes[i1].v[iy];
  x2 = nodes[i2].v[ix];
  y2 = nodes[i2].v[iy];
  findABC(x1, y1, x2, y2, a1, b1, c1);

  //  L O O P   A L O N G   E X T E R N A L   C O N T O U R S

  int icontour, iedge, k1, k2;
  double x3, y3, x4, y4, a2, b2, c2, d1, d2;
  for(icontour=0; icontour<external_contours.size(); icontour++) {
    iedge = external_contours[icontour];
    while(iedge > 0) {
      k1 = edges[iedge].i1;
      k2 = edges[iedge].i2;
      iedge = edges[iedge].inext;
      if (k1 == i1 || k2 == i1) continue;
      if (k1 == i2 || k2 == i2) continue;
      x3 = nodes[k1].v[ix];
      y3 = nodes[k1].v[iy];
      x4 = nodes[k2].v[ix];
      y4 = nodes[k2].v[iy];
      d1 = a1*x3 + b1*y3 + c1;
      d2 = a1*x4 + b1*y4 + c1;
      if (d1 >  del && d2 >  del) continue;
      if (d1 < -del && d2 < -del) continue;

      findABC(x3, y3, x4, y4, a2, b2, c2);
      d1 = a2*x1 + b2*y1 + c2;
      d2 = a2*x2 + b2*y2 + c2;
      if (d1 >  del && d2 >  del) continue;
      if (d1 < -del && d2 < -del) continue;
      return 1;
    }
  }

  //  L O O P   A L O N G   E X T E R N A L   C O N T O U R S

  for(icontour=0; icontour<internal_contours.size(); icontour++) {
    iedge = internal_contours[icontour];
    while(iedge > 0) {
      k1 = edges[iedge].i1;
      k2 = edges[iedge].i2;
      iedge = edges[iedge].inext;
      if (k1 == i1 || k2 == i1) continue;
      if (k1 == i2 || k2 == i2) continue;
      x3 = nodes[k1].v[ix];
      y3 = nodes[k1].v[iy];
      x4 = nodes[k2].v[ix];
      y4 = nodes[k2].v[iy];
      d1 = a1*x3 + b1*y3 + c1;
      d2 = a1*x4 + b1*y4 + c1;
      if (d1 >  del && d2 >  del) continue;
      if (d1 < -del && d2 < -del) continue;

      findABC(x3, y3, x4, y4, a2, b2, c2);
      d1 = a2*x1 + b2*y1 + c2;
      d2 = a2*x2 + b2*y2 + c2;
      if (d1 >  del && d2 >  del) continue;
      if (d1 < -del && d2 < -del) continue;
      return 1;
    }
  }
  return 0;
}

void BooleanProcessor::mergeContours(int ix, int iy, int kext, int kint)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::mergeContours             Date:    06.03.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Attemp to merge internal contour with external one        *
 *                                                                     *
 ***********************************************************************/
{
  int    i1ext, i2ext, i1int, i2int, i, k[6];
  double x[6], y[6];

  //   L O O P   A L O N G   E X T E R N A L   C O N T O U R

  i1ext = external_contours[kext];
  while (i1ext > 0) {
    i2ext = edges[i1ext].inext;
    if (i2ext == 0) i2ext = external_contours[kext];
    k[0] = edges[i1ext].i1;
    k[1] = edges[i1ext].i2;
    k[2] = edges[i2ext].i2; 
    for (i=0; i<3; i++) {
      x[i] = nodes[k[i]].v[ix];
      y[i] = nodes[k[i]].v[iy];
    }

    //   L O O P   A L O N G   I N T E R N A L   C O N T O U R

    i1int = internal_contours[kint];
    while (i1int > 0) {
      i2int = edges[i1int].inext;
      if (i2int == 0) i2int = internal_contours[kint];
      k[3] = edges[i1int].i1;
      k[4] = edges[i1int].i2;
      k[5] = edges[i2int].i2; 
      for (i=3; i<6; i++) {
        x[i] = nodes[k[i]].v[ix];
        y[i] = nodes[k[i]].v[iy];
      }

      //   T E S T   L I N E   K1 - K4
      //   I F   O K   T H E N   M E R G E   C O N T O U R S

      if (checkDirection(x, y) == 0) {
        if (checkIntersection(ix, iy, k[1], k[4]) == 0) {
          i = i1int;
          for(;;) {
            if (edges[i].inext == 0) {
              edges[i].inext = internal_contours[kint];
              internal_contours[kint] = 0;
              break;
            }else{
              i = edges[i].inext;
            }
          }
          i = edges[i1int].iface1;
          edges.push_back(ExtEdge(k[1], k[4], i, -(edges.size()+1), -1));
          edges.back().inext = i2int;
          edges.push_back(ExtEdge(k[4], k[1], i, -(edges.size()-1), -1));
          edges.back().inext = edges[i1ext].inext;
          edges[i1ext].inext = edges.size()-2;
          edges[i1int].inext = edges.size()-1;
          return;
        }
      }
      i1int = edges[i1int].inext;
    }
    i1ext = edges[i1ext].inext;
  }
}

int
BooleanProcessor::checkTriangle(int iedge1, int iedge2, int ix, int iy) const
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::checkTriangle             Date:    08.03.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Check triangle for correctness                            *
 *                                                                     *
 ***********************************************************************/
{
  int    k[3];
  double x[3], y[3];
  double a1, b1, c1;

  k[0] = edges[iedge1].i1;
  k[1] = edges[iedge1].i2;
  k[2] = edges[iedge2].i2;
  for (int i=0; i<3; i++) {
    x[i] = nodes[k[i]].v[ix];
    y[i] = nodes[k[i]].v[iy];
  }

  //  C H E C K   P R I N C I P A L   C O R R E C T N E S S  
    
  findABC(x[2], y[2], x[0], y[0], a1, b1, c1);
  if (a1*x[1]+b1*y[1]+c1 <= 0.1*del) return 1;

  //   C H E C K   T H A T   T H E R E   I S   N O   P O I N T S   I N S I D E

  int    inode, iedge;
  double a2, b2, c2, a3, b3, c3;

  findABC(x[0], y[0], x[1], y[1], a2, b2, c2);
  findABC(x[1], y[1], x[2], y[2], a3, b3, c3);
  iedge = iedge2;
  for (;;) {
    iedge = edges[iedge].inext;
    if (edges[iedge].inext == iedge1) return 0;
    inode = edges[iedge].i2;
    if (inode == k[0])                continue;
    if (inode == k[1])                continue;
    if (inode == k[2])                continue;
    x[1]  = nodes[inode].v[ix];
    y[1]  = nodes[inode].v[iy];
    if (a1*x[1]+b1*y[1]+c1 < -0.1*del)    continue;
    if (a2*x[1]+b2*y[1]+c2 < -0.1*del)    continue;
    if (a3*x[1]+b3*y[1]+c3 < -0.1*del)    continue;
    return 1;
  }
}

void BooleanProcessor::triangulateContour(int ix, int iy, int ihead)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::triangulateContour        Date:    06.03.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Triangulate external contour                              *
 *                                                                     *
 ***********************************************************************/
{
   
  //std::cout << "Next Countour" << std::endl;
  //int draw_flag = 0;
  //if (draw_flag) draw_contour(5, 3, ihead);
    
  //   C L O S E   C O N T O U R
  
  int ipnext = ihead, nnode = 1;
  for (;;) {
    if (edges[ipnext].inext > 0) {
      ipnext = edges[ipnext].inext;
      nnode++; 
    }else{
      edges[ipnext].inext = ihead;
      break;
    }
  }

  //   L O O P   A L O N G   C O N T O U R  

  int iedge1, iedge2, iedge3, istart = 0;
  for (;;) {
    iedge1 = edges[ipnext].inext;
    iedge2 = edges[iedge1].inext;
    if (istart == 0) {
      istart = iedge1;
      if (nnode <= 3) {
        iedge3 = edges[iedge2].inext;
        edges[iedge1].iface1 = faces.size();
        edges[iedge2].iface1 = faces.size();
        edges[iedge3].iface1 = faces.size();
        edges[iedge3].inext = 0;
        faces.push_back(ExtFace());
        faces.back().iold = iedge1;
        faces.back().inew = ORIGINAL_FACE;

  //if (draw_flag) draw_contour(4, 2, iedge1);

        break;
      }
    }else if (istart == iedge1) {
      processor_error = 1;
      std::cerr
        << "BooleanProcessor::triangulateContour : "
        << "could not generate a triangle (infinite loop)"
        << std::endl;
      break;
    }

    //   C H E C K   C O R E C T N E S S   O F   T H E   T R I A N G L E

    if (checkTriangle(iedge1,iedge2,ix,iy) != 0) {
      ipnext  = edges[ipnext].inext;
      continue; 
    }

    //   M O D I F Y   C O N T O U R  
    
    int i1 = edges[iedge1].i1;
    int i3 = edges[iedge2].i2;
    int iface1 = edges[iedge1].iface1;
    int iface2 = faces.size();

    edges[ipnext].inext = edges.size();
    edges.push_back(ExtEdge(i1, i3, iface1, -(edges.size()+1), -1));
    edges.back().inext = edges[iedge2].inext;

    //   A D D   N E W   T R I A N G L E   T O   T H E   L I S T

    edges[iedge2].inext = edges.size();
    edges.push_back(ExtEdge(i3, i1, iface2, -(edges.size()-1), -1));
    faces.push_back(ExtFace());
    faces.back().iold   = iedge1;
    faces.back().inew   = ORIGINAL_FACE;
    edges[iedge1].iface1 = iface2;
    edges[iedge2].iface1 = iface2;
    ipnext  = edges[ipnext].inext;
    istart = 0;
    nnode--; 

  //if (draw_flag)  draw_contour(4, 2, iedge1);

  }
}

void BooleanProcessor::modifyReference(int iface, int i1, int i2, int iref)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::modifyReference           Date:    13.03.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Modify reference to the neighbouring face                 *
 *                                                                     *
 ***********************************************************************/
{
  int iedge = faces[iface].iold;
  while (iedge > 0) {
    if (edges[iedge].i1 == i2 && edges[iedge].i2 == i1) {
      edges[iedge].iface2 = iref;
      return;
    }
    iedge = edges[iedge].inext;
  }
  processor_error = 1;
  std::cerr
    << "BooleanProcessor::modifyReference : could not find the edge, "
    << "iface=" << iface << ", i1,i2=" << i1 << "," << i2 << ", iref=" << iref
    << std::endl;
}

void BooleanProcessor::triangulateFace(int iface)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::triangulateFace           Date:    02.03.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Triangulation of an extended face                         *
 *                                                                     *
 ***********************************************************************/
{

  //   F I N D   M A X   C O M P O N E N T   O F   T H E   N O R M A L
  //   S E T  IX, IY, IZ

  Normal3D<double> normal = faces[iface].plane.normal();
  int ix, iy, iz = 0;
  if (std::abs(normal[1]) > std::abs(normal[iz])) iz = 1;
  if (std::abs(normal[2]) > std::abs(normal[iz])) iz = 2;
  if (normal[iz] > 0) {
    ix = (iz+1)%3; iy = (ix+1)%3;
  }else{
    iy = (iz+1)%3; ix = (iy+1)%3;
  }

  //   F I L L   L I S T S   O F   C O N T O U R S

  external_contours.clear();
  internal_contours.clear();
  double z;
  int    i1, i2, ifirst, iedge, icontour = faces[iface].iold;
  while (icontour > 0) { 
    iedge  = icontour;
    ifirst = edges[iedge].i1;
    z      = 0.0;
    for(;;) {
      if (iedge > 0) { 
        i1 = edges[iedge].i1;
        i2 = edges[iedge].i2;
        z += nodes[i1].v[ix]*nodes[i2].v[iy]-nodes[i2].v[ix]*nodes[i1].v[iy];
        if (ifirst != i2) {
          iedge = edges[iedge].inext;
          continue;
        }else{
          if (z > del*del) {
            external_contours.push_back(icontour);
          }else if (z < -del*del) {
            internal_contours.push_back(icontour);
          }else{  
            processor_error = 1;
            std::cerr
              << "BooleanProcessor::triangulateFace : too small contour"
              << std::endl;
          }
          icontour = edges[iedge].inext;
          edges[iedge].inext = 0;
          break;
        }
      }else{
        processor_error = 1;
        std::cerr
          << "BooleanProcessor::triangulateFace : broken contour"
          << std::endl;
        icontour = 0;
        break;
      }
    }
  }

  //   G E T   R I D   O F   I N T E R N A L   C O N T O U R S

  int kint, kext;
  for (kint=0; kint < internal_contours.size(); kint++) {
    for (kext=0; kext < external_contours.size(); kext++) {
      mergeContours(ix, iy, kext, kint); 
      if (internal_contours[kint] == 0) break;
    }
    if (kext == external_contours.size()) {
      processor_error = 1;
      std::cerr
        << "BooleanProcessor::triangulateFace : "
        << "could not merge internal contour " << kint 
        << std::endl;
    }      
  }

  //   T R I A N G U L A T E   C O N T O U R S

  int nface = faces.size();
  for (kext=0; kext < external_contours.size(); kext++) {
    triangulateContour(ix, iy, external_contours[kext]);
  }
  faces[iface].inew = UNSUITABLE_FACE;

  //   M O D I F Y   R E F E R E N C E S

  for (int ifa=nface; ifa<faces.size(); ifa++) {
    iedge = faces[ifa].iold;
    while (iedge > 0) {
      if (edges[iedge].iface1 != ifa) {
        processor_error = 1;
        std::cerr
          << "BooleanProcessor::triangulateFace : wrong reference to itself, "
          << "iface=" << ifa << ", iface1=" << edges[iedge].iface1
          << std::endl;
      }else if (edges[iedge].iface2 > 0) {
        modifyReference(edges[iedge].iface2,
                        edges[iedge].i1, edges[iedge].i2, ifa);
      }else if (edges[iedge].iface2 < 0) {
        edges[iedge].iface2 = edges[-edges[iedge].iface2].iface1;
      }
      iedge = edges[iedge].inext;
    }
  }
}

HepPolyhedron BooleanProcessor::createPolyhedron()
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::createPolyhedron()        Date:    14.03.00 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Create HepPolyhedron.                                     *
 *                                                                     *
 ***********************************************************************/
{
  int i, iedge, nnode = 0, nface = 0;

  //   R E N U M E R A T E   N O D E S   A N D   F A C E S

  for (i=1; i<nodes.size(); i++) nodes[i].s = 0;

  for (i=1; i<faces.size(); i++) {
    if (faces[i].inew == ORIGINAL_FACE) {
      faces[i].inew = ++nface;
      iedge = faces[i].iold;
      while (iedge > 0) {
        nodes[edges[iedge].i1].s = 1;
        iedge = edges[iedge].inext;
      }
    }else{
      faces[i].inew = 0;
    }
  } 

  for (i=1; i<nodes.size(); i++) {
    if (nodes[i].s == 1) nodes[i].s = ++nnode;
  }

  //   A L L O C A T E   M E M O R Y

  ExtPolyhedron polyhedron;
  if (nface == 0) return polyhedron;
  polyhedron.AllocateMemory(nnode, nface);

  //   S E T   N O D E S

  for (i=1; i<nodes.size(); i++) {
    if (nodes[i].s != 0)  polyhedron.pV[nodes[i].s] = nodes[i].v;
  }

  //   S E T   F A C E S

  int k, v[4], f[4];
  for (i=1; i<faces.size(); i++) {
    if (faces[i].inew == 0) continue;
    v[3] = f[3] = k = 0;
    iedge = faces[i].iold;
    while (iedge > 0) {
      if (k > 3) {
        std::cerr
          << "BooleanProcessor::createPolyhedron : too many edges"
          << std::endl;
        break;
      }
      v[k]  = nodes[edges[iedge].i1].s;
      if (edges[iedge].ivis < 0) v[k] = -v[k];
      f[k]  = faces[edges[iedge].iface2].inew;
      iedge = edges[iedge].inext;
      k++;
    }
    if (k < 3) {
      std::cerr
        << "BooleanProcessor::createPolyhedron : "
        << "face has only " << k << " edges"
        << std::endl;
    }
    polyhedron.pF[faces[i].inew] =
      G4Facet(v[0],f[0], v[1],f[1], v[2],f[2], v[3],f[3]);
  }
  return polyhedron;
}

HepPolyhedron BooleanProcessor::execute(int op,
                                        const HepPolyhedron & a,
                                        const HepPolyhedron & b)
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::execute                   Date:    10.12.99 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Execute boolean operation.                                *
 *                                                                     *
 ***********************************************************************/
{
  static int ishift = 0;
  static double shift[8][3] = {
    {  31,  23,  17},
    { -31, -23, -17},
    { -23,  17,  31},
    {  23, -17, -31},
    { -17, -31,  23},
    {  17,  31, -23},
    {  31, -23,  17},
    { -31,  23, -17}
  };           

  //   I N I T I A T E   P R O C E S S O R

  processor_error = 0;
  operation = op;
  nodes.clear(); nodes.push_back(CRAZY_POINT);
  edges.clear(); edges.push_back(ExtEdge());
  faces.clear(); faces.push_back(ExtFace());

  //   T A K E   P O L Y H E D R A

  ifaces1 = faces.size(); takePolyhedron(a,0,0,0);
  ifaces2 = faces.size(); takePolyhedron(b,0,0,0);

  if (processor_error) {             // corrapted polyhedron
    std::cerr
      << "BooleanProcessor: corrapted input polyhedron"
      << std::endl;
    return HepPolyhedron();
  }
  if (ifaces1 == ifaces2) {          // a is empty
    switch (operation) {
    case OP_UNION:
      return b;
    case OP_INTERSECTION:
      std::cerr 
        << "BooleanProcessor: intersection with empty polyhedron"
        << std::endl;
      return HepPolyhedron();
    case OP_SUBTRACTION:
      std::cerr
        << "BooleanProcessor: subtraction from empty polyhedron"
        << std::endl;
      return HepPolyhedron();
    }
  }
  if (ifaces2 == faces.size()) {     // b is empty
    switch (operation) {
    case OP_UNION:
      return a;
    case OP_INTERSECTION:
      std::cerr
        << "BooleanProcessor: intersection with empty polyhedron"
        << std::endl;
      return HepPolyhedron();
    case OP_SUBTRACTION:
      return a;
    }
  }

  //   S E T   I N I T I A L   M I N - M A X   A N D   T O L E R A N C E

  del = findMinMax();

  //   W O R K A R O U N D   T O   A V O I D   I E   A N D   E E
          
  double ddxx = del*shift[ishift][0];
  double ddyy = del*shift[ishift][1];
  double ddzz = del*shift[ishift][2];
  ishift++; if (ishift == 8) ishift = 0;

  operation = op;
  nodes.clear(); nodes.push_back(CRAZY_POINT);
  edges.clear(); edges.push_back(ExtEdge());
  faces.clear(); faces.push_back(ExtFace());

  ifaces1 = faces.size(); takePolyhedron(a,0,0,0);
  ifaces2 = faces.size(); takePolyhedron(b,ddxx,ddyy,ddzz);

  del = findMinMax();

  //   P R E S E L E C T   O U T S I D E   F A C E S

  iout1 = iout2 = 0;
  selectOutsideFaces(ifaces1, iout1);
  selectOutsideFaces(ifaces2, iout2);

  //   P R E S E L E C T   N O   I N T E R S E C T I O N   F A C E S

  int ifa1, ifa2;
  iunk1 = iunk2 = 0;
  if (iout1 != 0 || iout2 != 0) {
    for(;;) {
      ifa1 = iunk1;
      ifa2 = iunk2;
      selectOutsideFaces(ifaces1, iunk1);
      selectOutsideFaces(ifaces2, iunk2);
      if (iunk1 == ifa1 && iunk2 == ifa2) break;
      findMinMax();
    }
  }

#define PROCESSOR_ERROR \
std::cerr << "BooleanProcessor: boolean operation failed" << std::endl;\
return a;

  //   F I N D   N E W   E D G E S

  if (ifaces1 != 0 && ifaces2 != 0 ) {
    ifa1 = ifaces1;
    while (ifa1 > 0) {
      ifa2 = ifaces2;
      while (ifa2 > 0) {
        testFaceVsFace(ifa1, ifa2);
        ifa2 = faces[ifa2].inext;
      }
      ifa1 = faces[ifa1].inext;
    }
  }
  if (processor_error) { PROCESSOR_ERROR }

  //   C O N S T R U C T   N E W   F A C E S

  assembleNewFaces((operation == OP_INTERSECTION) ? 1 : 0, ifaces1);
  if (processor_error) { PROCESSOR_ERROR }
  assembleNewFaces((operation == OP_UNION) ? 0 : 1, ifaces2);
  if (processor_error) { PROCESSOR_ERROR }

  //   A S S E M B L E   S U I T A B L E   F A C E S

  initiateLists();
  for (;;) {
    assemblePolyhedra();
    if (unknown_faces.front() != 0) {
      processor_error = 1;
      std::cerr
        << "BooleanProcessor::execute : unknown faces !!!"
        << std::endl;
    }
    break;
  }
  if (processor_error) { PROCESSOR_ERROR }

  //   T R I A N G U L A T E   A C C E P T E D   F A C E S 

  ifa1 = result_faces.front();
  while (ifa1 > 0) {
    ifa2 = ifa1;
    ifa1 = faces[ifa2].inext;
    if (faces[ifa2].inew == NEW_FACE) triangulateFace(ifa2);
    if (processor_error) { PROCESSOR_ERROR }
  }

  //   C R E A T E   P O L Y H E D R O N
  
  return createPolyhedron();
}


//#include <cfortran.h>
//#include <higz.h>
//#include "zbuf.h"
//void BooleanProcessor::draw()
/***********************************************************************
 *                                                                     *
 * Name: BooleanProcessor::draw                      Date:    10.12.99 *
 * Author: E.Chernyaev                               Revised:          *
 *                                                                     *
 * Function: Draw                                                      *
 *                                                                     *
 ***********************************************************************/
/*
{
  int II;
  int   icol, i1, i2, iedge, iface, ilist[4];
  float p1[3], p2[3];

  ilist[0] = ifaces1;
  ilist[1] = ifaces2;
  ilist[2] = iout1;
  ilist[3] = iout2;

  for (int i=0; i<4; i++) { 

    if (i == 0) cout << "========= Ifaces_1" << endl;
    if (i == 1) cout << "========= Ifaces_2" << endl;
    if (i == 2) cout << "========= Iout_1" << endl;
    if (i == 3) cout << "========= Iout_2" << endl;

    icol = i+1;
    iface = ilist[i];
    while (iface > 0) {

      cout << "iface = " << iface << endl;
      cout << "--- iold" << endl;

      iedge = faces[iface].iold;
      icol = 2;

      while (iedge > 0) {

        cout << "  iegde = " << iedge
             << " i1,i2 =" << edges[iedge].i1 << "," << edges[iedge].i2
             << " iface1,iface2 = "
             << edges[iedge].iface1 << "," << edges[iedge].iface2
             << endl;

        i1 = edges[iedge].i1;
        p1[0] = nodes[i1].v.x();
        p1[1] = nodes[i1].v.y();
        p1[2] = nodes[i1].v.z();
        IHWTON(p1,p1);
        i2 = edges[iedge].i2;
        p2[0] = nodes[i2].v.x();
        p2[1] = nodes[i2].v.y();
        p2[2] = nodes[i2].v.z();
        IHWTON(p2,p2);
//        icol =  (edges[iedge].ivis > 0) ? 1 : 2;
        IHZLIN(icol,p1[0],p1[1],p1[2], p2[0],p2[1],p2[2]);
        iedge = edges[iedge].inext;
      }
      
      cout << "--- inew" << endl;

      iedge = faces[iface].inew;
      icol = 3;

      while (iedge > 0) {

        cout << "  iegde = " << iedge
             << " i1,i2 =" << edges[iedge].i1 << "," << edges[iedge].i2
             << " iface1,iface2 = "
             << edges[iedge].iface1 << "," << edges[iedge].iface2
             << endl;

        i1 = edges[iedge].i1;
        p1[0] = nodes[i1].v.x();
        p1[1] = nodes[i1].v.y();
        p1[2] = nodes[i1].v.z();
        IHWTON(p1,p1);
        i2 = edges[iedge].i2;
        p2[0] = nodes[i2].v.x();
        p2[1] = nodes[i2].v.y();
        p2[2] = nodes[i2].v.z();
        IHWTON(p2,p2);
//        icol =  (edges[iedge].ivis > 0) ? 1 : 2;
        IHZLIN(icol,p1[0],p1[1],p1[2], p2[0],p2[1],p2[2]);
        iedge = edges[iedge].inext;
      }
      iface = faces[iface].inext;

      IHZTOX(0,100,100);
      ixupdwi(0);
      cin >> II;
      ixclrwi();
      IHZCLE(0);
    }
  }
}
*/

/*
//--------------------------------------------------------------------
void
BooleanProcessor::draw_edge(int icol, int iedge) {
  int   i1, i2;
  float p1[3], p2[3];

  i1 = edges[iedge].i1;
  p1[0] = nodes[i1].v.x();
  p1[1] = nodes[i1].v.y();
  p1[2] = nodes[i1].v.z();
  IHWTON(p1,p1);
  i2 = edges[iedge].i2;
  p2[0] = nodes[i2].v.x();
  p2[1] = nodes[i2].v.y();
  p2[2] = nodes[i2].v.z();
  IHWTON(p2,p2);
  IHZLIN(icol,p1[0],p1[1],p1[2], p2[0],p2[1],p2[2]);
}

//--------------------------------------------------------------------
void
BooleanProcessor::draw_contour(int i1col, int i2col, int ihead) {
  int iedge, icol;
  iedge = ihead;
  while (iedge > 0) {
    icol = (edges[iedge].ivis > 0) ? i1col : i2col;
    draw_edge(icol, iedge);
    iedge = edges[iedge].inext;
  }

  IHZTOX(0,100,100);
  ixupdwi(0);

  int i; 
  std::cin >> i; 
}

//--------------------------------------------------------------------
void
BooleanProcessor::print_face(int iface) {
  cout.precision(3);
  cout << "\n====== Face N " << iface << endl;
  cout << "iedges[4] = "
       << faces[iface].iedges[0] << ", "
       << faces[iface].iedges[1] << ", "
       << faces[iface].iedges[2] << ", "
       << faces[iface].iedges[3] << endl;
  cout << "rmin[3] = "
       << faces[iface].rmin[0] << ", "
       << faces[iface].rmin[1] << ", "
       << faces[iface].rmin[2] << endl;
  cout << "rmax[3] = "
       << faces[iface].rmax[0] << ", "
       << faces[iface].rmax[1] << ", "
       << faces[iface].rmax[2] << endl;
  cout << "iprev,inext = "
       << faces[iface].iprev << ", "
       << faces[iface].inext << endl;
  cout << "iold = " << faces[iface].iold << endl;
  for(int i = faces[iface].iold; i != 0;) {
    print_edge(i);
    i = edges[abs(i)].inext;
  }

  cout << "inew = ";
  switch (faces[iface].inew) {
  case UNKNOWN_FACE:
    cout << "UNKNOWN_FACE" << endl;
    break;
  case ORIGINAL_FACE:
    cout << "ORIGINAL_FACE" << endl;
    break;
  case NEW_FACE:
    cout << "NEW_FACE" << endl;
    break;
  case UNSUITABLE_FACE:
    cout << "UNSUITABLE_FACE" << endl;
    break;
  case DEFECTIVE_FACE:
    cout << "DEFECTIVE_FACE" << endl;
    break;
  default:
    cout << faces[iface].inew << endl;
    for(int k = faces[iface].inew; k != 0;) {
      print_edge(k);
      k = edges[abs(k)].inext;
    }
  }
}

//--------------------------------------------------------------------
void
BooleanProcessor::print_edge(int iedge) {
  cout << "==== Edge N " << iedge << endl;
  int i = std::abs(iedge);
  int i1 = edges[i].i1;
  int i2 = edges[i].i2;
  cout << "node[" << i1 << "] = " 
       << nodes[i1].v.x() << ", "
       << nodes[i1].v.y() << ", "
       << nodes[i1].v.z() << endl;

  cout << "node[" << i2 << "] = " 
       << nodes[i2].v.x() << ", "
       << nodes[i2].v.y() << ", "
       << nodes[i2].v.z() << endl;

  cout << "iface1,iface2,ivis,inext = "
       << edges[i].iface1 << ", "
       << edges[i].iface2 << ", "
       << edges[i].ivis   << ", "
       << edges[i].inext  << endl;
}
*/