source: trunk/source/graphics_reps/src/BooleanProcessor.src@ 1228

Last change on this file since 1228 was 1140, checked in by garnier, 16 years ago

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