45#include "G4Plane3D.hh"
46#include "G4Point3D.hh"
47#include "G4Normal3D.hh"
54#define ExtNode HEPVis_ExtNode
55#define ExtEdge HEPVis_ExtEdge
56#define ExtFace HEPVis_ExtFace
57#define FaceList HEPVis_FaceList
58#define ExtPolyhedron HEPVis_ExtPolyhedron
59#define BooleanProcessor HEPVis_BooleanProcessor
61#define HepPolyhedron SbPolyhedron
62#define G4Facet SbFacet
82#define INITIAL_SIZE 200
83#define CRAZY_POINT HVPoint3D(-10.e+6, -10.e+6, -10.e+6)
85#define GRANULARITY 10.e+5
87#define SWAP(A,B) w = A; A = B; B = w
90#define OP_INTERSECTION 1
91#define OP_SUBTRACTION 2
97#define NON_PLANAR_FACE 4
100#define ORIGINAL_FACE -1
102#define UNSUITABLE_FACE -3
103#define DEFECTIVE_FACE -4
154 :
v(vertex),
s(status) {}
160 if (&
node ==
this)
return *
this;
177 ExtEdge(
int k1=0,
int k2=0,
int kface1=0,
int kface2=0,
int kvis=0) :
187 if (&edge ==
this)
return *
this;
218 ExtFace(std::vector<ExtEdge>& a_edges,
int iedge)
221 {
for (
int i=0; i<4; i++) {
iedges[i] = 0; }}
222 {
for (
int i=0; i<3; i++) {
rmin[i] = 0;
rmax[i] = 0; }}
237 if (&face ==
this)
return *
this;
307 return HepPolyhedron::operator = (from);
343 void divideEdge(
int & i1,
int & i2);
344 void insertEdge(
const ExtEdge & edge);
348 void testFaceVsFace(
int iface1,
int iface2);
349 void invertNewEdges(
int iface);
350 void checkDoubleEdges(
int iface);
351 void assembleFace(
int what,
int iface);
352 void assembleNewFaces(
int what,
int ihead);
353 void initiateLists();
354 void assemblePolyhedra();
355 void findABC(
double x1,
double y1,
double x2,
double y2,
356 double &
a,
double &b,
double &c)
const;
357 int checkDirection(
double *
x,
double *
y)
const;
358 int checkIntersection(
int ix,
int iy,
int i1,
int i2)
const;
359 void mergeContours(
int ix,
int iy,
int kext,
int kint);
360 int checkTriangle(
int iedge1,
int iedge2,
int ix,
int iy)
const;
361 void triangulateContour(
int ix,
int iy,
int ihead);
362 void modifyReference(
int iface,
int i1,
int i2,
int iref);
363 void triangulateFace(
int iface);
383 for (
int i = 0; i < 3; i++) {
405 static int get_shift();
406 static void set_shift(
int);
407 static int get_num_shift();
420 int iEprev, iEcur, iEnext;
422 iEprev = 0; iEcur =
iold;
433 iEprev = 0; iEcur =
inew;
452 double dx,
double dy,
double dz)
462 int i, k, nnode, iNodes[5], iVis[4], iFaces[4];
463 int dnode =
m_nodes.size() - 1;
464 int dface =
m_faces.size() - 1;
473 for (i=1; i <= p.GetNoVertices(); i++) {
475 ppp = p.GetVertex(i);
476 ppp.setX(ppp.x()+dx);
477 ppp.setY(ppp.y()+dy);
478 ppp.setZ(ppp.z()+dz);
480 ppp = p.GetVertexFast(i);
488 for (
int iface=1; iface <= p.GetNoFacets(); iface++) {
493 p.GetFacet(iface, nnode, iNodes, iVis, iFaces);
494 for (i=0; i<nnode; i++) {
498 if (iNodes[i] < 1 || iNodes[i] > p.GetNoVertices()) {
502 <<
"BooleanProcessor::takePolyhedron : problem 1."
506 if (iFaces[i] < 1 || iFaces[i] > p.GetNoFacets()) {
510 <<
"BooleanProcessor::takePolyhedron : problem 2."
520 iNodes[nnode] = iNodes[0];
522 for (i=0; i<nnode; i++) {
525 iface+dface, iFaces[i], iVis[i]));
533 for (i=0; i<3; i++) {
537 for (i=1; i<nnode; i++) {
539 for (k=0; k<3; k++) {
553 for (i=0; i<nnode; i++) { point +=
m_nodes[iNodes[i]].v; }
554 if (nnode > 1) point *= (1./nnode);
580 double rmin1[3], rmax1[3];
581 double rmin2[3], rmax2[3];
585 for (i=0; i<3; i++) {
595 for (i=0; i<3; i++) {
596 if (rmin1[i] > face.
rmin[i]) rmin1[i] = face.
rmin[i];
597 if (rmax1[i] < face.
rmax[i]) rmax1[i] = face.
rmax[i];
605 for (i=0; i<3; i++) {
606 if (rmin2[i] > face.
rmin[i]) rmin2[i] = face.
rmin[i];
607 if (rmax2[i] < face.
rmax[i]) rmax2[i] = face.
rmax[i];
614 for (i=0; i<3; i++) {
615 m_rmin[i] = (rmin1[i] > rmin2[i]) ? rmin1[i] : rmin2[i];
616 m_rmax[i] = (rmax1[i] < rmax2[i]) ? rmax1[i] : rmax2[i];
623 for (i=0; i<3; i++) {
624 if ((rmax1[i]-rmin1[i]) > del1) del1 = rmax1[i]-rmin1[i];
625 if ((rmax2[i]-rmin2[i]) > del2) del2 = rmax2[i]-rmin2[i];
627 return ((del1 < del2) ? del1 : del2) /
GRANULARITY;
640 int i, outflag, iface = ifaces, *prev;
656 for (i=0; i<3; i++) {
664 int npos = 0, nneg = 0;
666 for (i=0; i<8; i++) {
668 if (d > +
m_del) npos++;
669 if (d < -
m_del) nneg++;
671 if (npos == 8 || nneg == 8) outflag = 1;
700 int i, nnode, npos = 0, nneg = 0, nzer = 0;
705 nnode = (
m_faces[iface].iedges[3] == 0) ? 3 : 4;
706 for (i=0; i<nnode; i++) {
710 }
else if (dd[i] < -
m_del) {
719 if (npos == nnode || nneg == nnode)
return OUT_OF_PLANE;
727 int ie1 = 0, ie2 = 0, s1 = 0, s2 = 0, status, nint = 0;
728 enum { PLUS_MINUS, MINUS_PLUS, ZERO_ZERO, ZERO_PLUS, ZERO_MINUS };
731 for (i=0; i<nnode; i++) {
733 if (dd[i+1] >= 0)
continue;
735 }
else if (dd[i] < 0) {
736 if (dd[i+1] <= 0)
continue;
740 if (dd[i+1] > 0) status = ZERO_PLUS;
741 if (dd[i+1] < 0) status = ZERO_MINUS;
745 ie1 = i; s1 = status; nint++;
break;
747 ie2 = i; s2 = status; nint++;
break;
756 if (s1 != ZERO_ZERO && s2 != ZERO_ZERO) {
758 int iedge, i1 = 0, i2 = 0, ii[2];
759 double d1 = 0., d2 = 0., d3 = 0.;
760 ii[0] = ie1; ii[1] = ie2;
761 for (i=0; i<2; i++) {
762 iedge =
m_faces[iface].iedges[ii[i]];
771 }
else if (d1 < -
m_del) {
778 if (ii[i] == (
int)
m_nodes.size()) {
781 std::cerr<<
"d3 is zero in "<<__FILE__<<std::endl;
790 if (s1 == MINUS_PLUS || s1 == ZERO_PLUS) {
800 edge.
inext = (s1 == ZERO_ZERO) ? ie1+1 : ie2+1;
801 if (s1 == ZERO_PLUS || s2 == ZERO_MINUS) {
826 if (i1 == i2)
return;
827 if (
m_nodes[i1].s == 0 ||
m_nodes.back().s == 0) { i1 = i2;
return; }
830 if (i1 == ilast) { i1 = i2;
m_nodes.pop_back();
return; }
831 if (i2 == ilast) { i2 = i1; }
832 if (i3 == ilast) { i3 = i1; }
833 if (i4 == ilast) { i4 = i1; }
850 for (i=0; i<3; i++) {
853 if (d > dd) { dd = d; ii = i; }
859 if (t2-t1 < 0.) { t1 = -t1; t2 = -t2; t3 = -t3; t4 = -t4; }
861 if (t3 <= t1+m_del || t4 >= t2-
m_del)
return 0;
864 }
else if (t3 < t2-
m_del) {
869 }
else if (t4 > t1+
m_del) {
891 if (i1 < i2) { i2 = i1; }
892 else if (i1 > i2) { i1 = i2; }
893 else { iedges[1] = 0; }
894 if (iedges[0] == iedges[1])
return;
896 int ie1, ie2,
inode = i1;
898 for (
int i=0; i<2; i++) {
902 if ((ie1 = iedges[i]) == 0)
continue;
971 <<
"BooleanProcessor::caseIE : unimplemented case"
989 <<
"BooleanProcessor::caseEE : unimplemented case"
1011 for (
int i=0; i<3; i++) {
1067 int iedge =
m_faces[iface].inew;
1106#define INSERT_EDGE_TO_THE_LIST(A) \
1107*ilink = A; ilink = &m_edges[A].inext; *ilink = 0
1112 if (face.
inew == 0)
break;
1128 if (edge_i.
i1 == edge_cur.
i2)
break;
1135 if (edge_i.
i1 == edge_cur.
i2) {
1146 if (
m_edges[icur].i2 == ifirst) {
break; }
else {
continue; }
1151 <<
"BooleanProcessor::assembleFace(" << iface <<
") : "
1152 <<
"could not find next edge of the contour"
1165 if (what == 0 && ioldflag == 0 && iedge > 0) {
1167 if (
m_edges[iedge].inext > 0) {
1189 iface2 =
m_edges[iedge].iface2;
1208 if (
m_faces[iface].inew > 0) {
1333 int i, iedge, iface;
1342 iface =
m_edges[iedge].iface2;
1363 iface =
m_edges[iedge].iface2;
1406 double &
a,
double &b,
double &c)
const
1420 w = ::fabs(
a)+::fabs(b);
1436 double a1, b1, c1, a2, b2, c2, d1, d2;
1442 d1 = a1*
x[4] + b1*
y[4] + c1;
1443 d2 = a2*
x[4] + b2*
y[4] + c2;
1446 if ( a1*
x[2] + b1*
y[2] + c1 >= -
m_del)
return 1;
1453 d1 = a1*
x[1] + b1*
y[1] + c1;
1454 d2 = a2*
x[1] + b2*
y[1] + c2;
1457 if ( a1*
x[5] + b1*
y[5] + c1 >= -
m_del)
return 1;
1474 double x1, y1, x2, y2, a1, b1, c1;
1479 findABC(x1, y1, x2, y2, a1, b1, c1);
1483 int icontour, iedge, k1, k2;
1484 double x3, y3, x4, y4, a2, b2, c2, d1, d2;
1491 if (k1 == i1 || k2 == i1)
continue;
1492 if (k1 == i2 || k2 == i2)
continue;
1497 d1 = a1*x3 + b1*y3 + c1;
1498 d2 = a1*x4 + b1*y4 + c1;
1502 findABC(x3, y3, x4, y4, a2, b2, c2);
1503 d1 = a2*x1 + b2*y1 + c2;
1504 d2 = a2*x2 + b2*y2 + c2;
1519 if (k1 == i1 || k2 == i1)
continue;
1520 if (k1 == i2 || k2 == i2)
continue;
1525 d1 = a1*x3 + b1*y3 + c1;
1526 d2 = a1*x4 + b1*y4 + c1;
1530 findABC(x3, y3, x4, y4, a2, b2, c2);
1531 d1 = a2*x1 + b2*y1 + c2;
1532 d2 = a2*x2 + b2*y2 + c2;
1551 int i1ext, i2ext, i1int, i2int, i, k[6];
1563 for (i=0; i<3; i++) {
1577 for (i=3; i<6; i++) {
1630 for (
int i=0; i<3; i++) {
1638 if (a1*
x[1]+b1*
y[1]+c1 <= 0.1*
m_del)
return 1;
1643 double a2, b2, c2, a3, b3, c3;
1650 if (
m_edges[iedge].inext == iedge1)
return 0;
1652 if (
inode == k[0])
continue;
1653 if (
inode == k[1])
continue;
1654 if (
inode == k[2])
continue;
1657 if (a1*
x[1]+b1*
y[1]+c1 < -0.1*
m_del)
continue;
1658 if (a2*
x[1]+b2*
y[1]+c2 < -0.1*
m_del)
continue;
1659 if (a3*
x[1]+b3*
y[1]+c3 < -0.1*
m_del)
continue;
1682 int ipnext = ihead, nnode = 1;
1684 if (
m_edges[ipnext].inext > 0) {
1685 ipnext =
m_edges[ipnext].inext;
1688 m_edges[ipnext].inext = ihead;
1698 int iedge1, iedge2, iedge3, istart = 0;
1700 iedge1 =
m_edges[ipnext].inext;
1701 iedge2 =
m_edges[iedge1].inext;
1713 iedge3 =
m_edges[iedge2].inext;
1726 }
else if (istart == iedge1) {
1730 <<
"BooleanProcessor::triangulateContour : "
1731 <<
"could not generate a triangle (infinite loop)"
1740 ipnext =
m_edges[ipnext].inext;
1748 int iface1 =
m_edges[iedge1].iface1;
1762 m_edges[iedge1].iface1 = iface2;
1763 m_edges[iedge2].iface1 = iface2;
1764 ipnext =
m_edges[ipnext].inext;
1783 int iedge =
m_faces[iface].iold;
1794 <<
"BooleanProcessor::modifyReference : could not find the edge, "
1795 <<
"iface=" << iface <<
", i1,i2=" << i1 <<
"," << i2 <<
", iref=" << iref
1822 if (::fabs(normal[1]) > ::fabs(normal[iz])) iz = 1;
1823 if (::fabs(normal[2]) > ::fabs(normal[iz])) iz = 2;
1824 if (normal[iz] > 0) {
1825 ix = (iz+1)%3; iy = (ix+1)%3;
1827 iy = (iz+1)%3; ix = (iy+1)%3;
1835 int i1, i2, ifirst, iedge, icontour =
m_faces[iface].iold;
1836 while (icontour > 0) {
1846 z += node_1.
v[ix]*node_2.
v[iy]-node_2.
v[ix]*node_1.
v[iy];
1859 <<
"BooleanProcessor::triangulateFace : too small contour"
1863 icontour =
m_edges[iedge].inext;
1871 <<
"BooleanProcessor::triangulateFace : broken contour"
1892 <<
"BooleanProcessor::triangulateFace : "
1893 <<
"could not merge internal contour " << kint
1907 <<
"BooleanProcessor::triangulateFace : "
1908 <<
"triangulateContour failed."
1918 for (
int ifa=nface; ifa<(int)
m_faces.size(); ifa++) {
1921 if (
m_edges[iedge].iface1 != ifa) {
1925 <<
"BooleanProcessor::triangulateFace : wrong reference to itself, "
1926 <<
"iface=" << ifa <<
", iface1=" <<
m_edges[iedge].iface1
1929 }
else if (
m_edges[iedge].iface2 > 0) {
1932 }
else if (
m_edges[iedge].iface2 < 0) {
1950 int i, iedge, nnode = 0, nface = 0;
1956 for (i=1; i<(int)
m_faces.size(); i++) {
1969 for (i=1; i<(int)
m_nodes.size(); i++) {
1976 if (nface == 0)
return polyhedron;
1977 polyhedron.AllocateMemory(nnode, nface);
1981 for (i=1; i<(int)
m_nodes.size(); i++) {
1987 int k, v[4]={0}, f[4]={0};
1988 for (i=1; i<(int)
m_faces.size(); i++) {
1989 if (
m_faces[i].inew == 0)
continue;
1990 v[3] = f[3] = k = 0;
1994 std::cerr <<
"BooleanProcessor::createPolyhedron : too many edges" << std::endl;
1998 if (
m_edges[iedge].ivis < 0) v[k] = -v[k];
2004 std::cerr <<
"BooleanProcessor::createPolyhedron : "
2005 <<
"face has only " << k <<
" edges"
2008 polyhedron.m_pF[
m_faces[i].inew] =
2009 G4Facet(v[0],f[0], v[1],f[1], v[2],f[2], v[3],f[3]);
2015int BooleanProcessor::s_ishift = 0;
2069 <<
"BooleanProcessor: corrupted input polyhedron"
2081 <<
"BooleanProcessor: intersection with empty polyhedron"
2086 <<
"BooleanProcessor: subtraction from empty polyhedron"
2098 <<
"BooleanProcessor: intersection with empty polyhedron"
2118#define PROCESSOR_ERROR(a_what)
2120 unsigned int try_count = 1;
2196 <<
"BooleanProcessor::execute : unknown faces !!!"
2219 <<
"BooleanProcessor::execute : had converged."
2242 <<
"BooleanProcessor::execute : try another tilt..."
2249#undef PROCESSOR_ERROR
2478 std::cerr <<
"nodes : " <<
number << std::endl;
2481 std::cerr <<
" " <<
index
2482 <<
" x = " <<
node.v[0]
2483 <<
" y = " <<
node.v[1]
2484 <<
" z = " <<
node.v[2]
#define INSERT_EDGE_TO_THE_LIST(A)
#define PROCESSOR_ERROR(a_what)
HEPVis::SbPlane HVPlane3D
std::vector< int > m_internal_contours
std::vector< int > m_external_contours
void caseII(ExtEdge &edge1, ExtEdge &edge2)
void insertEdge(const ExtEdge &edge)
void draw_contour(int, int, int)
std::vector< ExtNode > m_nodes
FaceList m_unsuitable_faces
HepPolyhedron createPolyhedron()
void modifyReference(int iface, int i1, int i2, int iref)
void invertNewEdges(int iface)
std::vector< ExtFace > m_faces
void findABC(double x1, double y1, double x2, double y2, double &a, double &b, double &c) const
void draw_faces(int, int, int)
void triangulateFace(int iface)
void takePolyhedron(const HepPolyhedron &p, double, double, double)
void caseIE(ExtEdge &edge1, ExtEdge &edge2)
void assembleNewFaces(int what, int ihead)
int checkDirection(double *x, double *y) const
void selectOutsideFaces(int &ifaces, int &iout)
int testEdgeVsEdge(ExtEdge &edge1, ExtEdge &edge2)
void mergeContours(int ix, int iy, int kext, int kint)
void checkDoubleEdges(int iface)
HepPolyhedron execute(int op, const HepPolyhedron &a, const HepPolyhedron &b, int &err)
void caseEE(ExtEdge &edge1, ExtEdge &edge2)
FaceList m_suitable_faces
void renumberNodes(int &i1, int &i2, int &i3, int &i4)
static void set_shift(int)
std::vector< ExtEdge > m_edges
void triangulateContour(int ix, int iy, int ihead)
int checkTriangle(int iedge1, int iedge2, int ix, int iy) const
void assembleFace(int what, int iface)
int checkIntersection(int ix, int iy, int i1, int i2) const
int testFaceVsPlane(ExtEdge &edge)
void testFaceVsFace(int iface1, int iface2)
void divideEdge(int &i1, int &i2)
static int get_num_shift()
int get_processor_error() const
ExtEdge(const ExtEdge &edge)
ExtEdge(int k1=0, int k2=0, int kface1=0, int kface2=0, int kvis=0)
ExtEdge & operator=(const ExtEdge &edge)
ExtFace & operator=(const ExtFace &face)
ExtFace(std::vector< ExtEdge > &a_edges, int iedge)
ExtFace(const ExtFace &face)
std::vector< ExtEdge > & m_edges
ExtNode & operator=(const ExtNode &node)
ExtNode(const ExtNode &node)
ExtNode(HVPoint3D vertex=HVPoint3D(), int status=0)
friend class BooleanProcessor
FaceList(std::vector< ExtFace > &a_faces)
std::vector< ExtFace > & m_faces
double distance(const SbVec3d &point) const
std::string number(const double &d, const std::string &s)