30 template <
int tmplDim,
typename TmplVec>
37 template <
int tmplDim,
typename TmplVec>
43 template <
int tmplDim,
typename TmplVec>
49 template <
int tmplDim,
typename TmplVec>
53 for (
int i = 0; i < tmplDim; i++)
55 product *= (this->getMaxPt()[i] - this->getMinPt()[i]);
60 template <
int tmplDim,
typename TmplVec>
61 template <
typename TmplSphere>
64 double distSqrd = 0.0;
65 for (
int i = 0; i < tmplDim; i++)
67 if (sphere.getCentre()[i] < this->getMinPt()[i])
69 distSqrd +=
square(sphere.getCentre()[i] - this->getMinPt()[i]);
71 else if (sphere.getCentre()[i] > this->getMaxPt()[i])
73 distSqrd +=
square(sphere.getCentre()[i] - this->getMaxPt()[i]);
76 return (distSqrd <=
square(sphere.getRadius()));
79 template <
int tmplDim,
typename TmplVec>
82 for (
int i = 0; (i < tmplDim); i++)
84 if ((this->getMinPt()[i] > pt[i]) || (pt[i] > this->getMaxPt()[i]))
92 template <
int tmplDim,
typename TmplVec>
93 template <
typename TmplSphere>
96 for (
int i = 0; (i < tmplDim); i++)
99 pt[i] = sphere.getRadius();
102 !(intersectsWith(sphere.getCentre() + pt))
104 !(intersectsWith(sphere.getCentre() - pt))
113 template <
int tmplDim,
typename TmplVec>
117 for (
int i = 0; i < tmplDim; i++)
124 template <
int tmplDim,
typename TmplVec>
128 for (
int i = 0; i < tmplDim; i++)
135 template <
int tmplDim,
typename TmplVec>
140 template <
int tmplDim,
typename TmplVec>
144 m_invNormalNorm((1.0/norm(normal)))
148 template <
int tmplDim,
typename TmplVec>
150 : m_normal(plane.m_normal),
152 m_invNormalNorm(plane.m_invNormalNorm)
156 template <
int tmplDim,
typename TmplVec>
165 template <
int tmplDim,
typename TmplVec>
171 (
dot(m_normal, pt) -
dot(m_normal, m_pt))*m_invNormalNorm
175 template <
int tmplDim,
typename TmplVec>
178 return fabs(getSignedDistanceTo(pt));
181 template <
int tmplDim,
typename TmplVec>
187 template <
int tmplDim,
typename TmplVec>
190 template <
int tmplDim,
typename TmplVec>
193 template <
int tmplDim,
typename TmplVec>
200 template <
int tmplDim,
typename TmplVec>
202 : m_centre(centrePt),
207 template <
int tmplDim,
typename TmplVec>
209 : m_centre(sphere.m_centre),
210 m_radius(sphere.m_radius)
214 template <
int tmplDim,
typename TmplVec>
222 template <
int tmplDim,
typename TmplVec>
228 template <
int tmplDim,
typename TmplVec>
235 template <
int tmplDim,
typename TmplVec>
238 return (tmplDim == 2) ? M_PI*getRadius()*getRadius() : FOUR_THIRDS_PI*getRadius()*getRadius()*getRadius();
241 inline void checkDomain(
double r,
double x1,
double y1,
double x2,
double y2)
243 const double rSqrd = r*r;
244 const double x1Sqrd = x1*x1;
245 const double x2Sqrd = x2*x2;
246 const double y1Sqrd = y1*y1;
247 const double y2Sqrd = y2*y2;
250 ((rSqrd - x1Sqrd - y1Sqrd) < 0)
252 ((rSqrd - x1Sqrd - y2Sqrd) < 0)
254 ((rSqrd - x2Sqrd - y1Sqrd) < 0)
256 ((rSqrd - x2Sqrd - y2Sqrd) < 0)
259 std::stringstream msg;
261 <<
"Invalid rectangular domain for sphere integration, points in domain "
262 <<
"(" << x1 <<
"," << y1 <<
"), (" << x2 <<
"," << y2 <<
" lie outside "
263 <<
"sphere of radius " << r <<
" centred at the origin.";
264 throw std::runtime_error(msg.str());
268 template <
int tmplDim,
typename TmplVec>
272 if ((tmplDim == 2) || (tmplDim == 3))
274 if (minPt[dimX] != maxPt[dimX])
276 const double x1 = minPt[dimX] - getCentre()[dimX];
277 const double x2 = maxPt[dimX] - getCentre()[dimX];
278 const double r = getRadius();
282 const double rSqrd = r*r;
283 const double x1Sqrd = x1*x1;
284 const double x2Sqrd = x2*x2;
290 x2*sqrt(rSqrd-x2Sqrd)
294 x1*sqrt(rSqrd-x1Sqrd)
297 else if (tmplDim == 3)
299 if (minPt[dimY] != maxPt[dimY])
301 const double y1 = minPt[dimY] - getCentre()[dimY];
302 const double y2 = maxPt[dimY] - getCentre()[dimY];
318 const double t30 = y2*y2;
319 const double t31 = x2*x2;
320 const double t36 = r*r;
321 const double t59 = t31-t36;
322 const double t40 = sqrt(-t30-t59);
323 const double t10 = 1.0/t40;
324 const double t32 = x1*x1;
325 const double t54 = t32-t36;
326 const double t42 = sqrt(-t30-t54);
327 const double t14 = 1.0/t42;
328 const double t64 = -atan(t10*x2)+atan(t14*x1);
329 const double t27 = y1*y1;
330 const double t39 = sqrt(-t27-t59);
331 const double t9 = 1.0/t39;
332 const double t41 = sqrt(-t27-t54);
333 const double t12 = 1.0/t41;
334 const double t63 = atan(t12*x1)-atan(t9*x2);
335 const double t62 = -atan(y2*t14)+atan(y1*t12);
336 const double t61 = -atan(y1*t9)+atan(t10*y2);
337 const double t37 = sqrt(t31);
338 const double t21 = 1.0/t37;
339 const double t60 = t21*t9;
340 const double t38 = sqrt(t32);
341 const double t24 = 1.0/t38;
342 const double t58 = t24*t14;
343 const double t57 = t24*t12;
344 const double t56 = t37*t38;
345 const double t55 = t21*t10;
346 const double t53 = 2.0*x2;
347 const double t52 = 2.0*x1;
348 const double t51 = t42*t56;
349 const double t28 = t27*y1;
350 const double t50 = t28-t36*y1;
351 const double t34 = t30*y2;
352 const double t49 = t34-t36*y2;
353 const double t48 = t41*t51;
354 const double t35 = t36*r;
355 const double t33 = t31*x2;
356 const double t29 = t32*x1;
357 const double t26 = r*y2;
358 const double t25 = y1*r;
363 (-2.0*t33*y1-t50*t53)*t40*t48
366 (2.0*t33*y2+t49*t53)*t48
369 (2.0*t29*y1+t50*t52)*t51
372 (-2.0*t29*y2-t49*t52)*t56
387 -atan((-t26+t59)*t55)
396 (-t64*t34+t61*t33+t62*t29+t63*t28+3.0*(t64*y2-t63*y1-t61*x2-t62*x1)*t36)*t37
411 template <
int tmplDim,
typename TmplVec>
414 double distSqrd = 0.0;
415 for (
int i = 0; i < tmplDim; i++)
417 distSqrd +=
square(getCentre()[i] - pt[i]);
419 return (distSqrd <=
square(getRadius()));
422 template <
int tmplDim,
typename TmplVec>
426 if ((tmplDim == 2) || (tmplDim == 3))
429 const double d = fabs(signedD);
435 const double rSqrd = getRadius()*getRadius();
436 vol = rSqrd*acos(d/getRadius()) - d*sqrt(rSqrd - d*d);
438 else if (tmplDim == 3)
441 const double h = getRadius() - d;
442 vol = ONE_THIRD_PI*h*h*(3.0*getRadius()-h);
444 vol = ((signedD < 0) ? vol : getVolume() - vol);
450 template <
int tmplDim,
typename TmplVec>
459 template <
int tmplDim,
typename TmplVec>
468 template <
int tmplDim,
typename TmplVec>
475 template <
int tmplDim,
typename TmplVec>
484 template <
int tmplDim,
typename TmplVec>
489 m_volume(sphere.m_volume)
493 template <
int tmplDim,
typename TmplVec>
504 template <
int tmplDim,
typename TmplVec>
510 template <
int tmplDim,
typename TmplVec>
517 template <
int tmplDim,
typename TmplVec>
523 template <
int tmplDim,
typename TmplVec>
534 template <
int tmplDim,
typename TmplVec>
540 template <
int tmplDim,
typename TmplVec>
546 template <
int tmplDim,
typename TmplVec>
552 template <
int tmplDim,
typename TmplVec>
557 template <
int tmplDim,
typename TmplVec>
562 template <
int tmplDim,
typename TmplVec>
567 template <
int tmplDim,
typename TmplVec>
575 template <
int tmplDim,
typename TmplVec>
582 template <
int tmplDim,
typename TmplVec>
588 template <
int tmplDim,
typename TmplVec>
598 template <
int tmplDim,
typename TmplVec>
610 template <
int tmplDim,
typename TmplVec>
616 BasicBox::operator=(box);
617 for (
int i = 0; i < getNumVertices(); i++)
624 template <
int tmplDim,
typename TmplVec>
628 m_vertexArray[j].setPoint(this->getMinPt());
630 for (j++; i < tmplDim; i++, j++)
632 Vec pt = this->getMinPt();
633 pt[i] = this->getMaxPt()[i];
634 m_vertexArray[j].setPoint(pt);
637 m_vertexArray[j].setPoint(this->getMaxPt());
638 for (i = 0, j++; i < tmplDim && j < s_numVertices; i++, j++)
640 Vec pt = this->getMaxPt();
641 pt[i] = this->getMinPt()[i];
642 m_vertexArray[j] = pt;
646 template <
int tmplDim,
typename TmplVec>
650 return m_vertexArray[i];
653 template <
int tmplDim,
typename TmplVec>
656 return s_numVertices;
659 template <
int tmplDim,
typename TmplVec>
668 template <
int tmplDim,
typename TmplVec>
675 template <
int tmplDim,
typename TmplVec>
683 template <
int tmplDim,
typename TmplVec>
690 template <
int tmplDim,
typename TmplVec>
697 template <
int tmplDim,
typename TmplVec>
705 for (
int i = 0; i < tmplDim; i++)
707 m[i] = min(p1[i], p2[i]);
712 template <
int tmplDim,
typename TmplVec>
720 for (
int i = 0; i < tmplDim; i++)
722 m[i] = max(p1[i], p2[i]);
727 template <
int tmplDim,
typename TmplVec>
733 const Vec ¢rePt = getSphere().getCentre();
734 const Vec diag = (getSphere().getCentre() - pt)*2.0;
735 const Vec oppCorner = pt + diag;
738 componentMin(pt, oppCorner),
739 componentMax(pt, oppCorner)
742 const double sphVol = getSphere().getVolume();
746 for (
int i = 0; i < tmplDim; i++)
748 s[i] = getSphere().getSegmentVolume(
Plane(getNormal((i+1)%tmplDim), box.
getMaxPt()));
752 v[0] = 0.50*(sphVol - 2.0*s[0] - boxVol);
753 v[1] = 0.50*(sphVol - 2.0*s[1] - boxVol);
754 v[2] = 0.25*(sphVol - 2.0*v[0] - 2.0*v[1] - boxVol);
756 if (pt[0] <= centrePt[0])
758 if (pt[1] <= centrePt[1])
760 vol = boxVol + v[0] + v[1] + v[2];
769 if (pt[1] <= centrePt[1])
779 else if (tmplDim == 3)
783 2.0*getSphere().getVolume(
794 2.0*getSphere().getVolume(
805 2.0*getSphere().getVolume(
816 e[0] = 0.250*(sphVol - 2.0*s[1] - 2.0*v[0] - 2.0*v[1] - boxVol);
817 e[1] = 0.250*(sphVol - 2.0*s[2] - 2.0*v[1] - 2.0*v[2] - boxVol);
818 e[2] = 0.250*(sphVol - 2.0*s[0] - 2.0*v[0] - 2.0*v[2] - boxVol);
824 2.0*v[0] - 2.0*v[1] - 2.0*v[2]
826 4.0*e[0] - 4.0*e[1] - 4.0*e[2]
831 if (pt[0] <= centrePt[0])
833 if (pt[1] <= centrePt[1])
835 if (pt[2] <= centrePt[2])
837 vol = boxVol + v[0] + v[1] + v[2] + v[3] + e[0] + e[1] + e[2];
841 vol = v[2] + v[3] + e[1] + e[2];
846 if (pt[2] <= centrePt[2])
848 vol = v[1] + v[3] + e[0] + e[1];
858 if (pt[1] <= centrePt[1])
860 if (pt[2] <= centrePt[2])
862 vol = v[0] + v[3] + e[0] + e[2];
871 if (pt[2] <= centrePt[2])
886 template <
int tmplDim,
typename TmplVec>
892 const int ZERO = orientDim;
893 const int ONE = (orientDim+1) % tmplDim;
894 const int TWO = (orientDim+2) % tmplDim;
897 const double sphVol = getSphere().getVolume();
898 const Vec ¢rePt = getSphere().getCentre();
900 if ((
square(pt[ONE]-centrePt[ONE]) +
square(pt[TWO]-centrePt[TWO])) <
square(getSphere().getRadius()))
902 Plane plane[tmplDim];
903 plane[ZERO] =
Plane();
904 plane[ONE] =
Plane(getNormal(ONE), pt);
905 plane[TWO] =
Plane(getNormal(TWO), pt);
907 const double halfSphVol = sphVol*0.5;
910 s[ONE] = getSphere().getSegmentVolume(plane[ONE]);
911 s[TWO] = getSphere().getSegmentVolume(plane[TWO]);
912 s[ONE] = ((s[ONE] > halfSphVol) ? (sphVol - s[ONE]) : s[ONE]);
913 s[TWO] = ((s[TWO] > halfSphVol) ? (sphVol - s[TWO]) : s[TWO]);
915 Vec distVec(4.0*getSphere().getRadius());
918 const double coreVol =
919 2.0*getSphere().getVolume(
920 centrePt -
Vec(distVec[ONE], distVec[TWO], distVec[ZERO]),
921 centrePt +
Vec(distVec[ONE], distVec[TWO], distVec[ZERO])
924 v[ONE] = 0.50*(sphVol - 2.0*s[TWO] - coreVol);
925 v[TWO] = 0.50*(sphVol - 2.0*s[ONE] - coreVol);
926 v[ZERO] = 0.25*(sphVol - 2.0*v[ONE] - 2.0*v[TWO] - coreVol);
928 if (pt[ONE] <= centrePt[ONE])
930 if (pt[TWO] <= centrePt[TWO])
932 vol = coreVol + v[ONE] + v[TWO] + v[ZERO];
936 vol = v[TWO] + v[ZERO];
941 if (pt[TWO] <= centrePt[TWO])
943 vol = v[ONE] + v[ZERO];
953 if (pt[ONE] <= centrePt[ONE])
955 if (pt[TWO] <= centrePt[TWO])
960 getSphere().getSegmentVolume(
Plane(getNegNormal(ONE), pt))
962 getSphere().getSegmentVolume(
Plane(getNegNormal(TWO), pt));
966 vol = getSphere().getSegmentVolume(
Plane(getNormal(TWO), pt));
971 if (pt[TWO] <= centrePt[TWO])
973 vol = getSphere().getSegmentVolume(
Plane(getNormal(ONE), pt));
985 template <
int tmplDim,
typename TmplVec>
991 const double sphVol = getSphere().getVolume();
992 const Vec ¢rePt = getSphere().getCentre();
995 if (pt[0] <= centrePt[0])
997 if (pt[1] <= centrePt[1])
1002 getSphere().getSegmentVolume(
Plane(getNegNormal(0), pt))
1004 getSphere().getSegmentVolume(
Plane(getNegNormal(1), pt));
1008 vol = getSphere().getSegmentVolume(
Plane(getNormal(1), pt));
1013 if (pt[1] <= centrePt[1])
1015 vol = getSphere().getSegmentVolume(
Plane(getNormal(0), pt));
1023 else if (tmplDim == 3)
1025 const Vec diag = (centrePt - pt)*2.0;
1026 const Vec oppCorner = pt + diag;
1029 componentMin(pt, oppCorner),
1030 componentMax(pt, oppCorner)
1035 for (
int i = 0; i < tmplDim; i++)
1037 s[i] = getSphere().getSegmentVolume(
Plane(getNormal(i), box.
getMaxPt()));
1038 e[i] = getTwoPlaneVolume(box.
getMaxPt(), i);
1040 double v[tmplDim+1];
1041 v[0] = s[0] - 2.0*e[1] - 2.0*e[2];
1042 v[1] = s[1] - 2.0*e[0] - 2.0*e[2];
1043 v[2] = s[2] - 2.0*e[0] - 2.0*e[1];
1044 v[3] = sphVol - (4.0*e[0] + 4.0*e[1] + 4.0*e[2] + 2.0*v[0] + 2.0*v[1] + 2.0*v[2]);
1046 if (pt[0] <= centrePt[0])
1048 if (pt[1] <= centrePt[1])
1050 if (pt[2] <= centrePt[2])
1052 vol = v[0] + v[1] + v[2] + v[3] + e[0] + e[1] + e[2];
1056 vol = v[2] + e[0] + e[1];
1061 if (pt[2] <= centrePt[2])
1063 vol = v[1] + e[0] + e[2];
1073 if (pt[1] <= centrePt[1])
1075 if (pt[2] <= centrePt[2])
1077 vol = v[0] + e[1] + e[2];
1086 if (pt[2] <= centrePt[2])
1101 template <
int tmplDim,
typename TmplVec>
1107 if (getSphere().intersectsWith(vtx.
getPoint()))
1109 vol = getInsidePointVolume(vtx.
getPoint());
1113 vol = getOutsidePointVolume(vtx.
getPoint());
1118 template <
int tmplDim,
typename TmplVec>
1126 std::vector<double> vtxVol(getVertexBox().getNumVertices());
1127 for (
int i = 0; i < getVertexBox().getNumVertices(); i++)
1129 vtxVol[i] = getVolume(getVertexBox().getVertex(i));
1134 vol = vtxVol[0] - vtxVol[1] - vtxVol[2] + vtxVol[3];
1136 else if (tmplDim == 3)
1139 vtxVol[7] + vtxVol[6] + vtxVol[5] - vtxVol[4]
1141 vtxVol[3] - vtxVol[2] - vtxVol[1] + vtxVol[0];
1147 template <
int tmplDim,
typename TmplVec>
1152 for (
int i = 0; i < getVertexBox().getNumVertices(); i++)
1154 if (!sphere.
intersectsWith(getVertexBox().getVertex(i).getPoint()))
1162 template <
int tmplDim,
typename TmplVec>
1168 if (getBox().intersectsWith(sphere))
1170 if (sphereContainsBox(sphere))
1172 vol = getBox().getVolume();
1174 else if (getBox().contains(sphere))
1180 vol = getVertexVolume(sphere);