24 inline double square(
double val)
29 template <
int tmplDim,
typename TmplVec>
30 DimBasicBox<tmplDim,TmplVec>::DimBasicBox(
const Vec &minPt,
const Vec &maxPt)
36 template <
int tmplDim,
typename TmplVec>
37 const typename DimBasicBox<tmplDim,TmplVec>::Vec &DimBasicBox<tmplDim,TmplVec>::getMinPt()
const
42 template <
int tmplDim,
typename TmplVec>
43 const typename DimBasicBox<tmplDim,TmplVec>::Vec &DimBasicBox<tmplDim,TmplVec>::getMaxPt()
const
48 template <
int tmplDim,
typename TmplVec>
49 double DimBasicBox<tmplDim,TmplVec>::getVolume()
const
52 for (
int i = 0; i < tmplDim; i++)
54 product *= (this->getMaxPt()[i] - this->getMinPt()[i]);
59 template <
int tmplDim,
typename TmplVec>
60 template <
typename TmplSphere>
61 bool DimBasicBox<tmplDim,TmplVec>::intersectsWith(
const TmplSphere &sphere)
const
63 double distSqrd = 0.0;
64 for (
int i = 0; i < tmplDim; i++)
66 if (sphere.getCentre()[i] < this->getMinPt()[i])
68 distSqrd += square(sphere.getCentre()[i] - this->getMinPt()[i]);
70 else if (sphere.getCentre()[i] > this->getMaxPt()[i])
72 distSqrd += square(sphere.getCentre()[i] - this->getMaxPt()[i]);
75 return (distSqrd <= square(sphere.getRadius()));
78 template <
int tmplDim,
typename TmplVec>
79 bool DimBasicBox<tmplDim,TmplVec>::intersectsWith(
const Vec &pt)
const
81 for (
int i = 0; (i < tmplDim); i++)
83 if ((this->getMinPt()[i] > pt[i]) || (pt[i] > this->getMaxPt()[i]))
91 template <
int tmplDim,
typename TmplVec>
92 template <
typename TmplSphere>
93 bool DimBasicBox<tmplDim,TmplVec>::contains(
const TmplSphere &sphere)
const
95 for (
int i = 0; (i < tmplDim); i++)
98 pt[i] = sphere.getRadius();
101 !(intersectsWith(sphere.getCentre() + pt))
103 !(intersectsWith(sphere.getCentre() - pt))
112 template <
int tmplDim,
typename TmplVec>
113 double DimPlane<tmplDim,TmplVec>::norm(
const Vec &pt)
116 for (
int i = 0; i < tmplDim; i++)
123 template <
int tmplDim,
typename TmplVec>
124 double DimPlane<tmplDim,TmplVec>::dot(
const Vec &p1,
const Vec &p2)
127 for (
int i = 0; i < tmplDim; i++)
134 template <
int tmplDim,
typename TmplVec>
135 DimPlane<tmplDim,TmplVec>::DimPlane() : m_normal(), m_pt(), m_invNormalNorm(0.0)
139 template <
int tmplDim,
typename TmplVec>
140 DimPlane<tmplDim,TmplVec>::DimPlane(
const Vec &normal,
const Vec &pt)
143 m_invNormalNorm((1.0/norm(normal)))
147 template <
int tmplDim,
typename TmplVec>
148 DimPlane<tmplDim,TmplVec>::DimPlane(
const DimPlane &plane)
149 : m_normal(plane.m_normal),
151 m_invNormalNorm(plane.m_invNormalNorm)
155 template <
int tmplDim,
typename TmplVec>
156 DimPlane<tmplDim,TmplVec> &DimPlane<tmplDim,TmplVec>::operator=(
const DimPlane &plane)
158 m_normal = plane.m_normal;
160 m_invNormalNorm = plane.m_invNormalNorm;
164 template <
int tmplDim,
typename TmplVec>
165 double DimPlane<tmplDim,TmplVec>::getSignedDistanceTo(
const Vec &pt)
const
170 (dot(m_normal, pt) - dot(m_normal, m_pt))*m_invNormalNorm
174 template <
int tmplDim,
typename TmplVec>
175 double DimPlane<tmplDim,TmplVec>::getDistanceTo(
const Vec &pt)
const
177 return fabs(getSignedDistanceTo(pt));
180 template <
int tmplDim,
typename TmplVec>
181 const typename DimPlane<tmplDim,TmplVec>::Vec &DimPlane<tmplDim,TmplVec>::getNormal()
const
186 template <
int tmplDim,
typename TmplVec>
187 const double DimBasicSphere<tmplDim,TmplVec>::FOUR_THIRDS_PI = (4.0/3.0)*M_PI;
189 template <
int tmplDim,
typename TmplVec>
190 const double DimBasicSphere<tmplDim,TmplVec>::ONE_THIRD_PI = M_PI/3.0;
192 template <
int tmplDim,
typename TmplVec>
193 DimBasicSphere<tmplDim,TmplVec>::DimBasicSphere()
199 template <
int tmplDim,
typename TmplVec>
200 DimBasicSphere<tmplDim,TmplVec>::DimBasicSphere(
const Vec ¢rePt,
double radius)
201 : m_centre(centrePt),
206 template <
int tmplDim,
typename TmplVec>
207 DimBasicSphere<tmplDim,TmplVec>::DimBasicSphere(
const DimBasicSphere &sphere)
208 : m_centre(sphere.m_centre),
209 m_radius(sphere.m_radius)
213 template <
int tmplDim,
typename TmplVec>
214 DimBasicSphere<tmplDim,TmplVec> &DimBasicSphere<tmplDim,TmplVec>::operator=(
const DimBasicSphere &sphere)
216 m_centre = sphere.m_centre;
217 m_radius = sphere.m_radius;
221 template <
int tmplDim,
typename TmplVec>
222 double DimBasicSphere<tmplDim,TmplVec>::getRadius()
const
227 template <
int tmplDim,
typename TmplVec>
228 const typename DimBasicSphere<tmplDim,TmplVec>::Vec &
229 DimBasicSphere<tmplDim,TmplVec>::getCentre()
const
234 template <
int tmplDim,
typename TmplVec>
235 double DimBasicSphere<tmplDim,TmplVec>::getVolume()
const
237 return (tmplDim == 2) ? M_PI*getRadius()*getRadius() : FOUR_THIRDS_PI*getRadius()*getRadius()*getRadius();
240 inline void checkDomain(
double r,
double x1,
double y1,
double x2,
double y2)
242 const double rSqrd = r*r;
243 const double x1Sqrd = x1*x1;
244 const double x2Sqrd = x2*x2;
245 const double y1Sqrd = y1*y1;
246 const double y2Sqrd = y2*y2;
249 ((rSqrd - x1Sqrd - y1Sqrd) < 0)
251 ((rSqrd - x1Sqrd - y2Sqrd) < 0)
253 ((rSqrd - x2Sqrd - y1Sqrd) < 0)
255 ((rSqrd - x2Sqrd - y2Sqrd) < 0)
258 std::stringstream msg;
260 <<
"Invalid rectangular domain for sphere integration, points in domain "
261 <<
"(" << x1 <<
"," << y1 <<
"), (" << x2 <<
"," << y2 <<
" lie outside "
262 <<
"sphere of radius " << r <<
" centred at the origin.";
263 throw std::runtime_error(msg.str());
267 template <
int tmplDim,
typename TmplVec>
268 double DimBasicSphere<tmplDim,TmplVec>::getVolume(
const Vec &minPt,
const Vec &maxPt,
const int dimX,
const int dimY)
const
271 if ((tmplDim == 2) || (tmplDim == 3))
273 if (minPt[dimX] != maxPt[dimX])
275 const double x1 = minPt[dimX] - getCentre()[dimX];
276 const double x2 = maxPt[dimX] - getCentre()[dimX];
277 const double r = getRadius();
281 const double rSqrd = r*r;
282 const double x1Sqrd = x1*x1;
283 const double x2Sqrd = x2*x2;
289 x2*sqrt(rSqrd-x2Sqrd)
293 x1*sqrt(rSqrd-x1Sqrd)
296 else if (tmplDim == 3)
298 if (minPt[dimY] != maxPt[dimY])
300 const double y1 = minPt[dimY] - getCentre()[dimY];
301 const double y2 = maxPt[dimY] - getCentre()[dimY];
303 checkDomain(r, x1, y1, x2, y2);
317 const double t30 = y2*y2;
318 const double t31 = x2*x2;
319 const double t36 = r*r;
320 const double t59 = t31-t36;
321 const double t40 = sqrt(-t30-t59);
322 const double t10 = 1.0/t40;
323 const double t32 = x1*x1;
324 const double t54 = t32-t36;
325 const double t42 = sqrt(-t30-t54);
326 const double t14 = 1.0/t42;
327 const double t64 = -atan(t10*x2)+atan(t14*x1);
328 const double t27 = y1*y1;
329 const double t39 = sqrt(-t27-t59);
330 const double t9 = 1.0/t39;
331 const double t41 = sqrt(-t27-t54);
332 const double t12 = 1.0/t41;
333 const double t63 = atan(t12*x1)-atan(t9*x2);
334 const double t62 = -atan(y2*t14)+atan(y1*t12);
335 const double t61 = -atan(y1*t9)+atan(t10*y2);
336 const double t37 = sqrt(t31);
337 const double t21 = 1.0/t37;
338 const double t60 = t21*t9;
339 const double t38 = sqrt(t32);
340 const double t24 = 1.0/t38;
341 const double t58 = t24*t14;
342 const double t57 = t24*t12;
343 const double t56 = t37*t38;
344 const double t55 = t21*t10;
345 const double t53 = 2.0*x2;
346 const double t52 = 2.0*x1;
347 const double t51 = t42*t56;
348 const double t28 = t27*y1;
349 const double t50 = t28-t36*y1;
350 const double t34 = t30*y2;
351 const double t49 = t34-t36*y2;
352 const double t48 = t41*t51;
353 const double t35 = t36*r;
354 const double t33 = t31*x2;
355 const double t29 = t32*x1;
356 const double t26 = r*y2;
357 const double t25 = y1*r;
362 (-2.0*t33*y1-t50*t53)*t40*t48
365 (2.0*t33*y2+t49*t53)*t48
368 (2.0*t29*y1+t50*t52)*t51
371 (-2.0*t29*y2-t49*t52)*t56
386 -atan((-t26+t59)*t55)
395 (-t64*t34+t61*t33+t62*t29+t63*t28+3.0*(t64*y2-t63*y1-t61*x2-t62*x1)*t36)*t37
410 template <
int tmplDim,
typename TmplVec>
411 bool DimBasicSphere<tmplDim,TmplVec>::intersectsWith(
const Vec &pt)
const
413 double distSqrd = 0.0;
414 for (
int i = 0; i < tmplDim; i++)
416 distSqrd += square(getCentre()[i] - pt[i]);
418 return (distSqrd <= square(getRadius()));
421 template <
int tmplDim,
typename TmplVec>
422 double DimBasicSphere<tmplDim,TmplVec>::getSegmentVolume(
const Plane &plane)
const
425 if ((tmplDim == 2) || (tmplDim == 3))
427 const double signedD = plane.getSignedDistanceTo(getCentre());
428 const double d = fabs(signedD);
434 const double rSqrd = getRadius()*getRadius();
435 vol = rSqrd*acos(d/getRadius()) - d*sqrt(rSqrd - d*d);
437 else if (tmplDim == 3)
440 const double h = getRadius() - d;
441 vol = ONE_THIRD_PI*h*h*(3.0*getRadius()-h);
443 vol = ((signedD < 0) ? vol : getVolume() - vol);
449 template <
int tmplDim,
typename TmplVec>
450 typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec
451 IntersectionVolCalculator<tmplDim,TmplVec>::getNormal(
int dim)
458 template <
int tmplDim,
typename TmplVec>
459 typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec
460 IntersectionVolCalculator<tmplDim,TmplVec>::getNegNormal(
int dim)
467 template <
int tmplDim,
typename TmplVec>
468 IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::VolumeSphere()
474 template <
int tmplDim,
typename TmplVec>
475 IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::VolumeSphere(
476 const BasicSphere &sphere
479 m_volume(sphere.getVolume())
483 template <
int tmplDim,
typename TmplVec>
484 IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::VolumeSphere(
485 const VolumeSphere &sphere
487 : m_sphere(BasicSphere(sphere.getCentre(), sphere.getRadius())),
488 m_volume(sphere.m_volume)
492 template <
int tmplDim,
typename TmplVec>
493 typename IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere &
494 IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::operator=(
495 const VolumeSphere &sphere
498 m_sphere = sphere.m_sphere;
499 m_volume = sphere.m_volume;
503 template <
int tmplDim,
typename TmplVec>
504 double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getRadius()
const
506 return m_sphere.getRadius();
509 template <
int tmplDim,
typename TmplVec>
510 const typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec &
511 IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getCentre()
const
513 return m_sphere.getCentre();
516 template <
int tmplDim,
typename TmplVec>
517 double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getVolume()
const
522 template <
int tmplDim,
typename TmplVec>
523 double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getVolume(
530 return m_sphere.getVolume(minPt, maxPt, dimX, dimY);
533 template <
int tmplDim,
typename TmplVec>
534 double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::calcVolume()
const
536 return m_sphere.getVolume();
539 template <
int tmplDim,
typename TmplVec>
540 bool IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::intersectsWith(
const Vec &pt)
const
542 return m_sphere.intersectsWith(pt);
545 template <
int tmplDim,
typename TmplVec>
546 double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getSegmentVolume(
const Plane &plane)
const
548 return m_sphere.getSegmentVolume(plane);
551 template <
int tmplDim,
typename TmplVec>
552 IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::Vertex() : m_pt()
556 template <
int tmplDim,
typename TmplVec>
557 IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::Vertex(
const Vec &pt) : m_pt(pt)
561 template <
int tmplDim,
typename TmplVec>
562 IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::Vertex(
const Vertex &vtx) : m_pt(vtx.m_pt)
566 template <
int tmplDim,
typename TmplVec>
567 typename IntersectionVolCalculator<tmplDim,TmplVec>::Vertex &
568 IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::operator=(
const Vertex &vtx)
574 template <
int tmplDim,
typename TmplVec>
575 const typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec &
576 IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::getPoint()
const
581 template <
int tmplDim,
typename TmplVec>
582 void IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::setPoint(
const Vec &pt)
587 template <
int tmplDim,
typename TmplVec>
588 IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::VertexBox(
597 template <
int tmplDim,
typename TmplVec>
598 IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::VertexBox(
603 for (
int i = 0; i < getNumVertices(); i++)
605 m_vertexArray[i] = box.m_vertexArray[i];
609 template <
int tmplDim,
typename TmplVec>
610 typename IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox &
611 IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::operator=(
615 BasicBox::operator=(box);
616 for (
int i = 0; i < getNumVertices(); i++)
618 m_vertexArray[i] = box.m_vertexArray[i];
623 template <
int tmplDim,
typename TmplVec>
624 void IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::createVertices()
627 m_vertexArray[j].setPoint(this->getMinPt());
629 for (j++; i < tmplDim; i++, j++)
631 Vec pt = this->getMinPt();
632 pt[i] = this->getMaxPt()[i];
633 m_vertexArray[j].setPoint(pt);
636 m_vertexArray[j].setPoint(this->getMaxPt());
637 for (i = 0, j++; i < tmplDim && j < s_numVertices; i++, j++)
639 Vec pt = this->getMaxPt();
640 pt[i] = this->getMinPt()[i];
641 m_vertexArray[j] = pt;
645 template <
int tmplDim,
typename TmplVec>
646 const typename IntersectionVolCalculator<tmplDim,TmplVec>::Vertex &
647 IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::getVertex(
int i)
const
649 return m_vertexArray[i];
652 template <
int tmplDim,
typename TmplVec>
653 int IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::getNumVertices()
655 return s_numVertices;
658 template <
int tmplDim,
typename TmplVec>
659 IntersectionVolCalculator<tmplDim,TmplVec>::IntersectionVolCalculator(
662 : m_sphere(BasicSphere(Vec(), 1.0)),
667 template <
int tmplDim,
typename TmplVec>
668 const typename IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere &
669 IntersectionVolCalculator<tmplDim,TmplVec>::getSphere()
const
674 template <
int tmplDim,
typename TmplVec>
675 void IntersectionVolCalculator<tmplDim,TmplVec>::setSphere(
676 const BasicSphere &sphere
682 template <
int tmplDim,
typename TmplVec>
683 const typename IntersectionVolCalculator<tmplDim,TmplVec>::BasicBox &
684 IntersectionVolCalculator<tmplDim,TmplVec>::getBox()
const
689 template <
int tmplDim,
typename TmplVec>
690 const typename IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox &
691 IntersectionVolCalculator<tmplDim,TmplVec>::getVertexBox()
const
696 template <
int tmplDim,
typename TmplVec>
697 typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec
698 IntersectionVolCalculator<tmplDim,TmplVec>::componentMin(
704 for (
int i = 0; i < tmplDim; i++)
706 m[i] = min(p1[i], p2[i]);
711 template <
int tmplDim,
typename TmplVec>
712 typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec
713 IntersectionVolCalculator<tmplDim,TmplVec>::componentMax(
719 for (
int i = 0; i < tmplDim; i++)
721 m[i] = max(p1[i], p2[i]);
726 template <
int tmplDim,
typename TmplVec>
727 double IntersectionVolCalculator<tmplDim,TmplVec>::getInsidePointVolume(
732 const Vec ¢rePt = getSphere().getCentre();
733 const Vec diag = (getSphere().getCentre() - pt)*2.0;
734 const Vec oppCorner = pt + diag;
737 componentMin(pt, oppCorner),
738 componentMax(pt, oppCorner)
740 const double boxVol = box.getVolume();
741 const double sphVol = getSphere().getVolume();
745 for (
int i = 0; i < tmplDim; i++)
747 s[i] = getSphere().getSegmentVolume(
Plane(getNormal((i+1)%tmplDim), box.getMaxPt()));
751 v[0] = 0.50*(sphVol - 2.0*s[0] - boxVol);
752 v[1] = 0.50*(sphVol - 2.0*s[1] - boxVol);
753 v[2] = 0.25*(sphVol - 2.0*v[0] - 2.0*v[1] - boxVol);
755 if (pt[0] <= centrePt[0])
757 if (pt[1] <= centrePt[1])
759 vol = boxVol + v[0] + v[1] + v[2];
768 if (pt[1] <= centrePt[1])
778 else if (tmplDim == 3)
782 2.0*getSphere().getVolume(
793 2.0*getSphere().getVolume(
804 2.0*getSphere().getVolume(
815 e[0] = 0.250*(sphVol - 2.0*s[1] - 2.0*v[0] - 2.0*v[1] - boxVol);
816 e[1] = 0.250*(sphVol - 2.0*s[2] - 2.0*v[1] - 2.0*v[2] - boxVol);
817 e[2] = 0.250*(sphVol - 2.0*s[0] - 2.0*v[0] - 2.0*v[2] - boxVol);
823 2.0*v[0] - 2.0*v[1] - 2.0*v[2]
825 4.0*e[0] - 4.0*e[1] - 4.0*e[2]
830 if (pt[0] <= centrePt[0])
832 if (pt[1] <= centrePt[1])
834 if (pt[2] <= centrePt[2])
836 vol = boxVol + v[0] + v[1] + v[2] + v[3] + e[0] + e[1] + e[2];
840 vol = v[2] + v[3] + e[1] + e[2];
845 if (pt[2] <= centrePt[2])
847 vol = v[1] + v[3] + e[0] + e[1];
857 if (pt[1] <= centrePt[1])
859 if (pt[2] <= centrePt[2])
861 vol = v[0] + v[3] + e[0] + e[2];
870 if (pt[2] <= centrePt[2])
885 template <
int tmplDim,
typename TmplVec>
886 double IntersectionVolCalculator<tmplDim,TmplVec>::getTwoPlaneVolume(
891 const int ZERO = orientDim;
892 const int ONE = (orientDim+1) % tmplDim;
893 const int TWO = (orientDim+2) % tmplDim;
896 const double sphVol = getSphere().getVolume();
897 const Vec ¢rePt = getSphere().getCentre();
899 if ((square(pt[ONE]-centrePt[ONE]) + square(pt[TWO]-centrePt[TWO])) < square(getSphere().getRadius()))
901 Plane plane[tmplDim];
902 plane[ZERO] =
Plane();
903 plane[ONE] =
Plane(getNormal(ONE), pt);
904 plane[TWO] =
Plane(getNormal(TWO), pt);
906 const double halfSphVol = sphVol*0.5;
909 s[ONE] = getSphere().getSegmentVolume(plane[ONE]);
910 s[TWO] = getSphere().getSegmentVolume(plane[TWO]);
911 s[ONE] = ((s[ONE] > halfSphVol) ? (sphVol - s[ONE]) : s[ONE]);
912 s[TWO] = ((s[TWO] > halfSphVol) ? (sphVol - s[TWO]) : s[TWO]);
914 Vec distVec(4.0*getSphere().getRadius());
915 distVec[ONE] = plane[ONE].getDistanceTo(centrePt);
916 distVec[TWO] = plane[TWO].getDistanceTo(centrePt);
917 const double coreVol =
918 2.0*getSphere().getVolume(
919 centrePt - Vec(distVec[ONE], distVec[TWO], distVec[ZERO]),
920 centrePt + Vec(distVec[ONE], distVec[TWO], distVec[ZERO])
923 v[ONE] = 0.50*(sphVol - 2.0*s[TWO] - coreVol);
924 v[TWO] = 0.50*(sphVol - 2.0*s[ONE] - coreVol);
925 v[ZERO] = 0.25*(sphVol - 2.0*v[ONE] - 2.0*v[TWO] - coreVol);
927 if (pt[ONE] <= centrePt[ONE])
929 if (pt[TWO] <= centrePt[TWO])
931 vol = coreVol + v[ONE] + v[TWO] + v[ZERO];
935 vol = v[TWO] + v[ZERO];
940 if (pt[TWO] <= centrePt[TWO])
942 vol = v[ONE] + v[ZERO];
952 if (pt[ONE] <= centrePt[ONE])
954 if (pt[TWO] <= centrePt[TWO])
959 getSphere().getSegmentVolume(
Plane(getNegNormal(ONE), pt))
961 getSphere().getSegmentVolume(
Plane(getNegNormal(TWO), pt));
965 vol = getSphere().getSegmentVolume(
Plane(getNormal(TWO), pt));
970 if (pt[TWO] <= centrePt[TWO])
972 vol = getSphere().getSegmentVolume(
Plane(getNormal(ONE), pt));
984 template <
int tmplDim,
typename TmplVec>
985 double IntersectionVolCalculator<tmplDim,TmplVec>::getOutsidePointVolume(
990 const double sphVol = getSphere().getVolume();
991 const Vec ¢rePt = getSphere().getCentre();
994 if (pt[0] <= centrePt[0])
996 if (pt[1] <= centrePt[1])
1001 getSphere().getSegmentVolume(
Plane(getNegNormal(0), pt))
1003 getSphere().getSegmentVolume(
Plane(getNegNormal(1), pt));
1007 vol = getSphere().getSegmentVolume(
Plane(getNormal(1), pt));
1012 if (pt[1] <= centrePt[1])
1014 vol = getSphere().getSegmentVolume(
Plane(getNormal(0), pt));
1022 else if (tmplDim == 3)
1024 const Vec diag = (centrePt - pt)*2.0;
1025 const Vec oppCorner = pt + diag;
1028 componentMin(pt, oppCorner),
1029 componentMax(pt, oppCorner)
1034 for (
int i = 0; i < tmplDim; i++)
1036 s[i] = getSphere().getSegmentVolume(
Plane(getNormal(i), box.getMaxPt()));
1037 e[i] = getTwoPlaneVolume(box.getMaxPt(), i);
1039 double v[tmplDim+1];
1040 v[0] = s[0] - 2.0*e[1] - 2.0*e[2];
1041 v[1] = s[1] - 2.0*e[0] - 2.0*e[2];
1042 v[2] = s[2] - 2.0*e[0] - 2.0*e[1];
1043 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]);
1045 if (pt[0] <= centrePt[0])
1047 if (pt[1] <= centrePt[1])
1049 if (pt[2] <= centrePt[2])
1051 vol = v[0] + v[1] + v[2] + v[3] + e[0] + e[1] + e[2];
1055 vol = v[2] + e[0] + e[1];
1060 if (pt[2] <= centrePt[2])
1062 vol = v[1] + e[0] + e[2];
1072 if (pt[1] <= centrePt[1])
1074 if (pt[2] <= centrePt[2])
1076 vol = v[0] + e[1] + e[2];
1085 if (pt[2] <= centrePt[2])
1100 template <
int tmplDim,
typename TmplVec>
1101 double IntersectionVolCalculator<tmplDim,TmplVec>::getVolume(
1106 if (getSphere().intersectsWith(vtx.getPoint()))
1108 vol = getInsidePointVolume(vtx.getPoint());
1112 vol = getOutsidePointVolume(vtx.getPoint());
1117 template <
int tmplDim,
typename TmplVec>
1118 double IntersectionVolCalculator<tmplDim,TmplVec>::getVertexVolume(
1119 const BasicSphere &sphere
1122 m_sphere = VolumeSphere(sphere);
1125 double vtxVol[getVertexBox().getNumVertices()];
1126 for (
int i = 0; i < getVertexBox().getNumVertices(); i++)
1128 vtxVol[i] = getVolume(getVertexBox().getVertex(i));
1133 vol = vtxVol[0] - vtxVol[1] - vtxVol[2] + vtxVol[3];
1135 else if (tmplDim == 3)
1138 vtxVol[7] + vtxVol[6] + vtxVol[5] - vtxVol[4]
1140 vtxVol[3] - vtxVol[2] - vtxVol[1] + vtxVol[0];
1146 template <
int tmplDim,
typename TmplVec>
1147 bool IntersectionVolCalculator<tmplDim,TmplVec>::sphereContainsBox(
1148 const BasicSphere &sphere
1151 for (
int i = 0; i < getVertexBox().getNumVertices(); i++)
1153 if (!sphere.intersectsWith(getVertexBox().getVertex(i).getPoint()))
1161 template <
int tmplDim,
typename TmplVec>
1162 double IntersectionVolCalculator<tmplDim,TmplVec>::getVolume(
1163 const BasicSphere &sphere
1167 if (getBox().intersectsWith(sphere))
1169 if (sphereContainsBox(sphere))
1171 vol = getBox().getVolume();
1173 else if (getBox().contains(sphere))
1175 vol = sphere.getVolume();
1179 vol = getVertexVolume(sphere);