183 #ifdef BUILD_OCTREE_DLL
184 #define EGS_OCTREE_EXPORT __declspec(dllexport)
186 #define EGS_OCTREE_EXPORT __declspec(dllimport)
188 #define EGS_OCTREE_LOCAL
192 #ifdef HAVE_VISIBILITY
193 #define EGS_OCTREE_EXPORT __attribute__ ((visibility ("default")))
194 #define EGS_OCTREE_LOCAL __attribute__ ((visibility ("hidden")))
196 #define EGS_OCTREE_EXPORT
197 #define EGS_OCTREE_LOCAL
207 int ixmin, iymin, izmin;
208 int ixmax, iymax, izmax;
213 vmin(boxMin), vmax(boxMax), nx(bboxRes[0]), ny(bboxRes[1]), nz(bboxRes[2]) {
215 ixmin = iymin = izmin = 0;
216 ixmax = iymax = izmax = 0;
225 if (b2.vmin.
x < b1.vmin.
x) {
226 b1.vmin.
x = b2.vmin.
x;
228 if (b2.vmin.
y < b1.vmin.
y) {
229 b1.vmin.
y = b2.vmin.
y;
231 if (b2.vmin.
z < b1.vmin.
z) {
232 b1.vmin.
z = b2.vmin.
z;
234 if (b2.vmax.
x > b1.vmax.
x) {
235 b1.vmax.
x = b2.vmax.
x;
237 if (b2.vmax.
y > b1.vmax.
y) {
238 b1.vmax.
y = b2.vmax.
y;
240 if (b2.vmax.
z > b1.vmax.
z) {
241 b1.vmax.
z = b2.vmax.
z;
245 EGS_Float dx1 = (b1.vmax.
x-b1.vmin.
x)/b1.nx;
246 EGS_Float dy1 = (b1.vmax.
y-b1.vmin.
y)/b1.ny;
247 EGS_Float dz1 = (b1.vmax.
z-b1.vmin.
z)/b1.nz;
248 EGS_Float dx2 = (b2.vmax.
x-b2.vmin.
x)/b2.nx;
249 EGS_Float dy2 = (b2.vmax.
y-b2.vmin.
y)/b2.ny;
250 EGS_Float dz2 = (b2.vmax.
z-b2.vmin.
z)/b2.nz;
254 b1.nx = (int)((b1.vmax.
x - b1.vmin.
x + 0.5*dx2) / dx2);
257 b1.nx = (int)((b1.vmax.
x - b1.vmin.
x +0.5*dx1) / dx1);
260 b1.ny = (int)((b1.vmax.
y - b1.vmin.
y + 0.5*dy2) / dy2);
263 b1.ny = (int)((b1.vmax.
y - b1.vmin.
y +0.5*dy1) / dy1);
266 b1.nz = (int)((b1.vmax.
z - b1.vmin.
z + 0.5*dz2) / dz2);
269 b1.nz = (int)((b1.vmax.
z - b1.vmin.
z +0.5*dz1) / dz1);
279 if (b2.vmin.
x < vmin.
x) {
282 if (b2.vmin.
y < vmin.
y) {
285 if (b2.vmin.
z < vmin.
z) {
288 if (b2.vmax.
x > vmax.
x) {
291 if (b2.vmax.
y > vmax.
y) {
294 if (b2.vmax.
z > vmax.
z) {
299 EGS_Float dx1 = (vmax.
x-vmin.
x)/nx;
300 EGS_Float dy1 = (vmax.
y-vmin.
y)/ny;
301 EGS_Float dz1 = (vmax.
z-vmin.
z)/nz;
302 EGS_Float dx2 = (b2.vmax.
x-b2.vmin.
x)/b2.nx;
303 EGS_Float dy2 = (b2.vmax.
y-b2.vmin.
y)/b2.ny;
304 EGS_Float dz2 = (b2.vmax.
z-b2.vmin.
z)/b2.nz;
308 nx = (int)((vmax.
x - vmin.
x + 0.5*dx2) / dx2);
311 nx = (int)((vmax.
x - vmin.
x +0.5*dx1) / dx1);
314 ny = (int)((vmax.
y - vmin.
y + 0.5*dy2) / dy2);
317 ny = (int)((vmax.
y - vmin.
y +0.5*dy1) / dy1);
320 nz = (int)((vmax.
z - vmin.
z + 0.5*dz2) / dz2);
323 nz = (int)((vmax.
z - vmin.
z +0.5*dz1) / dz1);
352 void createChildren() {
356 egsFatal(
"EGS_Octree_node::createChildren(): Memory allocation error");
358 for (
int i=0; i<8; i++) {
360 child[i].ix = (ix << 1) | (i>>0 & 0x1);
361 child[i].iy = (iy << 1) | (i>>1 & 0x1);
362 child[i].
iz = (iz << 1) | (i>>2 & 0x1);
369 void deleteChildren() {
371 for (
int i=0; i<8; i++) {
372 child[i].deleteChildren();
380 int collapseChildren() {
382 for (
int i=0; i<8; i++) {
399 int shift = bbox.maxlevel -
level;
402 if (ii < bbox.ixmin) {
406 if (ii > bbox.ixmax) {
410 if (ii < bbox.iymin) {
414 if (ii > bbox.iymax) {
418 if (ii < bbox.izmin) {
422 if (ii > bbox.izmax) {
431 int shift = bbox.maxlevel -
level;
435 if (bbox.ixmin > iimin) {
438 iimax = ~(~ix<<shift);
439 if (bbox.ixmax < iimax) {
447 if (bbox.iymin > iimin) {
450 iimax = ~(~iy<<shift);
451 if (bbox.iymax < iimax) {
459 if (bbox.izmin > iimin) {
462 iimax = ~(~
iz<<shift);
463 if (bbox.izmax < iimax) {
523 EGS_Float bbxmin, bbymin, bbzmin;
524 EGS_Float bbxmax, bbymax, bbzmax;
525 EGS_Float xmin, ymin, zmin;
526 EGS_Float xmax, ymax, zmax;
527 EGS_Float dx, dy, dz, dxi, dyi, dzi;
528 int ixmin, iymin, izmin;
529 int ixmax, iymax, izmax;
533 long int nLeaf, nLeafMax;
534 vector<EGS_Octree_node *> tmp;
542 for (
int i=1; i<vBox.size(); i++) {
545 vBox.push_back(bbox);
553 dx = (bbox.vmax.
x - bbox.vmin.
x) / bbox.nx;
554 dy = (bbox.vmax.
y - bbox.vmin.
y) / bbox.ny;
555 dz = (bbox.vmax.
z - bbox.vmin.
z) / bbox.nz;
570 maxlevel = (int) ceil(log((EGS_Float)res)/0.6931471805599452862);
579 xmax = bbox.vmin.
x + dx*n;
580 ymax = bbox.vmin.
y + dy*n;
581 zmax = bbox.vmin.
z + dz*n;
585 for (
int i=0; i<vBox.size(); i++) {
587 box->ixmin = (int)((box->vmin.
x - xmin + 0.5*dx) * dxi);
588 box->iymin = (int)((box->vmin.
y - ymin + 0.5*dy) * dyi);
589 box->izmin = (int)((box->vmin.
z - zmin + 0.5*dz) * dzi);
590 box->nx = (int)((box->vmax.
x - box->vmin.
x + 0.5*dx) * dxi);
591 box->ny = (int)((box->vmax.
y - box->vmin.
y + 0.5*dy) * dyi);
592 box->nz = (int)((box->vmax.
z - box->vmin.
z + 0.5*dz) * dzi);
593 box->ixmax = box->ixmin + box->nx - 1;
594 box->iymax = box->iymin + box->ny - 1;
595 box->izmax = box->izmin + box->nz - 1;
614 box->nx = box->ixmax - box->ixmin + 1;
615 box->ny = box->iymax - box->iymin + 1;
616 box->nz = box->izmax - box->izmin + 1;
634 bbxmin = box->vmin.
x;
635 bbxmax = box->vmax.
x;
636 bbymin = box->vmin.
y;
637 bbymax = box->vmax.
y;
638 bbzmin = box->vmin.
z;
639 bbzmax = box->vmax.
z;
649 growOctree(root, vBox, pruneTree);
655 for (
int i=0; i<
nreg; i++) {
659 tmp.erase(tmp.begin(),tmp.end());
664 statOctree(root, vBox);
671 root->deleteChildren();
679 void statOctree(
EGS_Octree_node *node, vector<EGS_Octree_bbox> &vBox) {
681 for (
int i=0; i<8; i++) {
682 statOctree(node->
child+i, vBox);
685 else if (node->insideBBox(vBox[vBox.size()-1]) && node->
medium>=0) {
686 int shift = maxlevel - node->
level;
688 nLeafMax += 1<<(3*shift);
694 void growOctree(
EGS_Octree_node *node, vector<EGS_Octree_bbox> &vBox,
bool prune) {
697 bool refineNode =
false;
700 if (node->
level < maxlevel) {
705 for (
int i=0; i<vBox.size()-1; i++) {
706 int boxlevel = vBox[i].level;
707 if (node->insideBBox(vBox[i])) {
715 for (
int i=0; i<vBox.size()-1; i++) {
716 int boxlevel = vBox[i].level;
717 if (!node->insideBBox(vBox[i]) && node->intersectBBox(vBox[i])) {
718 if (boxlevel > level) {
726 if (node->
level < level) {
733 node->createChildren();
734 for (
int i=0; i<8; i++) {
735 growOctree(node->
child+i, vBox, prune);
738 if (node->collapseChildren()) {
739 tmp.resize(tmp.size()-8);
740 node->
region = tmp.size();
747 bool insideSomeBox =
false;
748 for (
int i=0; i<vBox.size()-1; i++) {
749 if (node->insideBBox(vBox[i])) {
750 insideSomeBox =
true;
756 int shift = maxlevel - node->
level;
757 int ix = node->ix << shift;
758 int iy = node->iy << shift;
759 int iz = node->
iz << shift;
760 int ireg = geom->isWhere(
EGS_Vector(xmin+(ix+0.5)*dx, ymin+(iy+0.5)*dy, zmin+(iz+0.5)*dz));
762 node->
medium = geom->medium(ireg);
765 node->
region = tmp.size();
775 if (ixn < ixmin || ixn > ixmax) {
780 int shift = maxlevel - node->
level;
781 if (iyn < node->iy<<shift) {
782 iyn = node->iy<<shift;
784 else if (iyn > ~(~node->iy<<shift)) {
785 iyn = ~(~node->iy<<shift);
787 if (izn < node->iz<<shift) {
788 izn = node->
iz<<shift;
790 else if (izn > ~(~node->
iz<<shift)) {
791 izn = ~(~node->
iz<<shift);
795 int diff = node->ix ^ (ixn>>shift);
797 while (diff & (1<<shift)) {
801 shift = maxlevel - node->
level;
802 while (node->
child) {
805 childIndex = (ixn>>shift & 0x1);
806 childIndex |= (iyn>>shift & 0x1) << 1;
807 childIndex |= (izn>>shift & 0x1) << 2;
808 node = node->
child + childIndex;
810 shift = maxlevel - node->
level;
820 if (iyn < iymin || iyn > iymax) {
825 int shift = maxlevel - node->
level;
826 if (ixn < node->ix<<shift) {
827 ixn = node->ix<<shift;
829 else if (ixn > ~(~node->ix<<shift)) {
830 ixn = ~(~node->ix<<shift);
832 if (izn < node->iz<<shift) {
833 izn = node->
iz<<shift;
835 else if (izn > ~(~node->
iz<<shift)) {
836 izn = ~(~node->
iz<<shift);
840 int diff = node->iy ^ (iyn>>shift);
842 while (diff & (1<<shift)) {
846 shift = maxlevel - node->
level;
847 while (node->
child) {
850 childIndex = (ixn>>shift & 0x1);
851 childIndex |= (iyn>>shift & 0x1) << 1;
852 childIndex |= (izn>>shift & 0x1) << 2;
853 node = node->
child + childIndex;
863 if (izn < izmin || izn > izmax) {
868 int shift = maxlevel - node->
level;
869 if (ixn < node->ix<<shift) {
870 ixn = node->ix<<shift;
872 else if (ixn > ~(~node->ix<<shift)) {
873 ixn = ~(~node->ix<<shift);
875 if (iyn < node->iy<<shift) {
876 iyn = node->iy<<shift;
878 else if (iyn > ~(~node->iy<<shift)) {
879 iyn = ~(~node->iy<<shift);
883 int diff = node->
iz ^ (izn>>shift);
885 while (diff & (1<<shift)) {
889 shift = maxlevel - node->
level;
890 while (node->
child) {
893 childIndex = (ixn>>shift & 0x1);
894 childIndex |= (iyn>>shift & 0x1) << 1;
895 childIndex |= (izn>>shift & 0x1) << 2;
896 node = node->
child + childIndex;
904 if ((ix<ixmin) || (ix>ixmax) || (iy<iymin) || (iy>iymax) || (iz<izmin) || (iz>izmax)) {
908 int shift = maxlevel;
909 while (node->
child) {
912 childIndex = (ix>>shift & 0x1);
913 childIndex |= (iy>>shift & 0x1) << 1;
914 childIndex |= (iz>>shift & 0x1) << 2;
915 node = node->
child + childIndex;
922 void setIndices(
const EGS_Vector &r,
int &ix,
int &iy,
int &iz) {
923 ix = (int)((r.
x-xmin)*dxi);
924 iy = (int)((r.
y-ymin)*dyi);
925 iz = (int)((r.
z-zmin)*dzi);
949 if (r.
x >= bbxmin && r.
x <= bbxmax &&
950 r.
y >= bbymin && r.
y <= bbymax &&
951 r.
z >= bbzmin && r.
z <= bbzmax) {
964 setIndices(r,ix,iy,iz);
965 return getNode(ix,iy,iz)->region;
981 setIndices(r, ix, iy, iz);
982 return getNode(ix,iy,iz)->region;
987 int medium(
int ireg)
const {
988 return nodeReg[ireg]->medium;
1003 int shift = maxlevel - node->
level;
1016 xBound = xmin + (ix+1)*dx;
1018 EGS_Float d = (xBound-r.
x) / u.
x;
1035 xBound = xmin + ix*dx;
1037 EGS_Float d = (xBound-r.
x) / u.
x;
1057 yBound = ymin + (iy+1)*dy;
1059 EGS_Float d = (yBound - r.
y) / u.
y;
1076 yBound = ymin + iy*dy;
1078 EGS_Float d = (yBound-r.
y) / u.
y;
1097 zBound = zmin + (iz+1)*dz;
1099 EGS_Float d = (zBound-r.
z) / u.
z;
1116 zBound = zmin + iz*dz;
1118 EGS_Float d = (zBound-r.
z) / u.
z;
1135 setIndices(ryz, tmp, iy, iz);
1136 node = getNeighborNodeX(node, ix, iy, iz);
1138 else if (crossed==1) {
1140 setIndices(rxz, ix, tmp, iz);
1141 node = getNeighborNodeY(node, ix, iy, iz);
1143 else if (crossed==2) {
1145 setIndices(rxy, ix, iy, tmp);
1146 node = getNeighborNodeZ(node, ix, iy, iz);
1162 int ix=0, iy=0, iz=0;
1163 EGS_Float d, tlong = 2*t;
1167 if (r.
x <= bbxmin && u.
x > 0) {
1169 d = (bbxmin-r.
x) / u.
x;
1171 else if (r.
x >= bbxmax && u.
x < 0) {
1173 d = (bbxmax-r.
x) / u.
x;
1179 EGS_Float yy = r.
y + u.
y*d;
1180 EGS_Float zz = r.
z + u.
z*d;
1181 if (yy >= bbymin && yy <= bbymax && zz >= bbzmin && zz <= bbzmax) {
1183 setIndices(rr, tmp, iy, iz);
1188 node = getNode(ix,iy,iz);
1194 if (r.
y <= bbymin && u.
y > 0) {
1196 d = (bbymin-r.
y) / u.
y;
1198 else if (r.
y >= bbymax && u.
y < 0) {
1200 d = (bbymax-r.
y) / u.
y;
1206 EGS_Float xx = r.
x + u.
x*d;
1207 EGS_Float zz = r.
z + u.
z*d;
1208 if (xx >= bbxmin && xx <= bbxmax && zz >= bbzmin && zz <= bbzmax) {
1210 setIndices(rr, ix, tmp, iz);
1215 node = getNode(ix,iy,iz);
1221 if (r.
z <= bbzmin && u.
z > 0) {
1223 d = (bbzmin-r.
z) / u.
z;
1225 else if (r.
z >= bbzmax && u.
z < 0) {
1227 d = (bbzmax-r.
z) / u.
z;
1233 EGS_Float xx = r.
x + u.
x*d;
1234 EGS_Float yy = r.
y + u.
y*d;
1235 if (xx >= bbxmin && xx <= bbxmax && yy >= bbymin && yy <= bbymax) {
1237 setIndices(rr, ix, iy, tmp);
1242 node = getNode(ix,iy,iz);
1258 inew = howfarOut(r, u, t, normal);
1261 inew = howfarIn(nodeReg[ireg], r, u, t, normal);
1265 if (inew>=0 && newmed) {
1266 *newmed = nodeReg[inew]->medium;
1273 EGS_Float hownearIn(
int ireg,
const EGS_Vector &r) {
1275 int shift = maxlevel - node->
level;
1276 EGS_Float t1, t2, tx, ty, tz;
1280 imin = node->ix << shift;
1281 imax = ~(~node->ix << shift);
1282 t1 = (r.
x-xmin)-dx*imin;
1283 t2 = dx*(imax-imin+1)-t1;
1284 tx = t1 < t2 ? t1 : t2;
1287 imin = node->iy << shift;
1288 imax = ~(~node->iy << shift);
1289 t1 = (r.
y-ymin)-dy*imin;
1290 t2 = dy*(imax-imin+1)-t1;
1291 ty = t1 < t2 ? t1 : t2;
1294 imin = node->
iz << shift;
1295 imax = ~(~node->
iz << shift);
1296 t1 = (r.
z-zmin)-dz*imin;
1297 t2 = dz*(imax-imin+1)-t1;
1298 tz = t1 < t2 ? t1 : t2;
1300 return tx<ty && tx<tz ? tx : ty<tz ? ty : tz;
1307 return hownearIn(ireg, r);
1310 EGS_Float s1=0, s2=0;
1311 if (r.
x < bbxmin || r.
x > bbxmax) {
1312 EGS_Float t = r.
x < bbxmin ? bbxmin-r.
x : r.
x-bbxmax;
1317 if (r.
y < bbymin || r.
y > bbymax) {
1318 EGS_Float t = r.
y < bbymin ? bbymin-r.
y : r.
y-bbymax;
1323 if (r.
z < bbzmin || r.
z > bbzmax) {
1324 EGS_Float t = r.
z < bbzmin ? bbzmin-r.
z : r.
z-bbzmax;
1329 return nc == 1 ? s1 : sqrt(s2);
1334 const string &
getType()
const {
int iz
octree indices, representing the binary location code of the node in x, y, z;
virtual const string & getType() const =0
Get the geometry type.
virtual bool isInside(const EGS_Vector &x)=0
Is the position x inside the geometry?
virtual int medium(int ireg) const
Returns the medium index in region ireg.
virtual void printInfo() const
Print information about this geometry.
EGS_BaseGeometry(const string &Name)
Construct a geometry named Name.
EGS_Octree_node * child
pointer to children nodes, NULL is there are no children
A class representing 3D vectors.
Global egspp functions header file.
virtual int isWhere(const EGS_Vector &x)=0
In which region is poisition x?
virtual int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u, EGS_Float &t, int *newmed=0, EGS_Vector *normal=0)=0
Calculate the distance to a boundary from x along the direction u.
short medium
medium index for the node
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
virtual int inside(const EGS_Vector &x)=0
Returns the region index, if inside, or -1 if outside (obsolete)
EGS_Octree_node * parent
pointer to the parent node (only root node can have parent set to NULL)
virtual EGS_Float hownear(int ireg, const EGS_Vector &x)=0
Calculate the distance to a boundary for position x in any direction.
unsigned short level
depth of the node (root node is level 0)
int region
region number of the node
EGS_BaseGeometry class header file.
int nreg
Number of local regions in this geometry.