ESyS-Particle  4.0.1
IntersectionVolCalculator.hpp
00001 
00002 //                                                         //
00003 // Copyright (c) 2003-2011 by The University of Queensland //
00004 // Earth Systems Science Computational Centre (ESSCC)      //
00005 // http://www.uq.edu.au/esscc                              //
00006 //                                                         //
00007 // Primary Business: Brisbane, Queensland, Australia       //
00008 // Licensed under the Open Software License version 3.0    //
00009 // http://www.opensource.org/licenses/osl-3.0.php          //
00010 //                                                         //
00012 
00013 
00014 
00015 #include <stdexcept>
00016 #include <sstream>
00017 
00018 namespace esys
00019 {
00020   namespace lsm
00021   {
00022     namespace impl
00023     {
00024       inline double square(double val)
00025       {
00026         return val*val;
00027       }
00028 
00029       template <int tmplDim, typename TmplVec>
00030       DimBasicBox<tmplDim,TmplVec>::DimBasicBox(const Vec &minPt, const Vec &maxPt)
00031         : m_minPt(minPt),
00032           m_maxPt(maxPt)
00033       {
00034       }
00035 
00036       template <int tmplDim, typename TmplVec>
00037       const typename DimBasicBox<tmplDim,TmplVec>::Vec &DimBasicBox<tmplDim,TmplVec>::getMinPt() const
00038       {
00039         return m_minPt;
00040       }
00041 
00042       template <int tmplDim, typename TmplVec>
00043       const typename DimBasicBox<tmplDim,TmplVec>::Vec &DimBasicBox<tmplDim,TmplVec>::getMaxPt() const
00044       {
00045         return m_maxPt;
00046       }
00047 
00048       template <int tmplDim, typename TmplVec>
00049       double DimBasicBox<tmplDim,TmplVec>::getVolume() const
00050       {
00051         double product = 1.0;
00052         for (int i = 0; i < tmplDim; i++)
00053         {
00054           product *= (this->getMaxPt()[i] - this->getMinPt()[i]);
00055         }
00056         return product;
00057       }
00058       
00059       template <int tmplDim, typename TmplVec>
00060       template <typename TmplSphere>
00061       bool DimBasicBox<tmplDim,TmplVec>::intersectsWith(const TmplSphere &sphere) const
00062       {
00063         double distSqrd = 0.0;
00064         for (int i = 0; i < tmplDim; i++)
00065         {
00066           if (sphere.getCentre()[i] < this->getMinPt()[i])
00067           {
00068             distSqrd += square(sphere.getCentre()[i] - this->getMinPt()[i]);
00069           }
00070           else if (sphere.getCentre()[i] > this->getMaxPt()[i])
00071           {
00072             distSqrd += square(sphere.getCentre()[i] - this->getMaxPt()[i]);
00073           }
00074         }
00075         return (distSqrd <= square(sphere.getRadius()));
00076       }
00077 
00078       template <int tmplDim, typename TmplVec>
00079       bool DimBasicBox<tmplDim,TmplVec>::intersectsWith(const Vec &pt) const
00080       {
00081         for (int i = 0; (i < tmplDim); i++)
00082         {
00083           if ((this->getMinPt()[i] > pt[i]) || (pt[i] > this->getMaxPt()[i]))
00084           {
00085             return false;
00086           }
00087         }
00088         return true;
00089       }
00090 
00091       template <int tmplDim, typename TmplVec>
00092       template <typename TmplSphere>
00093       bool DimBasicBox<tmplDim,TmplVec>::contains(const TmplSphere &sphere) const
00094       {
00095         for (int i = 0; (i < tmplDim); i++)
00096         {
00097           Vec pt = Vec(0.0);
00098           pt[i] = sphere.getRadius();
00099           if
00100           (
00101             !(intersectsWith(sphere.getCentre() + pt))
00102             ||
00103             !(intersectsWith(sphere.getCentre() - pt))
00104           )
00105           {
00106             return false;
00107           }
00108         }
00109         return true;
00110       }
00111 
00112       template <int tmplDim, typename TmplVec>
00113       double DimPlane<tmplDim,TmplVec>::norm(const Vec &pt)
00114       {
00115         double sum = 0.0;
00116         for (int i = 0; i < tmplDim; i++)
00117         {
00118           sum += pt[i]*pt[i];
00119         }
00120         return sqrt(sum);
00121       }
00122 
00123       template <int tmplDim, typename TmplVec>
00124       double DimPlane<tmplDim,TmplVec>::dot(const Vec &p1, const Vec &p2)
00125       {
00126         double sum = 0.0;
00127         for (int i = 0; i < tmplDim; i++)
00128         {
00129           sum += p1[i]*p2[i];
00130         }
00131         return sum;
00132       }
00133       
00134       template <int tmplDim, typename TmplVec>
00135       DimPlane<tmplDim,TmplVec>::DimPlane() : m_normal(), m_pt(), m_invNormalNorm(0.0)
00136       {
00137       }
00138 
00139       template <int tmplDim, typename TmplVec>
00140       DimPlane<tmplDim,TmplVec>::DimPlane(const Vec &normal, const Vec &pt)
00141         : m_normal(normal),
00142           m_pt(pt),
00143           m_invNormalNorm((1.0/norm(normal)))
00144       {
00145       }
00146 
00147       template <int tmplDim, typename TmplVec>
00148       DimPlane<tmplDim,TmplVec>::DimPlane(const DimPlane &plane)
00149         : m_normal(plane.m_normal),
00150           m_pt(plane.m_pt),
00151           m_invNormalNorm(plane.m_invNormalNorm)
00152       {
00153       }
00154 
00155       template <int tmplDim, typename TmplVec>
00156       DimPlane<tmplDim,TmplVec> &DimPlane<tmplDim,TmplVec>::operator=(const DimPlane &plane)
00157       {
00158         m_normal        = plane.m_normal;
00159         m_pt            = plane.m_pt;
00160         m_invNormalNorm = plane.m_invNormalNorm;
00161         return *this;
00162       }
00163 
00164       template <int tmplDim, typename TmplVec>
00165       double DimPlane<tmplDim,TmplVec>::getSignedDistanceTo(const Vec &pt) const
00166       {
00167         // http://mathworld.wolfram.com/Point-PlaneDistance.html
00168         return
00169           (
00170             (dot(m_normal, pt) - dot(m_normal, m_pt))*m_invNormalNorm
00171           );
00172       }
00173 
00174       template <int tmplDim, typename TmplVec>
00175       double DimPlane<tmplDim,TmplVec>::getDistanceTo(const Vec &pt) const
00176       {
00177         return fabs(getSignedDistanceTo(pt));
00178       }
00179 
00180       template <int tmplDim, typename TmplVec>
00181       const typename DimPlane<tmplDim,TmplVec>::Vec &DimPlane<tmplDim,TmplVec>::getNormal() const
00182       {
00183         return m_normal;
00184       }
00185 
00186       template <int tmplDim, typename TmplVec>
00187       const double DimBasicSphere<tmplDim,TmplVec>::FOUR_THIRDS_PI = (4.0/3.0)*M_PI;
00188 
00189       template <int tmplDim, typename TmplVec>
00190       const double DimBasicSphere<tmplDim,TmplVec>::ONE_THIRD_PI   = M_PI/3.0;
00191 
00192       template <int tmplDim, typename TmplVec>
00193       DimBasicSphere<tmplDim,TmplVec>::DimBasicSphere()
00194         : m_centre(),
00195           m_radius(0.0)
00196       {
00197       }
00198 
00199       template <int tmplDim, typename TmplVec>
00200       DimBasicSphere<tmplDim,TmplVec>::DimBasicSphere(const Vec &centrePt, double radius)
00201         : m_centre(centrePt),
00202           m_radius(radius)
00203       {
00204       }
00205 
00206       template <int tmplDim, typename TmplVec>
00207       DimBasicSphere<tmplDim,TmplVec>::DimBasicSphere(const DimBasicSphere &sphere)
00208         : m_centre(sphere.m_centre),
00209           m_radius(sphere.m_radius)
00210       {
00211       }
00212 
00213       template <int tmplDim, typename TmplVec>
00214       DimBasicSphere<tmplDim,TmplVec> &DimBasicSphere<tmplDim,TmplVec>::operator=(const DimBasicSphere &sphere)
00215       {
00216         m_centre = sphere.m_centre;
00217         m_radius = sphere.m_radius;
00218         return *this;
00219       }
00220 
00221       template <int tmplDim, typename TmplVec>
00222       double DimBasicSphere<tmplDim,TmplVec>::getRadius() const
00223       {
00224         return m_radius;
00225       }
00226 
00227       template <int tmplDim, typename TmplVec>
00228       const typename DimBasicSphere<tmplDim,TmplVec>::Vec &
00229       DimBasicSphere<tmplDim,TmplVec>::getCentre() const
00230       {
00231         return m_centre;
00232       }
00233 
00234       template <int tmplDim, typename TmplVec>
00235       double DimBasicSphere<tmplDim,TmplVec>::getVolume() const
00236       {
00237         return (tmplDim == 2) ? M_PI*getRadius()*getRadius() : FOUR_THIRDS_PI*getRadius()*getRadius()*getRadius();
00238       }
00239 
00240       inline void checkDomain(double r, double x1, double y1, double x2, double y2)
00241       {
00242         const double rSqrd = r*r;
00243         const double x1Sqrd = x1*x1;
00244         const double x2Sqrd = x2*x2;
00245         const double y1Sqrd = y1*y1;
00246         const double y2Sqrd = y2*y2;
00247         if
00248         (
00249           ((rSqrd - x1Sqrd - y1Sqrd) < 0)
00250           ||
00251           ((rSqrd - x1Sqrd - y2Sqrd) < 0)
00252           ||
00253           ((rSqrd - x2Sqrd - y1Sqrd) < 0)
00254           ||
00255           ((rSqrd - x2Sqrd - y2Sqrd) < 0)
00256         )
00257         {
00258           std::stringstream msg;
00259           msg 
00260             << "Invalid rectangular domain for sphere integration, points in domain "
00261             << "(" << x1 << "," << y1 << "), (" << x2 << "," << y2 << " lie outside "
00262             << "sphere of radius " << r << " centred at the origin.";
00263           throw std::runtime_error(msg.str());
00264         }
00265       }
00266 
00267       template <int tmplDim, typename TmplVec>
00268       double DimBasicSphere<tmplDim,TmplVec>::getVolume(const Vec &minPt, const Vec &maxPt, const int dimX, const int dimY) const
00269       {
00270         double vol = 0.0;
00271         if ((tmplDim == 2) || (tmplDim == 3))
00272         {
00273           if (minPt[dimX] != maxPt[dimX])
00274           {
00275             const double x1 = minPt[dimX] - getCentre()[dimX];
00276             const double x2 = maxPt[dimX] - getCentre()[dimX];
00277             const double r  = getRadius();
00278   
00279             if (tmplDim == 2)
00280             {
00281               const double rSqrd  = r*r;
00282               const double x1Sqrd = x1*x1;
00283               const double x2Sqrd = x2*x2;
00284 
00285               vol = 
00286                 (
00287                   rSqrd*asin(x2/r)
00288                   +
00289                   x2*sqrt(rSqrd-x2Sqrd)
00290                   -
00291                   rSqrd*asin(x1/r)
00292                   -
00293                   x1*sqrt(rSqrd-x1Sqrd)
00294                 )*0.5;
00295             }
00296             else if (tmplDim == 3)
00297             {
00298               if (minPt[dimY] != maxPt[dimY])
00299               {
00300                 const double y1 = minPt[dimY] - getCentre()[dimY];
00301                 const double y2 = maxPt[dimY] - getCentre()[dimY];
00302 
00303                 checkDomain(r, x1, y1, x2, y2);
00304 
00305                 //
00306                 // Matlab6/maple generated code:
00307                 //
00308                 // syms x y x1 x2 y1 y2 r real
00309                 // sphereIntegral = int(int('sqrt(r^2-x^2-y^2)', x, x1, x2), y, y1, y2)
00310                 // maple('readlib(codegen)')
00311                 // maple('readlib(optimize)')
00312                 // optimizedIntegral = maple('optimize', sphereIntegral, 'tryhard')
00313                 // maple('cost', optimizedIntegral)
00314                 //
00315                 //43*additions+82*multiplications+6*divisions+22*functions+42*assignments
00316 
00317                 const double t30 = y2*y2; //y2^2;
00318                 const double t31 = x2*x2; //x2^2;
00319                 const double t36 = r*r;   //r^2;
00320                 const double t59 = t31-t36;
00321                 const double t40 = sqrt(-t30-t59); //pow(-t30-t59,1/2);
00322                 const double t10 = 1.0/t40;
00323                 const double t32 = x1*x1; //x1^2;
00324                 const double t54 = t32-t36;
00325                 const double t42 = sqrt(-t30-t54); //pow(-t30-t54,1/2);
00326                 const double t14 = 1.0/t42;
00327                 const double t64 = -atan(t10*x2)+atan(t14*x1);
00328                 const double t27 = y1*y1; //y1^2;
00329                 const double t39 = sqrt(-t27-t59); //pow(-t27-t59,1/2);
00330                 const double t9  = 1.0/t39;
00331                 const double t41 = sqrt(-t27-t54); //pow(-t27-t54,1/2);
00332                 const double t12 = 1.0/t41;
00333                 const double t63 = atan(t12*x1)-atan(t9*x2);
00334                 const double t62 = -atan(y2*t14)+atan(y1*t12);
00335                 const double t61 = -atan(y1*t9)+atan(t10*y2);
00336                 const double t37 = sqrt(t31); //pow(t31,1/2);
00337                 const double t21 = 1.0/t37;
00338                 const double t60 = t21*t9;
00339                 const double t38 = sqrt(t32); //pow(t32,1/2);
00340                 const double t24 = 1.0/t38;
00341                 const double t58 = t24*t14;
00342                 const double t57 = t24*t12;
00343                 const double t56 = t37*t38;
00344                 const double t55 = t21*t10;
00345                 const double t53 = 2.0*x2;
00346                 const double t52 = 2.0*x1;
00347                 const double t51 = t42*t56;
00348                 const double t28 = t27*y1;
00349                 const double t50 = t28-t36*y1;
00350                 const double t34 = t30*y2;
00351                 const double t49 = t34-t36*y2;
00352                 const double t48 = t41*t51;
00353                 const double t35 = t36*r;
00354                 const double t33 = t31*x2;
00355                 const double t29 = t32*x1;
00356                 const double t26 = r*y2;
00357                 const double t25 = y1*r;
00358                 vol = 
00359                   (-1.0/6.0)
00360                   *
00361                   (
00362                     (-2.0*t33*y1-t50*t53)*t40*t48
00363                     +
00364                     (
00365                       (2.0*t33*y2+t49*t53)*t48
00366                       +
00367                       (
00368                         (2.0*t29*y1+t50*t52)*t51
00369                         +
00370                         (
00371                           (-2.0*t29*y2-t49*t52)*t56
00372                           +
00373                           (
00374                             (
00375                               atan((t25+t54)*t57)
00376                               +
00377                               atan((-t26+t54)*t58)
00378                               -
00379                               atan((t26+t54)*t58)
00380                               -
00381                               atan((-t25+t54)*t57)
00382                             )*t35*t37*x1
00383                             +
00384                             (
00385                               (
00386                                 -atan((-t26+t59)*t55)
00387                                 -
00388                                 atan((t25+t59)*t60)
00389                                 +
00390                                 atan((-t25+t59)*t60)
00391                                 +
00392                                 atan((t26+t59)*t55)
00393                               )*t35*x2
00394                               +
00395                               (-t64*t34+t61*t33+t62*t29+t63*t28+3.0*(t64*y2-t63*y1-t61*x2-t62*x1)*t36)*t37
00396                             )*t38
00397                           )*t42
00398                         )*t41
00399                       )*t40
00400                     )*t39
00401                   )*t14*t9*t55*t57;
00402               }
00403             }
00404           }
00405         }
00406         
00407         return vol;
00408       }
00409 
00410       template <int tmplDim, typename TmplVec>
00411       bool DimBasicSphere<tmplDim,TmplVec>::intersectsWith(const Vec &pt) const
00412       {
00413         double distSqrd = 0.0;
00414         for (int i = 0; i < tmplDim; i++)
00415         {
00416           distSqrd += square(getCentre()[i] - pt[i]);
00417         }
00418         return (distSqrd <= square(getRadius()));
00419       }
00420 
00421       template <int tmplDim, typename TmplVec>
00422       double DimBasicSphere<tmplDim,TmplVec>::getSegmentVolume(const Plane &plane) const
00423       {
00424         double vol = 0.0;
00425         if ((tmplDim == 2) || (tmplDim == 3))
00426         {
00427           const double signedD = plane.getSignedDistanceTo(getCentre());
00428           const double d = fabs(signedD);
00429           if (d < getRadius())
00430           {
00431             if (tmplDim == 2)
00432             {
00433               // http://mathworld.wolfram.com/CircularSegment.html
00434               const double rSqrd = getRadius()*getRadius();
00435               vol = rSqrd*acos(d/getRadius()) - d*sqrt(rSqrd - d*d);
00436             }
00437             else if (tmplDim == 3)
00438             {
00439               // http://mathworld.wolfram.com/SphericalCap.html
00440               const double h = getRadius() - d;
00441               vol = ONE_THIRD_PI*h*h*(3.0*getRadius()-h);
00442             }
00443             vol = ((signedD < 0) ? vol : getVolume() - vol);
00444           }
00445         }
00446         return vol;
00447       }
00448 
00449       template <int tmplDim, typename TmplVec>
00450       typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec
00451       IntersectionVolCalculator<tmplDim,TmplVec>::getNormal(int dim)
00452       {
00453         Vec n = Vec(0.0);
00454         n[dim] = 1.0;
00455         return n;
00456       }
00457 
00458       template <int tmplDim, typename TmplVec>
00459       typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec
00460       IntersectionVolCalculator<tmplDim,TmplVec>::getNegNormal(int dim)
00461       {
00462         Vec n = Vec(0.0);
00463         n[dim] = -1.0;
00464         return n;
00465       }
00466 
00467       template <int tmplDim, typename TmplVec>
00468       IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::VolumeSphere()
00469         : m_sphere(),
00470           m_volume(0.0)
00471       {
00472       }
00473 
00474       template <int tmplDim, typename TmplVec>
00475       IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::VolumeSphere(
00476         const BasicSphere &sphere
00477       )
00478         : m_sphere(sphere),
00479           m_volume(sphere.getVolume())
00480       {
00481       }
00482 
00483       template <int tmplDim, typename TmplVec>
00484       IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::VolumeSphere(
00485         const VolumeSphere &sphere
00486       )
00487         : m_sphere(BasicSphere(sphere.getCentre(), sphere.getRadius())),
00488           m_volume(sphere.m_volume)
00489       {
00490       }
00491 
00492       template <int tmplDim, typename TmplVec>
00493       typename IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere &
00494       IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::operator=(
00495         const VolumeSphere &sphere
00496       )
00497       {
00498         m_sphere = sphere.m_sphere;
00499         m_volume = sphere.m_volume;
00500         return *this;
00501       }
00502 
00503       template <int tmplDim, typename TmplVec>
00504       double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getRadius() const
00505       {
00506         return m_sphere.getRadius();
00507       }
00508 
00509       template <int tmplDim, typename TmplVec>
00510       const typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec &
00511       IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getCentre() const
00512       {
00513         return m_sphere.getCentre();
00514       }
00515 
00516       template <int tmplDim, typename TmplVec>
00517       double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getVolume() const
00518       {
00519         return m_volume;
00520       }
00521 
00522       template <int tmplDim, typename TmplVec>
00523       double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getVolume(
00524         const Vec &minPt,
00525         const Vec &maxPt,
00526         const int dimX,
00527         const int dimY
00528       ) const
00529       {
00530         return m_sphere.getVolume(minPt, maxPt, dimX, dimY);
00531       }
00532 
00533       template <int tmplDim, typename TmplVec>
00534       double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::calcVolume() const
00535       {
00536         return m_sphere.getVolume();
00537       }
00538 
00539       template <int tmplDim, typename TmplVec>
00540       bool IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::intersectsWith(const Vec &pt) const
00541       {
00542         return m_sphere.intersectsWith(pt);
00543       }
00544 
00545       template <int tmplDim, typename TmplVec>
00546       double IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere::getSegmentVolume(const Plane &plane) const
00547       {
00548         return m_sphere.getSegmentVolume(plane);
00549       }
00550 
00551       template <int tmplDim, typename TmplVec>
00552       IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::Vertex() : m_pt()
00553       {
00554       }
00555 
00556       template <int tmplDim, typename TmplVec>
00557       IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::Vertex(const Vec &pt) : m_pt(pt)
00558       {
00559       }
00560 
00561       template <int tmplDim, typename TmplVec>
00562       IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::Vertex(const Vertex &vtx) : m_pt(vtx.m_pt)
00563       {
00564       }
00565 
00566       template <int tmplDim, typename TmplVec>
00567       typename IntersectionVolCalculator<tmplDim,TmplVec>::Vertex &
00568       IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::operator=(const Vertex &vtx)
00569       {
00570         m_pt = vtx.m_pt;
00571         return *this;
00572       }
00573 
00574       template <int tmplDim, typename TmplVec>
00575       const typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec &
00576       IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::getPoint() const
00577       {
00578         return m_pt;
00579       }
00580 
00581       template <int tmplDim, typename TmplVec>
00582       void IntersectionVolCalculator<tmplDim,TmplVec>::Vertex::setPoint(const Vec &pt)
00583       {
00584         m_pt = pt;
00585       }
00586 
00587       template <int tmplDim, typename TmplVec>
00588       IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::VertexBox(
00589         const BasicBox &box
00590       )
00591         : BasicBox(box),
00592           m_vertexArray()
00593       {
00594         createVertices();
00595       }
00596 
00597       template <int tmplDim, typename TmplVec>
00598       IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::VertexBox(
00599         const VertexBox &box
00600       )
00601         : BasicBox(box)
00602       {
00603         for (int i = 0; i < getNumVertices(); i++)
00604         {
00605           m_vertexArray[i] = box.m_vertexArray[i];
00606         }
00607       }
00608 
00609       template <int tmplDim, typename TmplVec>
00610       typename IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox &
00611       IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::operator=(
00612         const VertexBox &box
00613       )
00614       {
00615         BasicBox::operator=(box);
00616         for (int i = 0; i < getNumVertices(); i++)
00617         {
00618           m_vertexArray[i] = box.m_vertexArray[i];
00619         }
00620         return *this;
00621       }
00622 
00623       template <int tmplDim, typename TmplVec>
00624       void IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::createVertices()
00625       {
00626         int j = 0;
00627         m_vertexArray[j].setPoint(this->getMinPt());
00628         int i = 0;
00629         for (j++; i < tmplDim; i++, j++)
00630         {
00631           Vec pt = this->getMinPt();
00632           pt[i] = this->getMaxPt()[i];
00633           m_vertexArray[j].setPoint(pt);
00634         }
00635 
00636         m_vertexArray[j].setPoint(this->getMaxPt());
00637         for (i = 0, j++; i < tmplDim && j < s_numVertices; i++, j++)
00638         {
00639           Vec pt = this->getMaxPt();
00640           pt[i] = this->getMinPt()[i];
00641           m_vertexArray[j] = pt;
00642         }
00643       }
00644 
00645       template <int tmplDim, typename TmplVec>
00646       const typename IntersectionVolCalculator<tmplDim,TmplVec>::Vertex &
00647       IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::getVertex(int i) const
00648       {
00649         return m_vertexArray[i];
00650       }
00651 
00652       template <int tmplDim, typename TmplVec>
00653       int IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox::getNumVertices()
00654       {
00655         return s_numVertices;
00656       }
00657 
00658       template <int tmplDim, typename TmplVec>
00659       IntersectionVolCalculator<tmplDim,TmplVec>::IntersectionVolCalculator(
00660         const BasicBox &box
00661       )
00662         : m_sphere(BasicSphere(Vec(), 1.0)),
00663           m_box(box)
00664       {
00665       }
00666 
00667       template <int tmplDim, typename TmplVec>
00668       const typename IntersectionVolCalculator<tmplDim,TmplVec>::VolumeSphere &
00669       IntersectionVolCalculator<tmplDim,TmplVec>::getSphere() const
00670       {
00671         return m_sphere;
00672       }
00673 
00674       template <int tmplDim, typename TmplVec>
00675       void IntersectionVolCalculator<tmplDim,TmplVec>::setSphere(
00676         const BasicSphere &sphere
00677       )
00678       {
00679         m_sphere = sphere;
00680       }
00681 
00682       template <int tmplDim, typename TmplVec>
00683       const typename IntersectionVolCalculator<tmplDim,TmplVec>::BasicBox &
00684       IntersectionVolCalculator<tmplDim,TmplVec>::getBox() const
00685       {
00686         return m_box;
00687       }
00688 
00689       template <int tmplDim, typename TmplVec>
00690       const typename IntersectionVolCalculator<tmplDim,TmplVec>::VertexBox &
00691       IntersectionVolCalculator<tmplDim,TmplVec>::getVertexBox() const
00692       {
00693         return m_box;
00694       }
00695 
00696       template <int tmplDim, typename TmplVec>
00697       typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec
00698       IntersectionVolCalculator<tmplDim,TmplVec>::componentMin(
00699         const Vec &p1,
00700         const Vec &p2
00701       )
00702       {
00703         Vec m;
00704         for (int i = 0; i < tmplDim; i++)
00705         {
00706           m[i] = min(p1[i], p2[i]);
00707         }
00708         return m;
00709       }
00710 
00711       template <int tmplDim, typename TmplVec>
00712       typename IntersectionVolCalculator<tmplDim,TmplVec>::Vec
00713       IntersectionVolCalculator<tmplDim,TmplVec>::componentMax(
00714         const Vec &p1,
00715         const Vec &p2
00716       )
00717       {
00718         Vec m;
00719         for (int i = 0; i < tmplDim; i++)
00720         {
00721           m[i] = max(p1[i], p2[i]);
00722         }
00723         return m;
00724       }
00725 
00726       template <int tmplDim, typename TmplVec>
00727       double IntersectionVolCalculator<tmplDim,TmplVec>::getInsidePointVolume(
00728         const Vec &pt
00729       ) const
00730       {
00731         double vol = 0.0;
00732         const Vec &centrePt = getSphere().getCentre();
00733         const Vec diag = (getSphere().getCentre() - pt)*2.0;
00734         const Vec oppCorner = pt + diag;
00735         BasicBox box =
00736           BasicBox(
00737             componentMin(pt, oppCorner),
00738             componentMax(pt, oppCorner)
00739           );
00740         const double boxVol  = box.getVolume();
00741         const double sphVol  = getSphere().getVolume();
00742 
00743         double s[tmplDim];
00744         double v[tmplDim+1];
00745         for (int i = 0; i < tmplDim; i++)
00746         {
00747           s[i] = getSphere().getSegmentVolume(Plane(getNormal((i+1)%tmplDim), box.getMaxPt()));
00748         }
00749         if (tmplDim == 2)
00750         {
00751           v[0] = 0.50*(sphVol - 2.0*s[0] - boxVol);
00752           v[1] = 0.50*(sphVol - 2.0*s[1] - boxVol);
00753           v[2] = 0.25*(sphVol - 2.0*v[0] - 2.0*v[1] - boxVol);
00754 
00755           if (pt[0] <= centrePt[0])
00756           {
00757             if (pt[1] <= centrePt[1])
00758             {
00759               vol = boxVol + v[0] + v[1] + v[2];
00760             }
00761             else
00762             {
00763               vol = v[1] + v[2];
00764             }
00765           }
00766           else
00767           {
00768             if (pt[1] <= centrePt[1])
00769             {
00770               vol = v[0] + v[2];
00771             }
00772             else
00773             {
00774               vol = v[2];
00775             }
00776           }
00777         }
00778         else if (tmplDim == 3)
00779         {
00780           v[0] = 
00781             0.500*(
00782               2.0*getSphere().getVolume(
00783                 box.getMinPt(),
00784                 box.getMaxPt(),
00785                 1,
00786                 2
00787               )
00788               -
00789               boxVol
00790             );
00791           v[1] = 
00792             0.500*(
00793               2.0*getSphere().getVolume(
00794                 box.getMinPt(),
00795                 box.getMaxPt(),
00796                 0,
00797                 2
00798               )
00799               -
00800               boxVol
00801             );
00802           v[2] = 
00803             0.500*(
00804               2.0*getSphere().getVolume(
00805                 box.getMinPt(),
00806                 box.getMaxPt(),
00807                 0,
00808                 1
00809               )
00810               -
00811               boxVol
00812             );
00813 
00814           double e[3];
00815           e[0] = 0.250*(sphVol - 2.0*s[1] - 2.0*v[0] - 2.0*v[1] - boxVol);
00816           e[1] = 0.250*(sphVol - 2.0*s[2] - 2.0*v[1] - 2.0*v[2] - boxVol);
00817           e[2] = 0.250*(sphVol - 2.0*s[0] - 2.0*v[0] - 2.0*v[2] - boxVol);
00818           
00819           v[3] = 
00820             0.125*(
00821               sphVol
00822               -
00823               2.0*v[0] - 2.0*v[1] - 2.0*v[2]
00824               -
00825               4.0*e[0] - 4.0*e[1] - 4.0*e[2]
00826               -
00827               boxVol
00828             );
00829 
00830           if (pt[0] <= centrePt[0])
00831           {
00832             if (pt[1] <= centrePt[1])
00833             {
00834               if (pt[2] <= centrePt[2])
00835               {
00836                 vol = boxVol + v[0] + v[1] + v[2] + v[3] + e[0] + e[1] + e[2];
00837               }
00838               else
00839               {
00840                 vol = v[2] + v[3] + e[1] + e[2];
00841               }
00842             }
00843             else
00844             {
00845               if (pt[2] <= centrePt[2])
00846               {
00847                 vol = v[1] + v[3] + e[0] + e[1];
00848               }
00849               else
00850               {
00851                 vol = v[3] + e[1];
00852               }
00853             }
00854           }
00855           else
00856           {
00857             if (pt[1] <= centrePt[1])
00858             {
00859               if (pt[2] <= centrePt[2])
00860               {
00861                 vol = v[0] + v[3] + e[0] + e[2];
00862               }
00863               else
00864               {
00865                 vol = v[3] + e[2];
00866               }
00867             }
00868             else
00869             {
00870               if (pt[2] <= centrePt[2])
00871               {
00872                 vol = v[3] + e[0];
00873               }
00874               else
00875               {
00876                 vol = v[3];
00877               }
00878             }
00879           }
00880         }
00881 
00882         return vol;
00883       }
00884 
00885       template <int tmplDim, typename TmplVec>
00886       double IntersectionVolCalculator<tmplDim,TmplVec>::getTwoPlaneVolume(
00887         const Vec &pt,
00888         const int orientDim
00889       ) const
00890       {
00891         const int ZERO = orientDim;
00892         const int ONE  = (orientDim+1) % tmplDim;
00893         const int TWO  = (orientDim+2) % tmplDim;
00894 
00895         double vol = 0.0;
00896         const double sphVol = getSphere().getVolume();
00897         const Vec &centrePt = getSphere().getCentre();
00898 
00899         if ((square(pt[ONE]-centrePt[ONE]) + square(pt[TWO]-centrePt[TWO])) < square(getSphere().getRadius()))
00900         {
00901           Plane plane[tmplDim];
00902           plane[ZERO] = Plane();
00903           plane[ONE]  = Plane(getNormal(ONE), pt);
00904           plane[TWO]  = Plane(getNormal(TWO), pt);
00905 
00906           const double halfSphVol = sphVol*0.5;
00907           double s[tmplDim];
00908           s[ZERO] = 0;
00909           s[ONE]  = getSphere().getSegmentVolume(plane[ONE]);
00910           s[TWO]  = getSphere().getSegmentVolume(plane[TWO]);
00911           s[ONE] = ((s[ONE] > halfSphVol) ? (sphVol - s[ONE]) : s[ONE]);
00912           s[TWO] = ((s[TWO] > halfSphVol) ? (sphVol - s[TWO]) : s[TWO]);
00913 
00914           Vec distVec(4.0*getSphere().getRadius());
00915           distVec[ONE] = plane[ONE].getDistanceTo(centrePt);
00916           distVec[TWO] = plane[TWO].getDistanceTo(centrePt);
00917           const double coreVol =
00918             2.0*getSphere().getVolume(
00919               centrePt - Vec(distVec[ONE], distVec[TWO], distVec[ZERO]),
00920               centrePt + Vec(distVec[ONE], distVec[TWO], distVec[ZERO])
00921             );
00922           double v[tmplDim];
00923           v[ONE]  = 0.50*(sphVol - 2.0*s[TWO] - coreVol);
00924           v[TWO]  = 0.50*(sphVol - 2.0*s[ONE] - coreVol);
00925           v[ZERO] = 0.25*(sphVol - 2.0*v[ONE] - 2.0*v[TWO] - coreVol);
00926 
00927           if (pt[ONE] <= centrePt[ONE])
00928           {
00929             if (pt[TWO] <= centrePt[TWO])
00930             {
00931               vol = coreVol + v[ONE] + v[TWO] + v[ZERO];
00932             }
00933             else
00934             {
00935               vol = v[TWO] + v[ZERO];
00936             }
00937           }
00938           else
00939           {
00940             if (pt[TWO] <= centrePt[TWO])
00941             {
00942               vol = v[ONE] + v[ZERO];
00943             }
00944             else
00945             {
00946               vol = v[ZERO];
00947             }
00948           }
00949         }
00950         else
00951         {
00952           if (pt[ONE] <= centrePt[ONE])
00953           {
00954             if (pt[TWO] <= centrePt[TWO])
00955             {
00956               vol = 
00957                 sphVol
00958                 -
00959                 getSphere().getSegmentVolume(Plane(getNegNormal(ONE), pt))
00960                 -
00961                 getSphere().getSegmentVolume(Plane(getNegNormal(TWO), pt));
00962             }
00963             else
00964             {
00965               vol = getSphere().getSegmentVolume(Plane(getNormal(TWO), pt));
00966             }
00967           }
00968           else
00969           {
00970             if (pt[TWO] <= centrePt[TWO])
00971             {
00972               vol = getSphere().getSegmentVolume(Plane(getNormal(ONE), pt));
00973             }
00974             else
00975             {
00976               vol = 0.0;
00977             }
00978           }
00979         }
00980 
00981         return vol;
00982       }
00983 
00984       template <int tmplDim, typename TmplVec>
00985       double IntersectionVolCalculator<tmplDim,TmplVec>::getOutsidePointVolume(
00986         const Vec &pt
00987       ) const
00988       {
00989         double vol = 0.0;
00990         const double sphVol = getSphere().getVolume();
00991         const Vec &centrePt = getSphere().getCentre();
00992         if (tmplDim == 2)
00993         {
00994           if (pt[0] <= centrePt[0])
00995           {
00996             if (pt[1] <= centrePt[1])
00997             {
00998               vol = 
00999                 sphVol
01000                 -
01001                 getSphere().getSegmentVolume(Plane(getNegNormal(0), pt))
01002                 -
01003                 getSphere().getSegmentVolume(Plane(getNegNormal(1), pt));
01004             }
01005             else
01006             {
01007               vol = getSphere().getSegmentVolume(Plane(getNormal(1), pt));
01008             }
01009           }
01010           else
01011           {
01012             if (pt[1] <= centrePt[1])
01013             {
01014               vol = getSphere().getSegmentVolume(Plane(getNormal(0), pt));
01015             }
01016             else
01017             {
01018               vol = 0.0;
01019             }
01020           }
01021         }
01022         else if (tmplDim == 3)
01023         {
01024           const Vec diag = (centrePt - pt)*2.0;
01025           const Vec oppCorner = pt + diag;
01026           BasicBox box =
01027             BasicBox(
01028               componentMin(pt, oppCorner),
01029               componentMax(pt, oppCorner)
01030             );
01031 
01032           double s[tmplDim];
01033           double e[tmplDim];
01034           for (int i = 0; i < tmplDim; i++)
01035           {
01036             s[i] = getSphere().getSegmentVolume(Plane(getNormal(i), box.getMaxPt()));
01037             e[i] = getTwoPlaneVolume(box.getMaxPt(), i);
01038           }
01039           double v[tmplDim+1];
01040           v[0] = s[0] - 2.0*e[1] - 2.0*e[2];
01041           v[1] = s[1] - 2.0*e[0] - 2.0*e[2];
01042           v[2] = s[2] - 2.0*e[0] - 2.0*e[1];
01043           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]);
01044 
01045           if (pt[0] <= centrePt[0])
01046           {
01047             if (pt[1] <= centrePt[1])
01048             {
01049               if (pt[2] <= centrePt[2])
01050               {
01051                 vol = v[0] + v[1] + v[2] + v[3] + e[0] + e[1] + e[2];
01052               }
01053               else
01054               {
01055                 vol = v[2] + e[0] + e[1];
01056               }
01057             }
01058             else
01059             {
01060               if (pt[2] <= centrePt[2])
01061               {
01062                 vol = v[1] + e[0] + e[2];
01063               }
01064               else
01065               {
01066                 vol = e[0];
01067               }
01068             }
01069           }
01070           else
01071           {
01072             if (pt[1] <= centrePt[1])
01073             {
01074               if (pt[2] <= centrePt[2])
01075               {
01076                 vol = v[0] + e[1] + e[2];
01077               }
01078               else
01079               {
01080                 vol = e[1];
01081               }
01082             }
01083             else
01084             {
01085               if (pt[2] <= centrePt[2])
01086               {
01087                 vol = e[2];
01088               }
01089               else
01090               {
01091                 vol = 0.0;
01092               }
01093             }
01094           }
01095         }
01096         return vol;
01097       }
01098 
01099       
01100       template <int tmplDim, typename TmplVec>
01101       double IntersectionVolCalculator<tmplDim,TmplVec>::getVolume(
01102         const Vertex &vtx
01103       )
01104       {
01105         double vol = 0;
01106         if (getSphere().intersectsWith(vtx.getPoint()))
01107         {
01108           vol = getInsidePointVolume(vtx.getPoint());
01109         }
01110         else
01111         {
01112           vol = getOutsidePointVolume(vtx.getPoint());
01113         }
01114         return vol;
01115       }
01116 
01117       template <int tmplDim, typename TmplVec>
01118       double IntersectionVolCalculator<tmplDim,TmplVec>::getVertexVolume(
01119         const BasicSphere &sphere
01120       )
01121       {
01122         m_sphere = VolumeSphere(sphere);
01123         double vol = 0.0;
01124 
01125         double vtxVol[getVertexBox().getNumVertices()];
01126         for (int i = 0; i < getVertexBox().getNumVertices(); i++)
01127         {
01128           vtxVol[i] = getVolume(getVertexBox().getVertex(i));
01129         }
01130 
01131         if (tmplDim == 2)
01132         {
01133           vol = vtxVol[0] - vtxVol[1] - vtxVol[2] + vtxVol[3];
01134         }
01135         else if (tmplDim == 3)
01136         {
01137           vol = 
01138             vtxVol[7] + vtxVol[6] + vtxVol[5] - vtxVol[4]
01139             -
01140             vtxVol[3] - vtxVol[2] - vtxVol[1] + vtxVol[0];
01141         }
01142 
01143         return vol;
01144       }
01145 
01146       template <int tmplDim, typename TmplVec>
01147       bool IntersectionVolCalculator<tmplDim,TmplVec>::sphereContainsBox(
01148         const BasicSphere &sphere
01149       ) const
01150       {
01151         for (int i = 0; i < getVertexBox().getNumVertices(); i++)
01152         {
01153           if (!sphere.intersectsWith(getVertexBox().getVertex(i).getPoint()))
01154           {
01155             return false;
01156           }
01157         }
01158         return true;
01159       }
01160 
01161       template <int tmplDim, typename TmplVec>
01162       double IntersectionVolCalculator<tmplDim,TmplVec>::getVolume(
01163         const BasicSphere &sphere
01164       )
01165       {
01166         double vol = 0.0;
01167         if (getBox().intersectsWith(sphere))
01168         {
01169           if (sphereContainsBox(sphere))
01170           {
01171             vol = getBox().getVolume();
01172           }
01173           else if (getBox().contains(sphere))
01174           {
01175             vol = sphere.getVolume();
01176           }
01177           else
01178           {
01179             vol = getVertexVolume(sphere);
01180           }
01181         }
01182         return vol;
01183       }
01184     }
01185   }
01186 }