g4tools  5.4.0
pbp.icc
Go to the documentation of this file.
1 // Copyright (C) 2010, Guy Barrand. All rights reserved.
2 // See the file tools.license for terms.
3 
4 /***********************************************************************
5  * *
6  * Name: BooleanProcessor Date: 10.12.99 *
7  * Author: E.Chernyaev Revised: *
8  * *
9  * Function: Internal class for executing boolean operations *
10  * on Polyhedra *
11  * *
12  ***********************************************************************/
13 
14 #include "../lina/plane"
15 
16 //#define TOOLS_HEP_BP_OUT_ERR
17 //#define TOOLS_HEP_BP_CHECK_INDEX
18 
19 //#define TOOLS_HEP_BP_VERBOSE
20 //#define TOOLS_HEP_BP_NOT_OPT
21 
22 #if defined(TOOLS_HEP_BP_OUT_ERR) || defined(TOOLS_HEP_BP_VERBOSE) || defined(TOOLS_HEP_BP_CHECK_INDEX)
23 #include <iostream>
24 #endif
25 
26 namespace tools {
27 namespace hep {
28 
29 typedef plane<vec3d> HVPlane3D;
30 
31 const double GRANULARITY = 10.e+5;
32 
33 // Operations
34 const int OP_UNION = 0;
35 const int OP_INTERSECTION = 1;
36 const int OP_SUBTRACTION = 2;
37 
38 // Face vs face statuses
39 const int OUT_OF_PLANE = 0;
40 const int ON_PLANE = 1;
41 const int INTERSECTION = 2;
42 const int EDGE = 3;
43 const int NON_PLANAR_FACE = 4;
44 
45 // Face statuses
46 const int UNKNOWN_FACE = 0;
47 const int ORIGINAL_FACE = -1;
48 const int NEW_FACE = -2;
49 const int UNSUITABLE_FACE = -3;
50 const int DEFECTIVE_FACE = -4;
51 
52 // -------------------------------------------- Simplified STL vector ---
53 //G.Barrand : begin
54 //#include <vector>
55 //G.Barrand : end
56 
57 // ---------------------------------------------------- Extended node ---
58 class ExtNode {
59 #ifdef TOOLS_MEM
60  TOOLS_SCLASS(tools::hep::ExtNode)
61 #endif
62  public:
63  HVPoint3D v;
64  int s;
65 
66  public:
67  ExtNode(HVPoint3D vertex=HVPoint3D(), int status=0)
68  : v(vertex), s(status) {
69 #ifdef TOOLS_MEM
70  mem::increment(s_class().c_str());
71 #endif
72  }
73  ~ExtNode() {
74 #ifdef TOOLS_MEM
75  mem::decrement(s_class().c_str());
76 #endif
77  }
78 
79  ExtNode(const ExtNode & node) : v(node.v), s(node.s) {
80 #ifdef TOOLS_MEM
81  mem::increment(s_class().c_str());
82 #endif
83  }
84 
85  ExtNode & operator=(const ExtNode & node) {
86  v = node.v;
87  s = node.s;
88  return *this;
89  }
90 };
91 
92 // ---------------------------------------------------- Extended edge ---
93 class ExtEdge {
94 #ifdef TOOLS_MEM
95  TOOLS_SCLASS(tools::hep::ExtEdge)
96 #endif
97  public:
98  int i1, i2; // end points
99  int iface1; // native face
100  int iface2; // neighbouring face
101  int ivis; // visibility: +1 (visible), -1 (invisible)
102  int inext; // index of next edge
103 
104  public:
105  ExtEdge(int k1=0, int k2=0, int kface1=0, int kface2=0, int kvis=0) :
106  i1(k1), i2(k2), iface1(kface1), iface2(kface2), ivis(kvis), inext(0) {
107 #ifdef TOOLS_MEM
108  mem::increment(s_class().c_str());
109 #endif
110  }
111 
112  ~ExtEdge() {
113 #ifdef TOOLS_MEM
114  mem::decrement(s_class().c_str());
115 #endif
116  }
117 
118  ExtEdge(const ExtEdge & edge) :
119  i1(edge.i1), i2(edge.i2), iface1(edge.iface1), iface2(edge.iface2),
120  ivis(edge.ivis), inext(edge.inext) {
121 #ifdef TOOLS_MEM
122  mem::increment(s_class().c_str());
123 #endif
124  }
125 
126  ExtEdge & operator=(const ExtEdge & edge) {
127  i1 = edge.i1;
128  i2 = edge.i2;
129  iface1 = edge.iface1;
130  iface2 = edge.iface2;
131  ivis = edge.ivis;
132  inext = edge.inext;
133  return *this;
134  }
135 
136  void invert() {
137  int w;
138  //#define SWAP(A,B) w = A; A = B; B = w
139  //SWAP(i1, i2);
140  w = i1;i1=i2;i2=w;
141  }
142 };
143 
144 // ---------------------------------------------------- Extended face ---
145 class ExtFace {
146  private:
147 #ifdef TOOLS_MEM
148  TOOLS_SCLASS(tools::hep::ExtFace)
149 #endif
150  std::vector<ExtEdge>& edges; //G.Barrand
151  public:
152  int iedges[4]; // indices of original edges
153  HVPlane3D plane; // face plane
154  double rmin[3], rmax[3]; // bounding box
155  int iold; // head of the list of the original edges
156  int inew; // head of the list of the new edges
157  int iprev; // index of previous face
158  int inext; // index of next face
159 
160  public:
161  //G.Barrand : ExtFace(int iedge=0) : iold(iedge), inew(0), iprev(iprev), inext(0) {}
162  ExtFace(std::vector<ExtEdge>& a_edges,int iedge)
163  : edges(a_edges), iold(iedge), inew(0), iprev(0), inext(0) {
164 #ifdef TOOLS_MEM
165  mem::increment(s_class().c_str());
166 #endif
167  //G.Barrand : initialize arrays to quiet valgrind.
168  {for (int i=0; i<4; i++) { iedges[i] = 0; }}
169  {for (int i=0; i<3; i++) { rmin[i] = 0; rmax[i] = 0; }}
170  }
171  ~ExtFace() {
172 #ifdef TOOLS_MEM
173  mem::decrement(s_class().c_str());
174 #endif
175  }
176 
177  ExtFace(const ExtFace & face) :
178  edges(face.edges), //G.Barrand
179  plane(face.plane), iold(face.iold), inew(face.inew),
180  iprev(face.iprev), inext(face.inext)
181  {
182 #ifdef TOOLS_MEM
183  mem::increment(s_class().c_str());
184 #endif
185  int i;
186  for (i=0; i<4; i++) { iedges[i] = face.iedges[i]; }
187  for (i=0; i<3; i++) { rmin[i] = face.rmin[i]; rmax[i] = face.rmax[i]; }
188  }
189 
190  ExtFace & operator=(const ExtFace & face) {
191  //FIXME : edges(face.edges) ???? //G.Barrand
192  int i;
193  for (i=0; i<4; i++) { iedges[i] = face.iedges[i]; }
194  plane = face.plane;
195  for (i=0; i<3; i++) { rmin[i] = face.rmin[i]; rmax[i] = face.rmax[i]; }
196  iold = face.iold;
197  inew = face.inew;
198  iprev = face.iprev;
199  inext = face.inext;
200  return *this;
201  }
202 
203  bool invert(); //G.Barrand : return bool.
204 
205  //G.Barrand
206 #ifdef TOOLS_HEP_BP_CHECK_INDEX
207  bool check_edge_index(const std::string& a_where,int a_index) {
208  if(a_index<0) {
209  std::cerr << "ExtFace::check_edge_index :"
210  << " " << a_where << " :"
211  << " i<0."
212  << std::endl;
213  return false;
214  }
215  if(a_index>=int(edges.size())) {
216  std::cerr << "ExtFace::check_edge_index :"
217  << " " << a_where << " :"
218  << " i>=sz."
219  << " i=" << a_index
220  << " sz=" << edges.size()
221  << " cap=" << edges.capacity()
222  << std::endl;
223  return false;
224  }
225  return true;
226  }
227 #endif
228 };
229 
230 // ---------------------------------------------------- Global arrays ---
231 //G.Barrand : MacIntel : crash with g++-4.0.1 with -O on some subtract.
232 // Anyway static of objects is proved to be not safe.
233 // We put the below vector as members of BooleanProcessor.
234 //GB static std::vector<ExtNode> nodes; // vector of nodes
235 //GB static std::vector<ExtEdge> edges; // vector of edges
236 //GB static std::vector<ExtFace> faces; // vector of faces
237 
238 // ---------------------------------------------------- List of faces ---
239 class FaceList {
240 #ifdef TOOLS_MEM
241  TOOLS_SCLASS(tools::hep::FaceList)
242 #endif
243  private:
244  std::vector<ExtFace>& faces; //G.Barrad : end
245  private:
246  int ihead;
247  int ilast;
248 
249  public:
250  //G.Barrand : FaceList() : ihead(0), ilast(0) {}
251  FaceList(std::vector<ExtFace>& a_faces) : faces(a_faces),ihead(0),ilast(0) {
252 #ifdef TOOLS_MEM
253  mem::increment(s_class().c_str());
254 #endif
255  }
256  ~FaceList() {
257 #ifdef TOOLS_MEM
258  mem::decrement(s_class().c_str());
259 #endif
260  }
261  protected:
262  FaceList(const FaceList& a_from)
263  :faces(a_from.faces),ihead(0),ilast(0){
264 #ifdef TOOLS_MEM
265  mem::increment(s_class().c_str());
266 #endif
267  }
268 
269  FaceList& operator=(const FaceList&) {return *this;}
270  public:
271  void clean() { ihead = 0; ilast = 0; }
272  int front() { return ihead; }
273 
274  void push_back(int i) {
275  if (ilast == 0) { ihead = i; } else { faces[ilast].inext = i; }
276  ExtFace& face = faces[i]; //G.Barrand : optimize.
277  face.iprev = ilast;
278  face.inext = 0;
279  ilast = i;
280  }
281 
282  void remove(int i) {
283  ExtFace& face = faces[i]; //G.Barrand : optimize.
284  if (ihead == i) {
285  ihead = face.inext;
286  }else{
287  faces[face.iprev].inext = face.inext;
288  }
289  if (ilast == i) {
290  ilast = face.iprev;
291  }else{
292  faces[face.inext].iprev = face.iprev;
293  }
294  face.iprev = 0;
295  face.inext = 0;
296  }
297 };
298 
299 // --------------------- Polyhedron with extended access to
300 // its members from the BooleanProcessor class ---
301 class ExtPolyhedron : public polyhedron {
302  friend class BooleanProcessor;
303  virtual polyhedron& operator = (const polyhedron& from) {
304  return polyhedron::operator = (from);
305  }
306 };
307 
308 // ----------------------------------------- Boolean processor class ---
309 class BooleanProcessor {
310 #ifdef TOOLS_MEM
311  TOOLS_SCLASS(tools::hep::BooleanProcessor)
312 #endif
313  private:
314  //static int ishift; //G.Barrand
315  std::vector<ExtNode> nodes; // vector of nodes //G.Barrand
316  std::vector<ExtEdge> edges; // vector of edges //G.Barrand
317  std::vector<ExtFace> faces; // vector of faces //G.Barrand
318  private:
319  //int processor_error; // is set in case of error
320  int operation; // 0 (union), 1 (intersection), 2 (subtraction)
321  int ifaces1, ifaces2; // lists of faces
322  int iout1, iout2; // lists of faces with status "out"
323  int iunk1, iunk2; // lists of faces with status "unknown"
324  double rmin[3], rmax[3]; // intersection of bounding boxes
325  double del; // precision (tolerance)
326 
327  FaceList result_faces; // list of accepted faces
328  FaceList suitable_faces; // list of suitable faces
329  FaceList unsuitable_faces; // list of unsuitable faces
330  FaceList unknown_faces; // list of unknown faces
331 
332  std::vector<int> external_contours; // heads of external contours
333  std::vector<int> internal_contours; // heads of internal contours
334 
335  public:
336  private:
337  bool takePolyhedron(const polyhedron & p, double, double, double);
338  double findMinMax();
339  void selectOutsideFaces(int & ifaces, int & iout);
340  int testFaceVsPlane(ExtEdge & edge);
341  void renumberNodes(int & i1, int & i2, int & i3, int & i4);
342  int testEdgeVsEdge(ExtEdge & edge1, ExtEdge & edge2);
343  void removeJunkNodes() { while(nodes.back().s != 0) nodes.pop_back(); }
344  void divideEdge(int & i1, int & i2);
345  void insertEdge(const ExtEdge & edge);
346  void caseII(ExtEdge & edge1, ExtEdge & edge2);
347  bool caseIE(ExtEdge & edge1, ExtEdge & edge2);
348  bool caseEE(ExtEdge & edge1, ExtEdge & edge2);
349  bool testFaceVsFace(int iface1, int iface2);
350  void invertNewEdges(int iface);
351  void checkDoubleEdges(int iface);
352  bool assembleFace(int what, int iface);
353  bool assembleNewFaces(int what, int ihead);
354  bool initiateLists();
355  bool assemblePolyhedra();
356  void findABC(double x1, double y1, double x2, double y2,
357  double &a, double &b, double &c) const;
358  int checkDirection(double *x, double *y) const;
359  int checkIntersection(int ix, int iy, int i1, int i2) const;
360  void mergeContours(int ix, int iy, int kext, int kint);
361  bool checkTriangle(int iedge1, int iedge2, int ix, int iy) const;
362  bool triangulateContour(int ix, int iy, int ihead);
363  bool modifyReference(int iface, int i1, int i2, int iref);
364  bool triangulateFace(int iface);
365  polyhedron createPolyhedron();
366 
367  public:
368  //G.Barrand : BooleanProcessor() {}
369  BooleanProcessor() //G.Barrand
370  //:processor_error(0) //G.Barrand
371  :result_faces(faces)
372  ,suitable_faces(faces)
373  ,unsuitable_faces(faces)
374  ,unknown_faces(faces)
375  {
376 #ifdef TOOLS_MEM
377  mem::increment(s_class().c_str());
378 #endif
379  }
380 
381  ~BooleanProcessor() {
382 #ifdef TOOLS_MEM
383  mem::decrement(s_class().c_str());
384 #endif
385  }
386 private:
387  BooleanProcessor(const BooleanProcessor&)
388  //:processor_error(0) //G.Barrand
389  :result_faces(faces)
390  ,suitable_faces(faces)
391  ,unsuitable_faces(faces)
392  ,unknown_faces(faces)
393  {
394 #ifdef TOOLS_MEM
395  mem::increment(s_class().c_str());
396 #endif
397  }
398  BooleanProcessor& operator=(const BooleanProcessor&){return *this;}
399 public:
400  polyhedron execute(int op,
401  const polyhedron &a,
402  const polyhedron &b,
403  int& err);
404 
405  //G.Barrand : begin : drawing disconnected.
406  //void draw();
407  //void draw_edge(int, int);
408  //void draw_contour(int, int, int);
409  //void draw_faces(int, int, int);
410  //void print_face(int);
411  //void print_edge(int);
412  //G.Barrand : end
413 
414  //int get_processor_error() const {return processor_error;}
415 
416  void dump(std::ostream&); //G.Barrand
417  static int& get_shift(); //G.Barrand
418  static void set_shift(int); //G.Barrand
419  static void inc_shift(); //G.Barrand
420  static int get_num_shift(); //G.Barrand
421 
422 private:
423  bool inc_try_count(int,const polyhedron&,const polyhedron&,int&); //GB
424 
425  static HVPoint3D CRAZY_POINT(){return HVPoint3D(-10.e+6, -10.e+6, -10.e+6);}
426 #ifdef TOOLS_HEP_BP_OUT_ERR
427  static void PROCESSOR_ERROR(int a_what){
428  std::cerr << "BooleanProcessor: boolean operation problem (" << a_what
429  << "). Try again with other shifts."
430  << std::endl;
431  }
432 #else
433  static void PROCESSOR_ERROR(int){}
434 #endif
435 
436  //G.Barrand :
437 #ifdef TOOLS_HEP_BP_CHECK_INDEX
438  bool check_face_index(const std::string& a_where,int a_index) {
439  if(a_index<0) {
440  std::cerr << "BooleanProcessor::check_face_index :"
441  << " " << a_where << " :"
442  << " i<0."
443  << std::endl;
444  return false;
445  }
446  if(a_index>=int(faces.size())) {
447  std::cerr << "BooleanProcessor::check_face_index :"
448  << " " << a_where << " :"
449  << " i>=sz."
450  << " i=" << a_index
451  << " sz=" << faces.size()
452  << " cap=" << faces.capacity()
453  << std::endl;
454  return false;
455  }
456  return true;
457  }
458  bool check_edge_index(const std::string& a_where,int a_index) {
459  if(a_index<0) {
460  std::cerr << "BooleanProcessor::check_edge_index :"
461  << " " << a_where << " :"
462  << " i<0."
463  << std::endl;
464  return false;
465  }
466  if(a_index>=int(edges.size())) {
467  std::cerr << "BooleanProcessor::check_edge_index :"
468  << " " << a_where << " :"
469  << " i>=sz."
470  << " i=" << a_index
471  << " sz=" << edges.size()
472  << " cap=" << edges.capacity()
473  << std::endl;
474  return false;
475  }
476  return true;
477  }
478 #endif
479 
480 };
481 
482 inline bool ExtFace::invert() //G.Barrand : return bool.
483 /***********************************************************************
484  * *
485  * Name: ExtFace::invert() Date: 28.02.00 *
486  * Author: E.Chernyaev Revised: *
487  * *
488  * Function: Invert face *
489  * *
490  ***********************************************************************/
491 {
492  int iEprev, iEcur, iEnext;
493 
494  iEprev = 0; iEcur = iold;
495  while (iEcur > 0) {
496 #ifdef TOOLS_HEP_BP_CHECK_INDEX
497  if(!check_edge_index("ExtFace::invert(1)",iEcur)) return false;
498 #endif
499  ExtEdge& edge = edges[iEcur]; //G.Barrand : optimize.
500  edge.invert();
501  iEnext = edge.inext;
502  edge.inext = iEprev;
503  iEprev = iEcur;
504  iEcur = iEnext;
505  }
506  if (iold > 0) iold = iEprev;
507 
508  iEprev = 0; iEcur = inew;
509  while (iEcur > 0) {
510 #ifdef TOOLS_HEP_BP_CHECK_INDEX
511  if(!check_edge_index("ExtFace::invert(2)",iEcur)) return false;
512 #endif
513  ExtEdge& edge = edges[iEcur]; //G.Barrand : optimize.
514  edge.invert();
515  iEnext = edge.inext;
516  edge.inext = iEprev;
517  iEprev = iEcur;
518  iEcur = iEnext;
519  }
520  if (inew > 0) inew = iEprev;
521 
522  //plane = HVPlane3D(-plane.a(), -plane.b(), -plane.c(), -plane.d());
523 
524  HVNormal3D n = plane.normal();
525  n *= -1;
526  plane.set(n, -plane.distance_from_origin());
527 
528  return true;
529 }
530 
531 inline
532 bool BooleanProcessor::takePolyhedron(const polyhedron & p,
533  double dx, double dy, double dz)
534 /***********************************************************************
535  * *
536  * Name: BooleanProcessor::takePolyhedron Date: 16.12.99 *
537  * Author: E.Chernyaev Revised: *
538  * *
539  * Function: Transfer Polyhedron to internal representation *
540  * *
541  ***********************************************************************/
542 {
543  int i, k, nnode, iNodes[5], iVis[4], iFaces[4];
544  int dnode = int(nodes.size()) - 1;
545  int dface = int(faces.size()) - 1;
546 
547  // S E T N O D E S
548 
549  // for (i=1; i <= p.GetNoVertices(); i++) {
550  // nodes.push_back(ExtNode(p.GetVertex(i)));
551  // }
552 
553  HVPoint3D ppp;
554  for (i=1; i <= p.GetNoVertices(); i++) {
555 #ifndef TOOLS_HEP_BP_NOT_OPT //G.Barrand
556  nodes.push_back(ExtNode());
557  ExtNode& node = nodes.back();
558  node.v = p.GetVertexFast(i);
559  node.v.add(dx,dy,dz);
560 #else
561  ppp = p.GetVertexFast(i);
562  ppp += HVPoint3D(dx,dy,dz);
563  nodes.push_back(ExtNode(ppp));
564 #endif
565  }
566 
567  // S E T F A C E S
568 
569  HVNormal3D plane_n; //G.Barrand
570  HVPoint3D plane_p; //G.Barrand
571 
572  for (int iface=1; iface <= p.GetNoFacets(); iface++) {
573  faces.push_back(ExtFace(edges,int(edges.size())));
574 
575  // S E T F A C E N O D E S
576 
577  p.GetFacet(iface, nnode, iNodes, iVis, iFaces);
578  for (i=0; i<nnode; i++) {
579  //if (iNodes[i] < 1 || iNodes[i] > p.GetNoVertices()) processor_error = 1;
580  //if (iFaces[i] < 1 || iFaces[i] > p.GetNoFacets()) processor_error = 1;
581 
582  if (iNodes[i] < 1 || iNodes[i] > p.GetNoVertices()) { //G.Barrand
583  //G.Barrand : processor_error = 1;
584 #ifdef TOOLS_HEP_BP_OUT_ERR
585  std::cerr
586  << "BooleanProcessor::takePolyhedron : problem 1."
587  << " nnode " << nnode
588  << " p.GetNoVertices() " << p.GetNoVertices()
589  << " p.GetNoFacets() " << p.GetNoFacets()
590  << " iface " << iface
591  << " i " << i
592  << " iNodes[i] " << iNodes[i]
593  << std::endl;
594 #endif
595  return false; //G.Barrand
596  }
597  if (iFaces[i] < 1 || iFaces[i] > p.GetNoFacets()) { //G.Barrand
598  //G.Barrand : processor_error = 1;
599 #ifdef TOOLS_HEP_BP_OUT_ERR
600  std::cerr
601  << "BooleanProcessor::takePolyhedron : problem 2."
602  << " p.GetNoVertices() " << p.GetNoVertices()
603  << " p.GetNoFacets() " << p.GetNoFacets()
604  << " iface " << iface
605  << " i " << i
606  << " iFaces[i] " << iFaces[i]
607  << std::endl;
608 #endif
609  return false; //G.Barrand
610  }
611  iNodes[i] += dnode;
612  iFaces[i] += dface;
613  }
614 
615  // S E T E D G E S
616 
617  iNodes[nnode] = iNodes[0];
618  faces.back().iedges[3] = 0;
619  for (i=0; i<nnode; i++) {
620  faces.back().iedges[i] = int(edges.size());
621  edges.push_back(ExtEdge(iNodes[i], iNodes[i+1],
622  iface+dface, iFaces[i], iVis[i]));
623  edges.back().inext = int(edges.size());
624  }
625  edges.back().inext = 0;
626 
627  // S E T F A C E M I N - M A X
628 
629  ExtFace& face = faces.back(); //G.Barrand : optimize.
630  for (i=0; i<3; i++) {
631  face.rmin[i] = nodes[iNodes[0]].v[i];
632  face.rmax[i] = nodes[iNodes[0]].v[i];
633  }
634  for (i=1; i<nnode; i++) {
635  ExtNode& node = nodes[iNodes[i]]; //G.Barrand : optimize.
636  for (k=0; k<3; k++) {
637  if (face.rmin[k] > node.v[k])
638  face.rmin[k] = node.v[k];
639  if (face.rmax[k] < node.v[k])
640  face.rmax[k] = node.v[k];
641  }
642  }
643 
644  // S E T F A C E P L A N E
645 
646  (nodes[iNodes[2]].v-nodes[iNodes[0]].v).cross(nodes[iNodes[3]].v-nodes[iNodes[1]].v,plane_n);
647 
648  plane_p.set_value(0,0,0);
649  for (i=0; i<nnode; i++) { plane_p += nodes[iNodes[i]].v; }
650  plane_p *= (1./nnode);
651 
652  //G.Barrand : faces.back().plane = HVPlane3D(n.unit(), p);
653  faces.back().plane.set(plane_n,plane_p); //G.Barrand
654 
655  // 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
656 
657  faces.back().inext = int(faces.size());
658  }
659  faces.back().inext = 0;
660 
661  return true; //G.Barrand
662 }
663 
664 inline
665 double BooleanProcessor::findMinMax()
666 /***********************************************************************
667  * *
668  * Name: BooleanProcessor::findMinMax Date: 16.12.99 *
669  * Author: E.Chernyaev Revised: *
670  * *
671  * Function: Find min-max (bounding) boxes for polyhedra *
672  * *
673  ***********************************************************************/
674 {
675  if (ifaces1 == 0 || ifaces2 == 0) return 0;
676 
677  int i, iface;
678  double rmin1[3], rmax1[3];
679  double rmin2[3], rmax2[3];
680 
681  // F I N D B O U N D I N G B O X E S
682 
683  for (i=0; i<3; i++) {
684  rmin1[i] = faces[ifaces1].rmin[i];
685  rmax1[i] = faces[ifaces1].rmax[i];
686  rmin2[i] = faces[ifaces2].rmin[i];
687  rmax2[i] = faces[ifaces2].rmax[i];
688  }
689 
690  iface = faces[ifaces1].inext;
691  while(iface > 0) {
692  ExtFace& face = faces[iface]; //G.Barrand
693  for (i=0; i<3; i++) {
694  if (rmin1[i] > face.rmin[i]) rmin1[i] = face.rmin[i];
695  if (rmax1[i] < face.rmax[i]) rmax1[i] = face.rmax[i];
696  }
697  iface = face.inext;
698  }
699 
700  iface = faces[ifaces2].inext;
701  while(iface > 0) {
702  ExtFace& face = faces[iface]; //G.Barrand
703  for (i=0; i<3; i++) {
704  if (rmin2[i] > face.rmin[i]) rmin2[i] = face.rmin[i];
705  if (rmax2[i] < face.rmax[i]) rmax2[i] = face.rmax[i];
706  }
707  iface = face.inext;
708  }
709 
710  // 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
711 
712  for (i=0; i<3; i++) {
713  rmin[i] = (rmin1[i] > rmin2[i]) ? rmin1[i] : rmin2[i];
714  rmax[i] = (rmax1[i] < rmax2[i]) ? rmax1[i] : rmax2[i];
715  }
716 
717  // F I N D T O L E R A N C E
718 
719  double del1 = 0;
720  double del2 = 0;
721  for (i=0; i<3; i++) {
722  if ((rmax1[i]-rmin1[i]) > del1) del1 = rmax1[i]-rmin1[i];
723  if ((rmax2[i]-rmin2[i]) > del2) del2 = rmax2[i]-rmin2[i];
724  }
725  return ((del1 < del2) ? del1 : del2) / GRANULARITY;
726 }
727 
728 inline
729 void BooleanProcessor::selectOutsideFaces(int & ifaces, int & iout)
730 /***********************************************************************
731  * *
732  * Name: BooleanProcessor::selectOutsideFaces Date: 10.01.00 *
733  * Author: E.Chernyaev Revised: *
734  * *
735  * Function: Preselection of outside faces *
736  * *
737  ***********************************************************************/
738 {
739  int i, outflag, iface = ifaces, *prev;
740  HVPoint3D mmbox[8] = { HVPoint3D(rmin[0],rmin[1],rmin[2]),
741  HVPoint3D(rmax[0],rmin[1],rmin[2]),
742  HVPoint3D(rmin[0],rmax[1],rmin[2]),
743  HVPoint3D(rmax[0],rmax[1],rmin[2]),
744  HVPoint3D(rmin[0],rmin[1],rmax[2]),
745  HVPoint3D(rmax[0],rmin[1],rmax[2]),
746  HVPoint3D(rmin[0],rmax[1],rmax[2]),
747  HVPoint3D(rmax[0],rmax[1],rmax[2]) };
748  prev = &ifaces;
749  while (iface > 0) {
750 
751  // B O U N D I N G B O X vs B O U N D I N G B O X
752 
753  outflag = 0;
754  ExtFace& face = faces[iface]; //G.Barrand : optimize.
755  for (i=0; i<3; i++) {
756  if (face.rmin[i] > rmax[i] + del) { outflag = 1; break; }
757  if (face.rmax[i] < rmin[i] - del) { outflag = 1; break; }
758  }
759 
760  // B O U N D I N G B O X vs P L A N E
761 
762  if (outflag == 0) {
763  int npos = 0, nneg = 0;
764  double d;
765  for (i=0; i<8; i++) {
766  d = face.plane.distance(mmbox[i]); //G.Barrand : optimize
767  if (d > +del) npos++;
768  if (d < -del) nneg++;
769  }
770  if (npos == 8 || nneg == 8) outflag = 1;
771  }
772 
773  // U P D A T E L I S T S
774 
775  if (outflag == 1) {
776  *prev = face.inext;
777  face.inext = iout;
778  iout = iface;
779  }else{
780  prev = &face.inext;
781  }
782 
783  iface = *prev;
784  }
785 }
786 
787 inline
788 int BooleanProcessor::testFaceVsPlane(ExtEdge & edge)
789 /***********************************************************************
790  * *
791  * Name: BooleanProcessor::testFaceVsPlane Date: 19.01.00 *
792  * Author: E.Chernyaev Revised: *
793  * *
794  * Function: Find result of intersection of face by plane *
795  * *
796  ***********************************************************************/
797 {
798  int iface = edge.iface1;
799  const HVPlane3D& plane = faces[edge.iface2].plane;
800  int i, nnode, npos = 0, nneg = 0, nzer = 0;
801  double dd[5];
802 
803  // F I N D D I S T A N C E S
804 
805  nnode = (faces[iface].iedges[3] == 0) ? 3 : 4;
806  for (i=0; i<nnode; i++) {
807  dd[i] = plane.distance(nodes[edges[faces[iface].iedges[i]].i1].v);
808  if (dd[i] > del) {
809  npos++;
810  }else if (dd[i] < -del) {
811  nneg++;
812  }else{
813  nzer++; dd[i] = 0;
814  }
815  }
816 
817  // 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 )
818 
819  if (npos == nnode || nneg == nnode) return OUT_OF_PLANE;
820  if (nzer == 1 && nneg == 0) return OUT_OF_PLANE;
821  if (nzer == 1 && npos == 0) return OUT_OF_PLANE;
822  if (nzer == nnode) return ON_PLANE;
823  if (nzer == 3) return NON_PLANAR_FACE;
824 
825  // F I N D I N T E R S E C T I O N
826 
827  int ie1 = 0, ie2 = 0, s1 = 0, s2 = 0, status, nint = 0;
828  enum { PLUS_MINUS, MINUS_PLUS, ZERO_ZERO, ZERO_PLUS, ZERO_MINUS };
829 
830  dd[nnode] = dd[0];
831  for (i=0; i<nnode; i++) {
832  if (dd[i] > 0) {
833  if (dd[i+1] >= 0) continue;
834  status = PLUS_MINUS;
835  }else if (dd[i] < 0) {
836  if (dd[i+1] <= 0) continue;
837  status = MINUS_PLUS;
838  }else{
839  status = ZERO_ZERO;
840  if (dd[i+1] > 0) status = ZERO_PLUS;
841  if (dd[i+1] < 0) status = ZERO_MINUS;
842  }
843  switch (nint) {
844  case 0:
845  ie1 = i; s1 = status; nint++; break;
846  case 1:
847  ie2 = i; s2 = status; nint++; break;
848  default:
849  return NON_PLANAR_FACE;
850  }
851  }
852  if (nint != 2) return NON_PLANAR_FACE;
853 
854  // F O R M I N T E R S E C T I O N S E G M E N T
855 
856  if (s1 != ZERO_ZERO && s2 != ZERO_ZERO) {
857  if (s1 == s2) return NON_PLANAR_FACE;
858  int iedge, i1 = 0, i2 = 0, ii[2];
859  double d1 = 0., d2 = 0., ddd ;
860  ii[0] = ie1; ii[1] = ie2;
861  for (i=0; i<2; i++) {
862  iedge = faces[iface].iedges[ii[i]];
863  while (iedge > 0) {
864  i1 = edges[iedge].i1;
865  i2 = edges[iedge].i2;
866 
867  d1 = plane.distance(nodes[i1].v);
868  d2 = plane.distance(nodes[i2].v);
869  if (d1 > del) {
870  if (d2 < -del) { ii[i] = int(nodes.size()); break; } // +-
871  }else if (d1 < -del) {
872  if (d2 > del) { ii[i] = int(nodes.size()); break; } // -+
873  }else{
874  ii[i] = i1; break; // 0+ or 0-
875  }
876  iedge = edges[iedge].inext;
877  }
878  if (ii[i] == (int)nodes.size()) {
879  ddd = d2-d1; d1 = d1/ddd; d2 = d2/ddd;
880  nodes.push_back(ExtNode(d2*nodes[i1].v-d1*nodes[i2].v, iedge));
881  }
882  }
883  edge.inext = 0;
884  if (s1 == MINUS_PLUS || s1 == ZERO_PLUS) {
885  edge.i1 = ii[1];
886  edge.i2 = ii[0];
887  }else{
888  edge.i1 = ii[0];
889  edge.i2 = ii[1];
890  }
891  return INTERSECTION;
892  }else{
893  if (npos == nneg) return NON_PLANAR_FACE;
894  edge.inext = (s1 == ZERO_ZERO) ? ie1+1 : ie2+1;
895  if (s1 == ZERO_PLUS || s2 == ZERO_MINUS) {
896  edge.i1 = edges[faces[iface].iedges[ie2]].i1;
897  edge.i2 = edges[faces[iface].iedges[ie1]].i1;
898  }else{
899  edge.i1 = edges[faces[iface].iedges[ie1]].i1;
900  edge.i2 = edges[faces[iface].iedges[ie2]].i1;
901  }
902  return EDGE;
903  }
904 }
905 
906 inline
907 void BooleanProcessor::renumberNodes(int & i1, int & i2, int & i3, int & i4)
908 /***********************************************************************
909  * *
910  * Name: BooleanProcessor::renumberNodes Date: 19.01.00 *
911  * Author: E.Chernyaev Revised: *
912  * *
913  * Function: Renumber nodes and remove last temporary node. *
914  * Remark: In principal this routine can be replaced just *
915  * with i1 = i2; *
916  * Removal of temporary nodes provides additional control *
917  * on number of nodes, that is very useful for debugging. *
918  * *
919  ***********************************************************************/
920 {
921  if (i1 == i2) return;
922  if (nodes[i1].s == 0 || nodes.back().s == 0) { i1 = i2; return; }
923 
924  int ilast = int(nodes.size())-1;
925  if (i1 == ilast) { i1 = i2; nodes.pop_back(); return; }
926  if (i2 == ilast) { i2 = i1; }
927  if (i3 == ilast) { i3 = i1; }
928  if (i4 == ilast) { i4 = i1; }
929  nodes[i1] = nodes.back(); i1 = i2; nodes.pop_back();
930 }
931 
932 inline
933 int BooleanProcessor::testEdgeVsEdge(ExtEdge & edge1, ExtEdge & edge2)
934 /***********************************************************************
935  * *
936  * Name: BooleanProcessor::testEdgeVsEdge Date: 19.01.00 *
937  * Author: E.Chernyaev Revised: *
938  * *
939  * Function: Find common part of two edges *
940  * *
941  ***********************************************************************/
942 {
943  int i, ii = 0;
944  double d, dd = 0.;
945 
946  for (i=0; i<3; i++) {
947  d = nodes[edge1.i1].v[i]-nodes[edge1.i2].v[i];
948  if (d < 0.) d = -d;
949  if (d > dd) { dd = d; ii = i; }
950  }
951  double t1 = nodes[edge1.i1].v[ii];
952  double t2 = nodes[edge1.i2].v[ii];
953  double t3 = nodes[edge2.i1].v[ii];
954  double t4 = nodes[edge2.i2].v[ii];
955  if (t2-t1 < 0.) { t1 = -t1; t2 = -t2; t3 = -t3; t4 = -t4; }
956 
957  if (t3 <= t1+del || t4 >= t2-del) return 0;
958  if (t3 > t2+del) {
959  renumberNodes(edge2.i1, edge1.i2, edge1.i1, edge2.i2);
960  }else if (t3 < t2-del) {
961  renumberNodes(edge1.i2, edge2.i1, edge1.i1, edge2.i2);
962  }
963  if (t4 < t1-del) {
964  renumberNodes(edge2.i2, edge1.i1, edge1.i2, edge2.i1);
965  }else if (t4 > t1+del) {
966  renumberNodes(edge1.i1, edge2.i2, edge1.i2, edge2.i1);
967  }
968  return 1;
969 }
970 
971 inline
972 void BooleanProcessor::divideEdge(int & i1, int & i2)
973 /***********************************************************************
974  * *
975  * Name: BooleanProcessor::divideEdge Date: 24.01.00 *
976  * Author: E.Chernyaev Revised: *
977  * *
978  * Function: Unify the nodes and divide edge on two parts by the node. *
979  * *
980  ***********************************************************************/
981 {
982  int iedges[2];
983  iedges[0] = nodes[i1].s;
984  iedges[1] = nodes[i2].s;
985 
986  // U N I F Y N O D E S
987 
988  if (i1 < i2) { i2 = i1; }
989  else if (i1 > i2) { i1 = i2; }
990  else { iedges[1] = 0; }
991  if (iedges[0] == iedges[1]) return;
992 
993  int ie1, ie2, inode = i1;
994  nodes[inode].s = 0;
995  for (int i=0; i<2; i++) {
996 
997  // F I N D C O R R E S P O N D I N G E D G E
998 
999  if ((ie1 = iedges[i]) == 0) continue;
1000  ie2 = faces[edges[ie1].iface2].iedges[0];
1001  while (ie2 > 0) {
1002  if (edges[ie2].i1 == edges[ie1].i2 &&
1003  edges[ie2].i2 == edges[ie1].i1) break;
1004  ie2 = edges[ie2].inext;
1005  }
1006 
1007  // D I V I D E E D G E S
1008 
1009  edges.push_back(edges[ie1]);
1010  edges[ie1].inext = int(edges.size()) - 1;
1011  edges[ie1].i2 = inode;
1012  edges.back().i1 = inode;
1013 
1014  edges.push_back(edges[ie2]);
1015  edges[ie2].inext = int(edges.size()) - 1;
1016  edges[ie2].i2 = inode;
1017  edges.back().i1 = inode;
1018  }
1019 }
1020 
1021 inline
1022 void BooleanProcessor::insertEdge(const ExtEdge & edge)
1023 /***********************************************************************
1024  * *
1025  * Name: BooleanProcessor::insertEdge Date: 24.01.00 *
1026  * Author: E.Chernyaev Revised: *
1027  * *
1028  * Function: Insert edge to the list of new edges *
1029  * *
1030  ***********************************************************************/
1031 {
1032  int iface = edge.iface1;
1033  edges.push_back(edge);
1034  edges.back().inext = faces[iface].inew;
1035  faces[iface].inew = int(edges.size()) - 1;
1036 }
1037 
1038 inline
1039 void BooleanProcessor::caseII(ExtEdge & edge1, ExtEdge & edge2)
1040 /***********************************************************************
1041  * *
1042  * Name: BooleanProcessor::caseII Date: 19.01.00 *
1043  * Author: E.Chernyaev Revised: *
1044  * *
1045  * Function: Intersection/Intersection case *
1046  * *
1047  ***********************************************************************/
1048 {
1049  divideEdge(edge1.i1, edge2.i2);
1050  divideEdge(edge1.i2, edge2.i1);
1051  edge1.ivis = 1;
1052  edge2.ivis = 1;
1053  insertEdge(edge1);
1054  insertEdge(edge2);
1055 }
1056 
1057 inline
1058 bool BooleanProcessor::caseIE(ExtEdge &, ExtEdge &)
1059 /***********************************************************************
1060  * *
1061  * Name: BooleanProcessor::caseIE Date: 19.01.00 *
1062  * Author: E.Chernyaev Revised: *
1063  * *
1064  * Function: Intersection/Edge-touch case *
1065  * *
1066  ***********************************************************************/
1067 {
1068  //G.Barrand : processor_error = 1;
1069 #ifdef TOOLS_HEP_BP_OUT_ERR
1070  std::cerr
1071  << "BooleanProcessor::caseIE : unimplemented case"
1072  << std::endl;
1073 #endif
1074  return false; //G.Barrand
1075 }
1076 
1077 inline
1078 bool BooleanProcessor::caseEE(ExtEdge &, ExtEdge &)
1079 /***********************************************************************
1080  * *
1081  * Name: BooleanProcessor::caseEE Date: 19.01.00 *
1082  * Author: E.Chernyaev Revised: *
1083  * *
1084  * Function: Edge-touch/Edge-touch case *
1085  * *
1086  ***********************************************************************/
1087 {
1088  //G.Barrand : processor_error = 1;
1089 #ifdef TOOLS_HEP_BP_OUT_ERR
1090  std::cerr
1091  << "BooleanProcessor::caseEE : unimplemented case"
1092  << std::endl;
1093 #endif
1094  return false;
1095 }
1096 
1097 inline
1098 bool BooleanProcessor::testFaceVsFace(int iface1, int iface2)
1099 /***********************************************************************
1100  * *
1101  * Name: BooleanProcessor::testFaceVsFace Date: 11.01.00 *
1102  * Author: E.Chernyaev Revised: *
1103  * *
1104  * Function: Find result (an edge) of intersection of two faces *
1105  * *
1106  ***********************************************************************/
1107 {
1108  ExtEdge edge1, edge2;
1109  int irep1, irep2;
1110 
1111  // M I N - M A X
1112 
1113  {const ExtFace& face_1 = faces[iface1]; //G.Barrand : optimize
1114  const ExtFace& face_2 = faces[iface2];
1115  for (int i=0; i<3; i++) {
1116  if (face_1.rmin[i] > face_2.rmax[i] + del) return true; //GB : true ?
1117  if (face_1.rmax[i] < face_2.rmin[i] - del) return true; //GB : true ?
1118  }}
1119 
1120  // F A C E - 1 vs P L A N E - 2
1121 
1122  edge1.iface1 = iface1;
1123  edge1.iface2 = iface2;
1124  irep1 = testFaceVsPlane(edge1);
1125  if (irep1 == OUT_OF_PLANE || irep1 == ON_PLANE) {
1126  removeJunkNodes();
1127  return true; //GB : true ?
1128  }
1129 
1130  // F A C E - 2 vs P L A N E - 1
1131 
1132  edge2.iface1 = iface2;
1133  edge2.iface2 = iface1;
1134  irep2 = testFaceVsPlane(edge2);
1135  if (irep2 == OUT_OF_PLANE || irep2 == ON_PLANE) {
1136  removeJunkNodes();
1137  return true; //GB : true ?
1138  }
1139 
1140  // C H E C K F O R N O N P L A N A R F A C E
1141 
1142  if (irep1 == NON_PLANAR_FACE || irep2 == NON_PLANAR_FACE) {
1143  removeJunkNodes();
1144  return true; //GB : true ?
1145  }
1146 
1147  // F I N D I N T E R S E C T I O N P A R T
1148 
1149  if (testEdgeVsEdge(edge1, edge2) == 0) {
1150  return true; //GB : true ?
1151  }
1152 
1153  // C O N S I D E R D I F F E R E N T C A S E S
1154 
1155  if (irep1 == INTERSECTION && irep2 == INTERSECTION) caseII(edge1, edge2);
1156  if (irep1 == INTERSECTION && irep2 == EDGE) {
1157  if(!caseIE(edge1, edge2)) {
1158  //G.Barrand : processor_error = 1;
1159  removeJunkNodes(); //G.Barrand
1160  return false; //G.Barrand
1161  }
1162  }
1163  if (irep1 == EDGE && irep2 == INTERSECTION) {
1164  if(!caseIE(edge2, edge1)) {
1165  //G.Barrand : processor_error = 1;
1166  removeJunkNodes(); //G.Barrand
1167  return false; //G.Barrand
1168  }
1169  }
1170  if (irep1 == EDGE && irep2 == EDGE) {
1171  if(!caseEE(edge1, edge2)) {
1172  //G.Barrand : processor_error = 1;
1173  removeJunkNodes(); //G.Barrand
1174  return false; //G.Barrand
1175  }
1176  }
1177  removeJunkNodes();
1178 
1179  return true; //G.Barrand
1180 }
1181 
1182 inline
1183 void BooleanProcessor::invertNewEdges(int iface)
1184 /***********************************************************************
1185  * *
1186  * Name: BooleanProcessor::invertNewEdges Date: 04.02.00 *
1187  * Author: E.Chernyaev Revised: *
1188  * *
1189  * Function: Invert direction of new edges *
1190  * *
1191  ***********************************************************************/
1192 {
1193  int iedge = faces[iface].inew;
1194  while (iedge > 0) {
1195  edges[iedge].invert();
1196  iedge = edges[iedge].inext;
1197  }
1198 }
1199 
1200 inline
1201 void BooleanProcessor::checkDoubleEdges(int)
1202 /***********************************************************************
1203  * *
1204  * Name: BooleanProcessor::checkDoubleEdges Date: 04.02.00 *
1205  * Author: E.Chernyaev Revised: *
1206  * *
1207  * Function: Eliminate duplication of edges *
1208  * *
1209  ***********************************************************************/
1210 {
1211 
1212 }
1213 
1214 inline
1215 bool BooleanProcessor::assembleFace(int what, int iface)
1216 /***********************************************************************
1217  * *
1218  * Name: BooleanProcessor::assembleFace Date: 19.02.00 *
1219  * Author: E.Chernyaev Revised: *
1220  * *
1221  * Function: Assemble face *
1222  * *
1223  ***********************************************************************/
1224 {
1225  // A S S E M B L E N E W F A C E
1226 
1227  int ihead; // head of the list of edges for new face
1228  int icur; // current edge in the list - last edge inserted to the list
1229  int *ilink; // pointer to the current link
1230  int ifirst; // first node of a contour
1231  int *i; // pointer to the index of the current edge in a loop
1232  int ioldflag=0; // is set if an edge from iold has been taken
1233 
1234  ExtFace& face = faces[iface]; //G.Barrand : optimize.
1235  ilink = &ihead;
1236  for(;;) {
1237  if (face.inew == 0) break;
1238 
1239  // S T A R T N E W C O N T O U R
1240 
1241  icur = face.inew;
1242  face.inew = edges[icur].inext;
1243 
1244  //INSERT_EDGE_TO_THE_LIST(icur); //G.Barrand
1245  *ilink = icur; ilink = &edges[icur].inext; *ilink = 0; //G.Barrand
1246 
1247  ifirst = edges[icur].i1;
1248 
1249  // C O N S T R U C T T H E C O N T O U R
1250 
1251  for (;;) {
1252  i = &face.inew;
1253  ExtEdge& edge_cur = edges[icur];
1254  while(*i > 0) {
1255  ExtEdge& edge_i = edges[*i];
1256  if (edge_i.i1 == edge_cur.i2) break;
1257  i = &edge_i.inext;
1258  }
1259  if (*i == 0) {
1260  i = &face.iold;
1261  while(*i > 0) {
1262  ExtEdge& edge_i = edges[*i];
1263  if (edge_i.i1 == edge_cur.i2) {
1264  ioldflag = 1;
1265  break;
1266  }
1267  i = &edge_i.inext;
1268  }
1269  }
1270  if (*i > 0) {
1271  icur = *i;
1272  *i = edges[icur].inext;
1273 
1274  //INSERT_EDGE_TO_THE_LIST(icur);
1275  *ilink = icur; ilink = &edges[icur].inext; *ilink = 0; //G.Barrand
1276 
1277  if (edges[icur].i2 == ifirst) { break; } else { continue; }
1278  }else{
1279  //G.Barrand : processor_error = 1;
1280 #ifdef TOOLS_HEP_BP_OUT_ERR
1281  std::cerr
1282  << "BooleanProcessor::assembleFace(" << iface << ") : "
1283  << "could not find next edge of the contour"
1284  << std::endl;
1285 #endif
1286  face.inew = DEFECTIVE_FACE;
1287  return false; //G.Barrand : false
1288  }
1289  }
1290  }
1291 
1292  // C H E C K O R I G I N A L C O N T O U R
1293 
1294  int iedge;
1295  iedge = face.iold;
1296  if (what == 0 && ioldflag == 0 && iedge > 0) {
1297  for (;;) {
1298  if (edges[iedge].inext > 0) {
1299  if (edges[iedge].i2 == edges[edges[iedge].inext].i1) {
1300  iedge = edges[iedge].inext;
1301  }else{
1302  break;
1303  }
1304  }else{
1305  if (edges[iedge].i2 == edges[face.iold].i1) {
1306  edges[iedge].inext = ihead; // set new face
1307  return true; //G.Barrand : true.
1308  }else{
1309  break;
1310  }
1311  }
1312  }
1313  }
1314 
1315  // 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
1316 
1317  int iface2;
1318  iedge = face.iold;
1319  while(iedge > 0) {
1320  iface2 = edges[iedge].iface2;
1321  if (faces[iface2].inew == 0) faces[iface2].inew = UNSUITABLE_FACE;
1322  iedge = edges[iedge].inext;
1323  }
1324  face.iold = ihead; // set new face
1325 
1326  return true;
1327 }
1328 
1329 inline
1330 bool BooleanProcessor::assembleNewFaces(int what, int ihead)
1331 /***********************************************************************
1332  * *
1333  * Name: BooleanProcessor::assembleNewFaces Date: 30.01.00 *
1334  * Author: E.Chernyaev Revised: *
1335  * *
1336  * Function: Assemble internal or external parts of faces *
1337  * *
1338  ***********************************************************************/
1339 {
1340  int iface = ihead;
1341  while(iface > 0) {
1342  if (faces[iface].inew > 0) {
1343  if (what != 0) invertNewEdges(iface);
1344  checkDoubleEdges(iface);
1345  if(!assembleFace(what, iface)) return false; //G.Barrand : test.
1346  faces[iface].inew =
1347  (faces[iface].iold == 0) ? UNSUITABLE_FACE : NEW_FACE;
1348  }
1349  iface = faces[iface].inext;
1350  }
1351  return true; //G.Barrand
1352 }
1353 
1354 inline
1355 bool BooleanProcessor::initiateLists() //G.Barrand : return bool.
1356 /***********************************************************************
1357  * *
1358  * Name: BooleanProcessor::initiateLists Date: 28.02.00 *
1359  * Author: E.Chernyaev Revised: *
1360  * *
1361  * Function: Initiate lists of faces. *
1362  * *
1363  ***********************************************************************/
1364 {
1365  int i, iface;
1366 
1367  // R E S E T L I S T S O F F A C E S
1368 
1369  result_faces.clean();
1370  suitable_faces.clean();
1371  unsuitable_faces.clean();
1372  unknown_faces.clean();
1373 
1374  // I N I T I A T E T H E L I S T S
1375 
1376  iface = iout1;
1377  while (iface > 0) {
1378  i = iface;
1379 #ifdef TOOLS_HEP_BP_CHECK_INDEX
1380  if(!check_face_index("initiateLists(1)",i)) return false;
1381 #endif
1382  iface = faces[i].inext;
1383  if (operation == OP_INTERSECTION) {
1384  unsuitable_faces.push_back(i);
1385  faces[i].inew = UNSUITABLE_FACE;
1386  }else{
1387  suitable_faces.push_back(i);
1388  faces[i].inew = ORIGINAL_FACE;
1389  }
1390  }
1391 
1392  iface = iout2;
1393  while (iface > 0) {
1394  i = iface;
1395 #ifdef TOOLS_HEP_BP_CHECK_INDEX
1396  if(!check_face_index("initiateLists(2)",i)) return false;
1397 #endif
1398  iface = faces[i].inext;
1399  if (operation == OP_UNION) {
1400  suitable_faces.push_back(i);
1401  faces[i].inew = ORIGINAL_FACE;
1402  }else{
1403  unsuitable_faces.push_back(i);
1404  faces[i].inew = UNSUITABLE_FACE;
1405  }
1406  }
1407 
1408  iface = iunk1;
1409  while (iface > 0) {
1410  i = iface;
1411 #ifdef TOOLS_HEP_BP_CHECK_INDEX
1412  if(!check_face_index("initiateLists(3)",i)) return false;
1413 #endif
1414  iface = faces[i].inext;
1415  unknown_faces.push_back(i);
1416  }
1417  iface = iunk2;
1418  while (iface > 0) {
1419  i = iface;
1420 #ifdef TOOLS_HEP_BP_CHECK_INDEX
1421  if(!check_face_index("initiateLists(4)",i)) return false;
1422 #endif
1423  iface = faces[i].inext;
1424  if (operation == OP_SUBTRACTION) {
1425  if(!faces[i].invert()) return false; //G.Barrand
1426  }
1427  unknown_faces.push_back(i);
1428  }
1429 
1430  iface = ifaces1;
1431  while (iface > 0) {
1432  i = iface;
1433 #ifdef TOOLS_HEP_BP_CHECK_INDEX
1434  if(!check_face_index("initiateLists(5)",i)) return false;
1435 #endif
1436  iface = faces[i].inext;
1437  switch(faces[i].inew) {
1438  case UNKNOWN_FACE:
1439  unknown_faces.push_back(i);
1440  break;
1441  case ORIGINAL_FACE: case NEW_FACE:
1442  suitable_faces.push_back(i);
1443  break;
1444  case UNSUITABLE_FACE:
1445  unsuitable_faces.push_back(i);
1446  break;
1447  default:
1448  faces[i].iprev = 0;
1449  faces[i].inext = 0;
1450  break;
1451  }
1452  }
1453 
1454  iface = ifaces2;
1455  while (iface > 0) {
1456  i = iface;
1457  iface = faces[i].inext;
1458 #ifdef TOOLS_HEP_BP_CHECK_INDEX
1459  if(!check_face_index("initiateLists(6)",i)) return false;
1460 #endif
1461 
1462  if (operation == OP_SUBTRACTION) {
1463  if(!faces[i].invert()) return false;
1464  }
1465 
1466  switch(faces[i].inew) {
1467  case UNKNOWN_FACE:
1468  unknown_faces.push_back(i);
1469  break;
1470  case ORIGINAL_FACE: case NEW_FACE:
1471  suitable_faces.push_back(i);
1472  break;
1473  case UNSUITABLE_FACE:
1474  unsuitable_faces.push_back(i);
1475  break;
1476  default:
1477  faces[i].iprev = 0;
1478  faces[i].inext = 0;
1479  break;
1480  }
1481  }
1482  ifaces1 = ifaces2 = iout1 = iout2 = iunk1 = iunk2 = 0;
1483  return true; //G.Barrand.
1484 }
1485 
1486 inline
1487 bool BooleanProcessor::assemblePolyhedra() //G.Barrand : return bool.
1488 /***********************************************************************
1489  * *
1490  * Name: BooleanProcessor::assemblePolyhedra() Date: 10.12.99 *
1491  * Author: E.Chernyaev Revised: *
1492  * *
1493  * Function: Collect suitable faces and remove unsuitable ones. *
1494  * *
1495  ***********************************************************************/
1496 {
1497  int i, iedge, iface;
1498 
1499  // L O O P A L O N G S U I T A B L E F A C E S
1500 
1501  iface = suitable_faces.front();
1502  while(iface > 0) {
1503  i = iface;
1504 #ifdef TOOLS_HEP_BP_CHECK_INDEX
1505  if(!check_face_index("assemblePolyhedra(1)",i)) return false;
1506 #endif
1507  iedge = faces[i].iold;
1508  while(iedge > 0) {
1509 #ifdef TOOLS_HEP_BP_CHECK_INDEX
1510  if(!check_edge_index("assemblePolyhedra(2)",iedge)) return false;
1511 #endif
1512  iface = edges[iedge].iface2;
1513 #ifdef TOOLS_HEP_BP_CHECK_INDEX
1514  if(!check_face_index("assemblePolyhedra(3)",iface)) return false;
1515 #endif
1516  if (faces[iface].inew == UNKNOWN_FACE) {
1517  unknown_faces.remove(iface);
1518  suitable_faces.push_back(iface);
1519  faces[iface].inew = ORIGINAL_FACE;
1520  }
1521  iedge = edges[iedge].inext;
1522  }
1523 #ifdef TOOLS_HEP_BP_CHECK_INDEX
1524  if(!check_face_index("assemblePolyhedra(4)",i)) return false;
1525 #endif
1526  iface = faces[i].inext;
1527  suitable_faces.remove(i);
1528  result_faces.push_back(i);
1529  }
1530  if (unknown_faces.front() == 0) return true;
1531 
1532  // L O O P A L O N G U N S U I T A B L E F A C E S
1533 
1534  iface = unsuitable_faces.front();
1535  while(iface > 0) {
1536  i = iface;
1537 #ifdef TOOLS_HEP_BP_CHECK_INDEX
1538  if(!check_face_index("assemblePolyhedra(5)",i)) return false;
1539 #endif
1540  iedge = faces[i].iold;
1541  while(iedge > 0) {
1542 #ifdef TOOLS_HEP_BP_CHECK_INDEX
1543  if(!check_edge_index("assemblePolyhedra(6)",iedge)) return false;
1544 #endif
1545  iface = edges[iedge].iface2;
1546 #ifdef TOOLS_HEP_BP_CHECK_INDEX
1547  if(!check_face_index("assemblePolyhedra(7)",iface)) return false;
1548 #endif
1549  if (faces[iface].inew == UNKNOWN_FACE) {
1550  unknown_faces.remove(iface);
1551  unsuitable_faces.push_back(iface);
1552  faces[iface].inew = UNSUITABLE_FACE;
1553  }
1554  iedge = edges[iedge].inext;
1555  }
1556 #ifdef TOOLS_HEP_BP_CHECK_INDEX
1557  if(!check_face_index("assemblePolyhedra(8)",i)) return false;
1558 #endif
1559  iface = faces[i].inext;
1560  unsuitable_faces.remove(i);
1561  }
1562 
1563  //G.Barrand : begin
1564  /* From S.Ponce
1565  At last, there is a problem in the assemblePolyhedra method. At least, I
1566  think it is there. The problem deals with boolean operations on solids,
1567  when one of the two contains entirely the other one. It has no sense for
1568  intersection and union but still has sense for subtraction. In this
1569  case, faces from the inner solid are stored in the unknown_faces
1570  FaceList. And an error occurs in the execute method. This may be correct
1571  for intersection and union but in the case of subtraction, one should do
1572  that in assemblePolyhedra :
1573  */
1574  // Unknown faces are actually suitable face !!!
1575  iface = unknown_faces.front();
1576  while(iface > 0) {
1577  i = iface;
1578 #ifdef TOOLS_HEP_BP_CHECK_INDEX
1579  if(!check_face_index("assemblePolyhedra(9)",i)) return false;
1580 #endif
1581  faces[i].inew = ORIGINAL_FACE;
1582  iface = faces[i].inext;
1583  unknown_faces.remove(i);
1584  result_faces.push_back(i);
1585  }
1586  /*
1587  Otherwise, the inner hole that the second solid was building in the
1588  first one does not exist. I'm not very clear on what to do for unions
1589  and intersections. I think this kind of situation should be detected and
1590  one of the solid should simply be ignored.
1591  */
1592  //G.Barrand : end
1593 
1594  return true;
1595 }
1596 
1597 inline void
1598 BooleanProcessor::findABC(double x1, double y1, double x2, double y2,
1599  double &a, double &b, double &c) const
1600 /***********************************************************************
1601  * *
1602  * Name: BooleanProcessor::findABC Date: 07.03.00 *
1603  * Author: E.Chernyaev Revised: *
1604  * *
1605  * Function: Find line equation Ax+By+C=0 *
1606  * *
1607  ***********************************************************************/
1608 {
1609  double w;
1610  a = y1 - y2;
1611  b = x2 - x1;
1612  //G.Barrand : w = std::abs(a)+std::abs(b);
1613  w = ::fabs(a)+::fabs(b); //G.Barrand
1614  a /= w;
1615  b /= w;
1616  c = -(a*x2 + b*y2);
1617 }
1618 
1619 inline
1620 int BooleanProcessor::checkDirection(double *x, double *y) const
1621 /***********************************************************************
1622  * *
1623  * Name: BooleanProcessor::checkDirection Date: 06.03.00 *
1624  * Author: E.Chernyaev Revised: *
1625  * *
1626  * Function: Check direction of line 1-4 *
1627  * *
1628  ***********************************************************************/
1629 {
1630  double a1, b1, c1, a2, b2, c2, d1, d2;
1631 
1632  // 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
1633 
1634  findABC(x[0], y[0], x[1], y[1], a1, b1, c1);
1635  findABC(x[1], y[1], x[2], y[2], a2, b2, c2);
1636  d1 = a1*x[4] + b1*y[4] + c1;
1637  d2 = a2*x[4] + b2*y[4] + c2;
1638  if (d1 <= del && d2 <= del) return 1;
1639  if (! (d1 > del && d2 > del)) {
1640  if ( a1*x[2] + b1*y[2] + c1 >= -del) return 1;
1641  }
1642 
1643  // 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
1644 
1645  findABC(x[3], y[3], x[4], y[4], a1, b1, c1);
1646  findABC(x[4], y[4], x[5], y[5], a2, b2, c2);
1647  d1 = a1*x[1] + b1*y[1] + c1;
1648  d2 = a2*x[1] + b2*y[1] + c2;
1649  if (d1 <= del && d2 <= del) return 1;
1650  if (!(d1 > del && d2 > del)) {
1651  if ( a1*x[5] + b1*y[5] + c1 >= -del) return 1;
1652  }
1653  return 0;
1654 }
1655 
1656 inline
1657 int BooleanProcessor::checkIntersection(int ix, int iy, int i1, int i2) const
1658 /***********************************************************************
1659  * *
1660  * Name: BooleanProcessor::checkDirection Date: 06.03.00 *
1661  * Author: E.Chernyaev Revised: *
1662  * *
1663  * Function: Check line i1-i2 on intersection with contours *
1664  * *
1665  ***********************************************************************/
1666 {
1667  // F I N D L I N E E Q U A T I O N
1668 
1669  double x1, y1, x2, y2, a1, b1, c1;
1670  x1 = nodes[i1].v[ix];
1671  y1 = nodes[i1].v[iy];
1672  x2 = nodes[i2].v[ix];
1673  y2 = nodes[i2].v[iy];
1674  findABC(x1, y1, x2, y2, a1, b1, c1);
1675 
1676  // 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
1677 
1678  int icontour, iedge, k1, k2;
1679  double x3, y3, x4, y4, a2, b2, c2, d1, d2;
1680  for(icontour=0; icontour<(int)external_contours.size(); icontour++) {
1681  iedge = external_contours[icontour];
1682  while(iedge > 0) {
1683  k1 = edges[iedge].i1;
1684  k2 = edges[iedge].i2;
1685  iedge = edges[iedge].inext;
1686  if (k1 == i1 || k2 == i1) continue;
1687  if (k1 == i2 || k2 == i2) continue;
1688  x3 = nodes[k1].v[ix];
1689  y3 = nodes[k1].v[iy];
1690  x4 = nodes[k2].v[ix];
1691  y4 = nodes[k2].v[iy];
1692  d1 = a1*x3 + b1*y3 + c1;
1693  d2 = a1*x4 + b1*y4 + c1;
1694  if (d1 > del && d2 > del) continue;
1695  if (d1 < -del && d2 < -del) continue;
1696 
1697  findABC(x3, y3, x4, y4, a2, b2, c2);
1698  d1 = a2*x1 + b2*y1 + c2;
1699  d2 = a2*x2 + b2*y2 + c2;
1700  if (d1 > del && d2 > del) continue;
1701  if (d1 < -del && d2 < -del) continue;
1702  return 1;
1703  }
1704  }
1705 
1706  // 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
1707 
1708  for(icontour=0; icontour<(int)internal_contours.size(); icontour++) {
1709  iedge = internal_contours[icontour];
1710  while(iedge > 0) {
1711  k1 = edges[iedge].i1;
1712  k2 = edges[iedge].i2;
1713  iedge = edges[iedge].inext;
1714  if (k1 == i1 || k2 == i1) continue;
1715  if (k1 == i2 || k2 == i2) continue;
1716  x3 = nodes[k1].v[ix];
1717  y3 = nodes[k1].v[iy];
1718  x4 = nodes[k2].v[ix];
1719  y4 = nodes[k2].v[iy];
1720  d1 = a1*x3 + b1*y3 + c1;
1721  d2 = a1*x4 + b1*y4 + c1;
1722  if (d1 > del && d2 > del) continue;
1723  if (d1 < -del && d2 < -del) continue;
1724 
1725  findABC(x3, y3, x4, y4, a2, b2, c2);
1726  d1 = a2*x1 + b2*y1 + c2;
1727  d2 = a2*x2 + b2*y2 + c2;
1728  if (d1 > del && d2 > del) continue;
1729  if (d1 < -del && d2 < -del) continue;
1730  return 1;
1731  }
1732  }
1733  return 0;
1734 }
1735 
1736 inline
1737 void BooleanProcessor::mergeContours(int ix, int iy, int kext, int kint)
1738 /***********************************************************************
1739  * *
1740  * Name: BooleanProcessor::mergeContours Date: 06.03.00 *
1741  * Author: E.Chernyaev Revised: *
1742  * *
1743  * Function: Attemp to merge internal contour with external one *
1744  * *
1745  ***********************************************************************/
1746 {
1747  int i1ext, i2ext, i1int, i2int, i, k[6];
1748  double x[6], y[6];
1749 
1750  // L O O P A L O N G E X T E R N A L C O N T O U R
1751 
1752  i1ext = external_contours[kext];
1753  while (i1ext > 0) {
1754  i2ext = edges[i1ext].inext;
1755  if (i2ext == 0) i2ext = external_contours[kext];
1756  k[0] = edges[i1ext].i1;
1757  k[1] = edges[i1ext].i2;
1758  k[2] = edges[i2ext].i2;
1759  for (i=0; i<3; i++) {
1760  x[i] = nodes[k[i]].v[ix];
1761  y[i] = nodes[k[i]].v[iy];
1762  }
1763 
1764  // L O O P A L O N G I N T E R N A L C O N T O U R
1765 
1766  i1int = internal_contours[kint];
1767  while (i1int > 0) {
1768  i2int = edges[i1int].inext;
1769  if (i2int == 0) i2int = internal_contours[kint];
1770  k[3] = edges[i1int].i1;
1771  k[4] = edges[i1int].i2;
1772  k[5] = edges[i2int].i2;
1773  for (i=3; i<6; i++) {
1774  x[i] = nodes[k[i]].v[ix];
1775  y[i] = nodes[k[i]].v[iy];
1776  }
1777 
1778  // T E S T L I N E K1 - K4
1779  // I F O K T H E N M E R G E C O N T O U R S
1780 
1781  if (checkDirection(x, y) == 0) {
1782  if (checkIntersection(ix, iy, k[1], k[4]) == 0) {
1783  i = i1int;
1784  for(;;) {
1785  if (edges[i].inext == 0) {
1786  edges[i].inext = internal_contours[kint];
1787  internal_contours[kint] = 0;
1788  break;
1789  }else{
1790  i = edges[i].inext;
1791  }
1792  }
1793  i = edges[i1int].iface1;
1794  edges.push_back(ExtEdge(k[1], k[4], i, -1*(int(edges.size())+1), -1));
1795  edges.back().inext = i2int;
1796  edges.push_back(ExtEdge(k[4], k[1], i, -1*(int(edges.size())-1), -1));
1797  edges.back().inext = edges[i1ext].inext;
1798  edges[i1ext].inext = int(edges.size())-2;
1799  edges[i1int].inext = int(edges.size())-1;
1800  return;
1801  }
1802  }
1803  i1int = edges[i1int].inext;
1804  }
1805  i1ext = edges[i1ext].inext;
1806  }
1807 }
1808 
1809 inline bool
1810 BooleanProcessor::checkTriangle(int iedge1, int iedge2, int ix, int iy) const
1811 /***********************************************************************
1812  * *
1813  * Name: BooleanProcessor::checkTriangle Date: 08.03.00 *
1814  * Author: E.Chernyaev Revised: *
1815  * *
1816  * Function: Check triangle for correctness *
1817  * *
1818  ***********************************************************************/
1819 {
1820  int k[3];
1821  double x[3], y[3];
1822  double a1, b1, c1;
1823 
1824  k[0] = edges[iedge1].i1;
1825  k[1] = edges[iedge1].i2;
1826  k[2] = edges[iedge2].i2;
1827  for (int i=0; i<3; i++) {
1828  x[i] = nodes[k[i]].v[ix];
1829  y[i] = nodes[k[i]].v[iy];
1830  }
1831 
1832  // 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
1833 
1834  findABC(x[2], y[2], x[0], y[0], a1, b1, c1);
1835  if (a1*x[1]+b1*y[1]+c1 <= 0.1*del) return true;
1836 
1837  // 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
1838 
1839  int inode, iedge;
1840  double a2, b2, c2, a3, b3, c3;
1841 
1842  findABC(x[0], y[0], x[1], y[1], a2, b2, c2);
1843  findABC(x[1], y[1], x[2], y[2], a3, b3, c3);
1844  iedge = iedge2;
1845  for (;;) {
1846  iedge = edges[iedge].inext;
1847  if (edges[iedge].inext == iedge1) return false;
1848  inode = edges[iedge].i2;
1849  if (inode == k[0]) continue;
1850  if (inode == k[1]) continue;
1851  if (inode == k[2]) continue;
1852  x[1] = nodes[inode].v[ix];
1853  y[1] = nodes[inode].v[iy];
1854  if (a1*x[1]+b1*y[1]+c1 < -0.1*del) continue;
1855  if (a2*x[1]+b2*y[1]+c2 < -0.1*del) continue;
1856  if (a3*x[1]+b3*y[1]+c3 < -0.1*del) continue;
1857  return true;
1858  }
1859 }
1860 
1861 inline
1862 bool BooleanProcessor::triangulateContour(int ix, int iy, int ihead)
1863 /***********************************************************************
1864  * *
1865  * Name: BooleanProcessor::triangulateContour Date: 06.03.00 *
1866  * Author: E.Chernyaev Revised: *
1867  * *
1868  * Function: Triangulate external contour *
1869  * *
1870  ***********************************************************************/
1871 {
1872 
1873  //int draw_flag = 0;
1874  //if (draw_flag) draw_contour(5, 3, ihead);
1875 
1876  // C L O S E C O N T O U R
1877 
1878  int ipnext = ihead, nnode = 1;
1879  for (;;) {
1880  if (edges[ipnext].inext > 0) {
1881  ipnext = edges[ipnext].inext;
1882  nnode++;
1883  }else{
1884  edges[ipnext].inext = ihead;
1885  break;
1886  }
1887  }
1888 
1889  // L O O P A L O N G C O N T O U R
1890 
1891  //std::cerr << "debug : contour : begin : =================" << std::endl;
1892  //dump();//debug
1893 
1894  int iedge1, iedge2, iedge3, istart = 0;
1895  for (;;) {
1896  iedge1 = edges[ipnext].inext;
1897  iedge2 = edges[iedge1].inext;
1898 
1899  //std::cerr << "debug :"
1900  // << " ipnext " << ipnext
1901  // << " iedge1 " << iedge1
1902  // << " iedge2 " << iedge2
1903  // << " : istart " << istart
1904  // << " , nnode " << nnode
1905  // << std::endl;
1906 
1907  if (istart == 0) {
1908  istart = iedge1;
1909  if (nnode <= 3) {
1910  iedge3 = edges[iedge2].inext;
1911  edges[iedge1].iface1 = int(faces.size());
1912  edges[iedge2].iface1 = int(faces.size());
1913  edges[iedge3].iface1 = int(faces.size());
1914  edges[iedge3].inext = 0;
1915  faces.push_back(ExtFace(edges,0)); //G.Barrand : ok ?
1916  faces.back().iold = iedge1;
1917  faces.back().inew = ORIGINAL_FACE;
1918 
1919  //if (draw_flag) draw_contour(4, 2, iedge1);
1920 
1921  break;
1922  }
1923  }else if (istart == iedge1) {
1924  //G.Barrand processor_error = 1; //use returned status.
1925 #ifdef TOOLS_HEP_BP_OUT_ERR
1926  std::cerr
1927  << "BooleanProcessor::triangulateContour : "
1928  << "could not generate a triangle (infinite loop)"
1929  << std::endl;
1930 #endif
1931  return false; //G.Barrand
1932  }
1933 
1934  // 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
1935 
1936  if(checkTriangle(iedge1,iedge2,ix,iy)) {
1937  ipnext = edges[ipnext].inext;
1938  continue;
1939  }
1940 
1941  // M O D I F Y C O N T O U R
1942 
1943  int i1 = edges[iedge1].i1;
1944  int i3 = edges[iedge2].i2;
1945  int iface1 = edges[iedge1].iface1;
1946  int iface2 = int(faces.size());
1947 
1948  edges[ipnext].inext = int(edges.size());
1949  edges.push_back(ExtEdge(i1, i3, iface1, -1*(int(edges.size())+1), -1));
1950  edges.back().inext = edges[iedge2].inext;
1951 
1952  // A D D N E W T R I A N G L E T O T H E L I S T
1953 
1954  edges[iedge2].inext = int(edges.size());
1955  edges.push_back(ExtEdge(i3, i1, iface2, -1*(int(edges.size())-1), -1));
1956  faces.push_back(ExtFace(edges,0)); //G.Barrand : ok ?
1957  faces.back().iold = iedge1;
1958  faces.back().inew = ORIGINAL_FACE;
1959  edges[iedge1].iface1 = iface2;
1960  edges[iedge2].iface1 = iface2;
1961  ipnext = edges[ipnext].inext;
1962  istart = 0;
1963  nnode--;
1964 
1965  //if (draw_flag) draw_contour(4, 2, iedge1);
1966 
1967  }
1968 
1969  return true;
1970 }
1971 
1972 inline
1973 bool BooleanProcessor::modifyReference(int iface, int i1, int i2, int iref)
1974 /***********************************************************************
1975  * *
1976  * Name: BooleanProcessor::modifyReference Date: 13.03.00 *
1977  * Author: E.Chernyaev Revised: *
1978  * *
1979  * Function: Modify reference to the neighbouring face *
1980  * *
1981  ***********************************************************************/
1982 {
1983  int iedge = faces[iface].iold;
1984  while (iedge > 0) {
1985  if (edges[iedge].i1 == i2 && edges[iedge].i2 == i1) {
1986  edges[iedge].iface2 = iref;
1987  return true;
1988  }
1989  iedge = edges[iedge].inext;
1990  }
1991  //G.Barrand : processor_error = 1;
1992 #ifdef TOOLS_HEP_BP_OUT_ERR
1993  std::cerr
1994  << "BooleanProcessor::modifyReference : could not find the edge, "
1995  << "iface=" << iface << ", i1,i2=" << i1 << "," << i2 << ", iref=" << iref
1996  << std::endl;
1997 #endif
1998  return false; //G.Barrand
1999 }
2000 
2001 inline
2002 bool BooleanProcessor::triangulateFace(int iface)
2003 /***********************************************************************
2004  * *
2005  * Name: BooleanProcessor::triangulateFace Date: 02.03.00 *
2006  * Author: E.Chernyaev Revised: *
2007  * *
2008  * Function: Triangulation of an extended face *
2009  * *
2010  ***********************************************************************/
2011 {
2012 
2013  // 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
2014  // S E T IX, IY, IZ
2015 
2016  //HVNormal3D normal = faces[iface].plane.normal();
2017  const HVNormal3D& normal = faces[iface].plane.normal();
2018 
2019  int ix, iy, iz = 0;
2020  //G.Barrand : if (std::abs(normal[1]) > std::abs(normal[iz])) iz = 1;
2021  //G.Barrand : if (std::abs(normal[2]) > std::abs(normal[iz])) iz = 2;
2022  if (::fabs(normal[1]) > ::fabs(normal[iz])) iz = 1; //G.Barrand
2023  if (::fabs(normal[2]) > ::fabs(normal[iz])) iz = 2; //G.Barrand
2024  if (normal[iz] > 0) {
2025  ix = (iz+1)%3; iy = (ix+1)%3;
2026  }else{
2027  iy = (iz+1)%3; ix = (iy+1)%3;
2028  }
2029 
2030  // F I L L L I S T S O F C O N T O U R S
2031 
2032  external_contours.clear();
2033  internal_contours.clear();
2034  double z;
2035  int i1, i2, ifirst, iedge, icontour = faces[iface].iold;
2036  while (icontour > 0) {
2037  iedge = icontour;
2038  ifirst = edges[iedge].i1;
2039  z = 0.0;
2040  for(;;) {
2041  if (iedge > 0) {
2042  i1 = edges[iedge].i1;
2043  i2 = edges[iedge].i2;
2044  ExtNode& node_1 = nodes[i1];
2045  ExtNode& node_2 = nodes[i2];
2046  z += node_1.v[ix]*node_2.v[iy]-node_2.v[ix]*node_1.v[iy];
2047  if (ifirst != i2) {
2048  iedge = edges[iedge].inext;
2049  continue;
2050  }else{
2051  if (z > del*del) {
2052  external_contours.push_back(icontour);
2053  }else if (z < -del*del) {
2054  internal_contours.push_back(icontour);
2055  }else{
2056  //G.Barrand : processor_error = 1; //use returned value.
2057 #ifdef TOOLS_HEP_BP_OUT_ERR
2058  std::cerr
2059  << "BooleanProcessor::triangulateFace : too small contour"
2060  << " z " << z
2061  << " del*del " << del*del
2062  << std::endl;
2063 #endif
2064  return false; //G.Barrand.
2065  }
2066  icontour = edges[iedge].inext;
2067  edges[iedge].inext = 0;
2068  break;
2069  }
2070  }else{
2071  //G.Barrand : processor_error = 1; //use returned value.
2072 #ifdef TOOLS_HEP_BP_OUT_ERR
2073  std::cerr
2074  << "BooleanProcessor::triangulateFace : broken contour"
2075  << std::endl;
2076 #endif
2077  icontour = 0;
2078  //G.Barrand : break;
2079  return false;
2080  }
2081  }
2082  }
2083 
2084  // G E T R I D O F I N T E R N A L C O N T O U R S
2085 
2086  int kint, kext;
2087  for (kint=0; kint < (int)internal_contours.size(); kint++) {
2088  for (kext=0; kext < (int)external_contours.size(); kext++) {
2089  mergeContours(ix, iy, kext, kint);
2090  if (internal_contours[kint] == 0) break;
2091  }
2092  if (kext == (int)external_contours.size()) {
2093  //G.Barrand : processor_error = 1; //use returned value.
2094 #ifdef TOOLS_HEP_BP_OUT_ERR
2095  std::cerr
2096  << "BooleanProcessor::triangulateFace : "
2097  << "could not merge internal contour " << kint
2098  << std::endl;
2099 #endif
2100  return false; //G.Barrand
2101  }
2102  }
2103 
2104  // T R I A N G U L A T E C O N T O U R S
2105 
2106  int nface = int(faces.size());
2107  for (kext=0; kext < (int)external_contours.size(); kext++) {
2108  if(!triangulateContour(ix, iy, external_contours[kext])) {
2109  //G.Barrand : processor_error = 1; //use returned value.
2110 #ifdef TOOLS_HEP_BP_OUT_ERR
2111  std::cerr
2112  << "BooleanProcessor::triangulateFace : "
2113  << "triangulateContour failed."
2114  << std::endl;
2115 #endif
2116  //break; //G.Barrand : ok ?
2117  //return; //G.Barrand : ok ?
2118  return false; //G.Barrand
2119  }
2120  }
2121  faces[iface].inew = UNSUITABLE_FACE;
2122 
2123  // M O D I F Y R E F E R E N C E S
2124 
2125  for (int ifa=nface; ifa<(int)faces.size(); ifa++) {
2126  iedge = faces[ifa].iold;
2127  while (iedge > 0) {
2128  if (edges[iedge].iface1 != ifa) {
2129  //G.Barrand : processor_error = 1; //use returned value.
2130 #ifdef TOOLS_HEP_BP_OUT_ERR
2131  std::cerr
2132  << "BooleanProcessor::triangulateFace : wrong reference to itself, "
2133  << "iface=" << ifa << ", iface1=" << edges[iedge].iface1
2134  << std::endl;
2135 #endif
2136  return false; //G.Barrand
2137  }else if (edges[iedge].iface2 > 0) {
2138  if(!modifyReference(edges[iedge].iface2,
2139  edges[iedge].i1, edges[iedge].i2, ifa)) { //G.Barrand
2140  return false; //G.Barrand
2141  }
2142  }else if (edges[iedge].iface2 < 0) {
2143  edges[iedge].iface2 = edges[-edges[iedge].iface2].iface1;
2144  }
2145  iedge = edges[iedge].inext;
2146  }
2147  }
2148 
2149  return true; //G.Barrand.
2150 }
2151 
2152 inline
2153 polyhedron BooleanProcessor::createPolyhedron()
2154 /***********************************************************************
2155  * *
2156  * Name: BooleanProcessor::createPolyhedron() Date: 14.03.00 *
2157  * Author: E.Chernyaev Revised: *
2158  * *
2159  * Function: Create polyhedron. *
2160  * *
2161  ***********************************************************************/
2162 {
2163  int i, iedge, nnode = 0, nface = 0;
2164 
2165  // R E N U M E R A T E N O D E S A N D F A C E S
2166 
2167  for (i=1; i<(int)nodes.size(); i++) nodes[i].s = 0;
2168 
2169  for (i=1; i<(int)faces.size(); i++) {
2170  if (faces[i].inew == ORIGINAL_FACE) {
2171  faces[i].inew = ++nface;
2172  iedge = faces[i].iold;
2173  while (iedge > 0) {
2174  nodes[edges[iedge].i1].s = 1;
2175  iedge = edges[iedge].inext;
2176  }
2177  }else{
2178  faces[i].inew = 0;
2179  }
2180  }
2181 
2182  for (i=1; i<(int)nodes.size(); i++) {
2183  if (nodes[i].s == 1) nodes[i].s = ++nnode;
2184  }
2185 
2186  // A L L O C A T E M E M O R Y
2187 
2188  ExtPolyhedron polyhedron;
2189  if (nface == 0) return polyhedron;
2190  polyhedron.AllocateMemory(nnode, nface);
2191 
2192  // S E T N O D E S
2193 
2194  for (i=1; i<(int)nodes.size(); i++) {
2195  if (nodes[i].s != 0) polyhedron.pV[nodes[i].s] = nodes[i].v;
2196  }
2197 
2198  // S E T F A C E S
2199 
2200  int k, v[4], f[4];
2201  for (i=1; i<(int)faces.size(); i++) {
2202  if (faces[i].inew == 0) continue;
2203  v[3] = f[3] = k = 0;
2204  iedge = faces[i].iold;
2205  while (iedge > 0) {
2206  if (k > 3) {
2207 #ifdef TOOLS_HEP_BP_VERBOSE
2208  std::cerr
2209  << "BooleanProcessor::createPolyhedron : too many edges"
2210  << std::endl;
2211 #endif
2212  break;
2213  }
2214  v[k] = nodes[edges[iedge].i1].s;
2215  if (edges[iedge].ivis < 0) v[k] = -v[k];
2216  f[k] = faces[edges[iedge].iface2].inew;
2217  iedge = edges[iedge].inext;
2218  k++;
2219  }
2220  if (k < 3) {
2221 #ifdef TOOLS_HEP_BP_VERBOSE
2222  std::cerr
2223  << "BooleanProcessor::createPolyhedron : "
2224  << "face has only " << k << " edges"
2225  << std::endl;
2226 #endif
2227  }
2228  //polyhedron.pF[faces[i].inew] =
2229  // SbFacet(v[0],f[0], v[1],f[1], v[2],f[2], v[3],f[3]);
2230  polyhedron.pF[faces[i].inew].set( //G.Barrand
2231  v[0],f[0], v[1],f[1], v[2],f[2], v[3],f[3]);
2232  }
2233  return polyhedron;
2234 }
2235 
2236 //G.Barrand : begin
2237 //int BooleanProcessor::ishift = 0; //does not work with pure header code.
2238 inline int& BooleanProcessor::get_shift() { //& not thread safe.
2239  static int s_ishift;
2240  return s_ishift;
2241 }
2242 
2243 inline void BooleanProcessor::set_shift(int a_shift) {
2244  get_shift() = a_shift;
2245 }
2246 inline void BooleanProcessor::inc_shift() {
2247  int i = get_shift();
2248  i++;
2249  if(i>=get_num_shift()) i = 0;
2250  get_shift() = i;
2251 }
2252 
2253 inline int BooleanProcessor::get_num_shift() { return 8;}
2254 //G.Barrand : end
2255 
2256 //G.Barrand :
2257 inline
2258 bool BooleanProcessor::inc_try_count(
2259 #ifdef TOOLS_HEP_BP_VERBOSE
2260  int a_op,const polyhedron& a,const polyhedron& b
2261 #else
2262  int,const polyhedron&,const polyhedron&
2263 #endif
2264 ,int& a_try_count){
2265  if(a_try_count>get_num_shift()) {
2266 #ifdef TOOLS_HEP_BP_VERBOSE
2267  std::cerr << "BooleanProcessor: "
2268  << " all shifts tried. Boolean operation (" << a_op << ") failure."
2269  << " a name \"" << a.getName() << "\""
2270  << " b name \"" << b.getName() << "\""
2271  << std::endl;
2272 #endif
2273  return false; //stop
2274  }
2275 
2276 #ifdef TOOLS_HEP_BP_VERBOSE
2277  std::cerr
2278  << "BooleanProcessor::execute : try another tilt..."
2279  << std::endl;
2280 #endif
2281 
2282  a_try_count++;
2283  return true;
2284 }
2285 
2286 inline
2287 polyhedron BooleanProcessor::execute(int op,
2288  const polyhedron & a,
2289  const polyhedron & b,
2290  int& a_err) //G.Barrand
2291 /***********************************************************************
2292  * *
2293  * Name: BooleanProcessor::execute Date: 10.12.99 *
2294  * Author: E.Chernyaev Revised: *
2295  * *
2296  * Function: Execute boolean operation. *
2297  * *
2298  ***********************************************************************/
2299 {
2300  //static int ishift = 0; //G.Barrand
2301  //static double shift[8][3] = {
2302  static double shift[/*NUM_SHIFT*/ 8][3] = { //G.Barrand
2303  { 31, 23, 17},
2304  { -31, -23, -17},
2305  { -23, 17, 31},
2306  { 23, -17, -31},
2307  { -17, -31, 23},
2308  { 17, 31, -23},
2309  { 31, -23, 17},
2310  { -31, 23, -17}
2311  };
2312 
2313 
2314  //std::cerr << "BooleanProcessor::execute : ++++++++++++++++++++++"
2315  // << a.getName()
2316  // << b.getName()
2317  // << std::endl;
2318 
2319  // I N I T I A T E P R O C E S S O R
2320 
2321  operation = op;
2322  nodes.clear(); nodes.push_back(CRAZY_POINT());
2323  edges.clear(); edges.push_back(ExtEdge());
2324  faces.clear(); faces.push_back(ExtFace(edges,0)); //G.Barrand : ok ?
2325 
2326  // T A K E P O L Y H E D R A
2327 
2328  ifaces1 = int(faces.size());
2329  if(!takePolyhedron(a,0,0,0)){ //G.Barrand : use returned value.
2330  // corrapted polyhedron
2331 #ifdef TOOLS_HEP_BP_OUT_ERR
2332  std::cerr
2333  << "BooleanProcessor: a corrapted input polyhedron"
2334  << std::endl;
2335 #endif
2336  a_err = 1; //G.Barrand
2337  return polyhedron();
2338  }
2339 
2340  ifaces2 = int(faces.size());
2341  if(!takePolyhedron(b,0,0,0)){
2342  // corrapted polyhedron
2343 #ifdef TOOLS_HEP_BP_OUT_ERR
2344  std::cerr
2345  << "BooleanProcessor: b corrapted input polyhedron"
2346  << std::endl;
2347 #endif
2348  a_err = 1; //G.Barrand
2349  return polyhedron();
2350  }
2351 
2352  if (ifaces1 == ifaces2) { // a is empty
2353  a_err = 0; //G.Barrand
2354  switch (operation) {
2355  case OP_UNION:
2356  return b;
2357  case OP_INTERSECTION:
2358 #ifdef TOOLS_HEP_BP_VERBOSE
2359  std::cerr
2360  << "BooleanProcessor: intersection with empty polyhedron"
2361  << std::endl;
2362 #endif
2363  return polyhedron();
2364  case OP_SUBTRACTION:
2365 #ifdef TOOLS_HEP_BP_VERBOSE
2366  std::cerr
2367  << "BooleanProcessor: subtraction from empty polyhedron"
2368  << std::endl;
2369 #endif
2370  return polyhedron();
2371  }
2372  }
2373  if (ifaces2 == (int)faces.size()) { // b is empty
2374  a_err = 0; //G.Barrand
2375  switch (operation) {
2376  case OP_UNION:
2377  return a;
2378  case OP_INTERSECTION:
2379 #ifdef TOOLS_HEP_BP_VERBOSE
2380  std::cerr
2381  << "BooleanProcessor: intersection with empty polyhedron"
2382  << std::endl;
2383 #endif
2384  return polyhedron();
2385  case OP_SUBTRACTION:
2386  return a;
2387  }
2388  }
2389 
2390  // 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
2391 
2392  del = findMinMax();
2393 
2394  // W O R K A R O U N D T O A V O I D I E A N D E E
2395 
2396  int try_count = 1;
2397  while(true) { //G.Barrand
2398 
2399  double ddxx = del*shift[get_shift()][0];
2400  double ddyy = del*shift[get_shift()][1];
2401  double ddzz = del*shift[get_shift()][2];
2402  //ishift++; if (ishift == get_num_shift()) ishift = 0;
2403  inc_shift();
2404 
2405  operation = op;
2406  nodes.clear(); nodes.push_back(CRAZY_POINT());
2407  edges.clear(); edges.push_back(ExtEdge());
2408  faces.clear(); faces.push_back(ExtFace(edges,0)); //G.Barrand : ok ?
2409 
2410  ifaces1 = int(faces.size());
2411  if(!takePolyhedron(a,0,0,0)) {
2412 #ifdef TOOLS_HEP_BP_OUT_ERR
2413  std::cerr
2414  << "BooleanProcessor: shifted a corrapted input polyhedron"
2415  << std::endl;
2416 #endif
2417  if(!inc_try_count(op,a,b,try_count)) {a_err = 1;return a;}
2418  continue; //try another shift.
2419  }
2420 
2421  ifaces2 = int(faces.size());
2422  if(!takePolyhedron(b,ddxx,ddyy,ddzz)){
2423 #ifdef TOOLS_HEP_BP_OUT_ERR
2424  std::cerr
2425  << "BooleanProcessor: shifted b corrapted input polyhedron"
2426  << std::endl;
2427 #endif
2428  if(!inc_try_count(op,a,b,try_count)) {a_err = 1;return a;}
2429  continue; //try another shift.
2430  }
2431 
2432  del = findMinMax();
2433 
2434  // P R E S E L E C T O U T S I D E F A C E S
2435 
2436  iout1 = iout2 = 0;
2437  selectOutsideFaces(ifaces1, iout1);
2438  selectOutsideFaces(ifaces2, iout2);
2439 
2440  // 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
2441 
2442  int ifa1, ifa2;
2443  iunk1 = iunk2 = 0;
2444  if (iout1 != 0 || iout2 != 0) {
2445  for(;;) {
2446  ifa1 = iunk1;
2447  ifa2 = iunk2;
2448  selectOutsideFaces(ifaces1, iunk1);
2449  selectOutsideFaces(ifaces2, iunk2);
2450  if (iunk1 == ifa1 && iunk2 == ifa2) break;
2451  findMinMax();
2452  }
2453  }
2454 
2455  // F I N D N E W E D G E S
2456 
2457  {int processor_error = 0; //G.Barrand
2458  if (ifaces1 != 0 && ifaces2 != 0 ) {
2459  ifa1 = ifaces1;
2460  while (ifa1 > 0) {
2461  ifa2 = ifaces2;
2462  while (ifa2 > 0) {
2463  if(!testFaceVsFace(ifa1, ifa2)) {
2464  processor_error = 1; //G.Barrand
2465  break; //G.Barrand
2466  }
2467  ifa2 = faces[ifa2].inext;
2468  }
2469  if(processor_error) break; //G.Barrand
2470  ifa1 = faces[ifa1].inext;
2471  }
2472  }
2473  if (processor_error) {
2474  PROCESSOR_ERROR(4);
2475  if(!inc_try_count(op,a,b,try_count)) {a_err = 1;return a;}
2476  continue; //try another shift.
2477  }} //G.Barrand
2478 
2479  // C O N S T R U C T N E W F A C E S
2480 
2481  if(!assembleNewFaces((operation == OP_INTERSECTION) ? 1 : 0, ifaces1)){
2482  //G.Barrand
2483  PROCESSOR_ERROR(5);
2484  if(!inc_try_count(op,a,b,try_count)) {a_err = 1;return a;}
2485  continue; //try another shift.
2486  }
2487  if(!assembleNewFaces((operation == OP_UNION) ? 0 : 1, ifaces2)){
2488  //G.Barrand
2489  PROCESSOR_ERROR(6);
2490  if(!inc_try_count(op,a,b,try_count)) {a_err = 1;return a;}
2491  continue; //try another shift.
2492  }
2493  // A S S E M B L E S U I T A B L E F A C E S
2494 
2495  if(!initiateLists()) { //G.Barrand : return status.
2496 #ifdef TOOLS_HEP_BP_OUT_ERR
2497  std::cerr
2498  << "BooleanProcessor: initiateLists failed."
2499  << std::endl;
2500 #endif
2501  if(!inc_try_count(op,a,b,try_count)) {a_err = 1;return a;}
2502  continue; //try another shift.
2503  }
2504 
2505  {int processor_error = 0;
2506  for (;;) {
2507  if(!assemblePolyhedra()) { //G.Barrand : handle returned value.
2508  processor_error = 1;
2509 #ifdef TOOLS_HEP_BP_OUT_ERR
2510  std::cerr
2511  << "BooleanProcessor: assemblePolyhedra failed."
2512  << std::endl;
2513 #endif
2514  break;
2515  }
2516  if (unknown_faces.front() != 0) {
2517  processor_error = 1;
2518 #ifdef TOOLS_HEP_BP_OUT_ERR
2519  std::cerr
2520  << "BooleanProcessor::execute : unknown faces !!!"
2521  << std::endl;
2522 #endif
2523  break;
2524  }
2525  break;
2526  }
2527  if (processor_error) {
2528  PROCESSOR_ERROR(7);
2529  if(!inc_try_count(op,a,b,try_count)) {a_err = 1;return a;}
2530  continue; //try another shift.
2531  }}
2532 
2533  // T R I A N G U L A T E A C C E P T E D F A C E S
2534 
2535  {int processor_error = 0;
2536  ifa1 = result_faces.front();
2537  while (ifa1 > 0) {
2538  ifa2 = ifa1;
2539  ifa1 = faces[ifa2].inext;
2540  if(faces[ifa2].inew == NEW_FACE) {
2541  if(!triangulateFace(ifa2)) { //G.Barrand : use returned status.
2542  processor_error = 1;
2543  break; //G.Barrand
2544  }
2545  }
2546  }
2547  if (processor_error) {
2548  PROCESSOR_ERROR(8);
2549  if(!inc_try_count(op,a,b,try_count)) {a_err = 1;return a;}
2550  continue; //try another shift.
2551  }}
2552 
2553 #ifdef TOOLS_HEP_BP_VERBOSE
2554  if(try_count!=1) {
2555  std::cerr
2556  << "BooleanProcessor::execute : had converged."
2557  << std::endl;
2558  }
2559 #endif
2560  break;
2561 
2562  } //G.Barrand : end while shift.
2563 
2564  // C R E A T E P O L Y H E D R O N
2565 
2566  a_err = 0;
2567  return createPolyhedron();
2568 }
2569 
2570 
2571 //#include <cfortran.h>
2572 //#include <higz.h>
2573 //#include "zbuf.h"
2574 //void BooleanProcessor::draw()
2575 /***********************************************************************
2576  * *
2577  * Name: BooleanProcessor::draw Date: 10.12.99 *
2578  * Author: E.Chernyaev Revised: *
2579  * *
2580  * Function: Draw *
2581  * *
2582  ***********************************************************************/
2583 /*
2584 {
2585  int II;
2586  int icol, i1, i2, iedge, iface, ilist[4];
2587  float p1[3], p2[3];
2588 
2589  ilist[0] = ifaces1;
2590  ilist[1] = ifaces2;
2591  ilist[2] = iout1;
2592  ilist[3] = iout2;
2593 
2594  for (int i=0; i<4; i++) {
2595 
2596  if (i == 0) cout << "========= Ifaces_1" << endl;
2597  if (i == 1) cout << "========= Ifaces_2" << endl;
2598  if (i == 2) cout << "========= Iout_1" << endl;
2599  if (i == 3) cout << "========= Iout_2" << endl;
2600 
2601  icol = i+1;
2602  iface = ilist[i];
2603  while (iface > 0) {
2604 
2605  cout << "iface = " << iface << endl;
2606  cout << "--- iold" << endl;
2607 
2608  iedge = faces[iface].iold;
2609  icol = 2;
2610 
2611  while (iedge > 0) {
2612 
2613  cout << " iegde = " << iedge
2614  << " i1,i2 =" << edges[iedge].i1 << "," << edges[iedge].i2
2615  << " iface1,iface2 = "
2616  << edges[iedge].iface1 << "," << edges[iedge].iface2
2617  << endl;
2618 
2619  i1 = edges[iedge].i1;
2620  p1[0] = nodes[i1].v.x();
2621  p1[1] = nodes[i1].v.y();
2622  p1[2] = nodes[i1].v.z();
2623  IHWTON(p1,p1);
2624  i2 = edges[iedge].i2;
2625  p2[0] = nodes[i2].v.x();
2626  p2[1] = nodes[i2].v.y();
2627  p2[2] = nodes[i2].v.z();
2628  IHWTON(p2,p2);
2629 // icol = (edges[iedge].ivis > 0) ? 1 : 2;
2630  IHZLIN(icol,p1[0],p1[1],p1[2], p2[0],p2[1],p2[2]);
2631  iedge = edges[iedge].inext;
2632  }
2633 
2634  cout << "--- inew" << endl;
2635 
2636  iedge = faces[iface].inew;
2637  icol = 3;
2638 
2639  while (iedge > 0) {
2640 
2641  cout << " iegde = " << iedge
2642  << " i1,i2 =" << edges[iedge].i1 << "," << edges[iedge].i2
2643  << " iface1,iface2 = "
2644  << edges[iedge].iface1 << "," << edges[iedge].iface2
2645  << endl;
2646 
2647  i1 = edges[iedge].i1;
2648  p1[0] = nodes[i1].v.x();
2649  p1[1] = nodes[i1].v.y();
2650  p1[2] = nodes[i1].v.z();
2651  IHWTON(p1,p1);
2652  i2 = edges[iedge].i2;
2653  p2[0] = nodes[i2].v.x();
2654  p2[1] = nodes[i2].v.y();
2655  p2[2] = nodes[i2].v.z();
2656  IHWTON(p2,p2);
2657 // icol = (edges[iedge].ivis > 0) ? 1 : 2;
2658  IHZLIN(icol,p1[0],p1[1],p1[2], p2[0],p2[1],p2[2]);
2659  iedge = edges[iedge].inext;
2660  }
2661  iface = faces[iface].inext;
2662 
2663  IHZTOX(0,100,100);
2664  ixupdwi(0);
2665  cin >> II;
2666  ixclrwi();
2667  IHZCLE(0);
2668  }
2669  }
2670 }
2671 */
2672 
2673 /*
2674 //--------------------------------------------------------------------
2675 inline void
2676 BooleanProcessor::draw_edge(int icol, int iedge) {
2677  int i1, i2;
2678  float p1[3], p2[3];
2679 
2680  i1 = edges[iedge].i1;
2681  p1[0] = nodes[i1].v.x();
2682  p1[1] = nodes[i1].v.y();
2683  p1[2] = nodes[i1].v.z();
2684  IHWTON(p1,p1);
2685  i2 = edges[iedge].i2;
2686  p2[0] = nodes[i2].v.x();
2687  p2[1] = nodes[i2].v.y();
2688  p2[2] = nodes[i2].v.z();
2689  IHWTON(p2,p2);
2690  IHZLIN(icol,p1[0],p1[1],p1[2], p2[0],p2[1],p2[2]);
2691 }
2692 
2693 //--------------------------------------------------------------------
2694 inline void
2695 BooleanProcessor::draw_contour(int i1col, int i2col, int ihead) {
2696  int iedge, icol;
2697  iedge = ihead;
2698  while (iedge > 0) {
2699  icol = (edges[iedge].ivis > 0) ? i1col : i2col;
2700  draw_edge(icol, iedge);
2701  iedge = edges[iedge].inext;
2702  }
2703 
2704  IHZTOX(0,100,100);
2705  ixupdwi(0);
2706 
2707  int i;
2708  std::cin >> i;
2709 }
2710 
2711 //--------------------------------------------------------------------
2712 inline void
2713 BooleanProcessor::print_face(int iface) {
2714  cout.precision(3);
2715  cout << "\n====== Face N " << iface << endl;
2716  cout << "iedges[4] = "
2717  << faces[iface].iedges[0] << ", "
2718  << faces[iface].iedges[1] << ", "
2719  << faces[iface].iedges[2] << ", "
2720  << faces[iface].iedges[3] << endl;
2721  cout << "rmin[3] = "
2722  << faces[iface].rmin[0] << ", "
2723  << faces[iface].rmin[1] << ", "
2724  << faces[iface].rmin[2] << endl;
2725  cout << "rmax[3] = "
2726  << faces[iface].rmax[0] << ", "
2727  << faces[iface].rmax[1] << ", "
2728  << faces[iface].rmax[2] << endl;
2729  cout << "iprev,inext = "
2730  << faces[iface].iprev << ", "
2731  << faces[iface].inext << endl;
2732  cout << "iold = " << faces[iface].iold << endl;
2733  for(int i = faces[iface].iold; i != 0;) {
2734  print_edge(i);
2735  i = edges[abs(i)].inext;
2736  }
2737 
2738  cout << "inew = ";
2739  switch (faces[iface].inew) {
2740  case UNKNOWN_FACE:
2741  cout << "UNKNOWN_FACE" << endl;
2742  break;
2743  case ORIGINAL_FACE:
2744  cout << "ORIGINAL_FACE" << endl;
2745  break;
2746  case NEW_FACE:
2747  cout << "NEW_FACE" << endl;
2748  break;
2749  case UNSUITABLE_FACE:
2750  cout << "UNSUITABLE_FACE" << endl;
2751  break;
2752  case DEFECTIVE_FACE:
2753  cout << "DEFECTIVE_FACE" << endl;
2754  break;
2755  default:
2756  cout << faces[iface].inew << endl;
2757  for(int k = faces[iface].inew; k != 0;) {
2758  print_edge(k);
2759  k = edges[abs(k)].inext;
2760  }
2761  }
2762 }
2763 
2764 //--------------------------------------------------------------------
2765 inline void
2766 BooleanProcessor::print_edge(int iedge) {
2767  cout << "==== Edge N " << iedge << endl;
2768  int i = std::abs(iedge);
2769  int i1 = edges[i].i1;
2770  int i2 = edges[i].i2;
2771  cout << "node[" << i1 << "] = "
2772  << nodes[i1].v.x() << ", "
2773  << nodes[i1].v.y() << ", "
2774  << nodes[i1].v.z() << endl;
2775 
2776  cout << "node[" << i2 << "] = "
2777  << nodes[i2].v.x() << ", "
2778  << nodes[i2].v.y() << ", "
2779  << nodes[i2].v.z() << endl;
2780 
2781  cout << "iface1,iface2,ivis,inext = "
2782  << edges[i].iface1 << ", "
2783  << edges[i].iface2 << ", "
2784  << edges[i].ivis << ", "
2785  << edges[i].inext << endl;
2786 }
2787 */
2788 
2789 inline void BooleanProcessor::dump(std::ostream& a_out) {//G.Barrand
2790  size_t number = nodes.size();
2791  a_out << "nodes : " << number << std::endl;
2792  for(size_t index=0;index<number;index++) {
2793  const ExtNode& node = nodes[index];
2794  a_out << " " << index
2795  << " x = " << node.v[0]
2796  << " y = " << node.v[1]
2797  << " z = " << node.v[2]
2798  << std::endl;
2799  }
2800 }
2801 
2802 }}
Generated by doxygen 1.8.20