1 // This is brl/bbas/bvgl/bvgl_volume_of_intersection.hxx 2 #ifndef bvgl_volume_of_intersection_hxx_ 3 #define bvgl_volume_of_intersection_hxx_ 4 //: 5 // \file 6 // \author Andrew Miller 7 8 #include <limits> 9 #include <cmath> 10 #include <iostream> 11 #include <algorithm> 12 #include <vector> 13 #include "bvgl_volume_of_intersection.h" 14 15 #include <cassert> 16 #ifdef _MSC_VER 17 # include <vcl_msvc_warnings.h> 18 #endif 19 #include <vnl/vnl_math.h> 20 #include <vgl/vgl_sphere_3d.h> 21 22 //tolerance constants/functions 23 static double bvgl_eps = 1.0e-8; // tolerance for intersections bvgl_near_zero(double x)24inline bool bvgl_near_zero(double x) { return x < bvgl_eps && x > -bvgl_eps; } bvgl_near_eq(double x,double y)25inline bool bvgl_near_eq(double x, double y) { return bvgl_near_zero(x-y); } 26 27 //: calculate the volume of intersection between these two spheres 28 template <typename T> bvgl_volume_of_intersection(vgl_sphere_3d<T> const & A,vgl_sphere_3d<T> const & B)29T bvgl_volume_of_intersection(vgl_sphere_3d<T> const& A, vgl_sphere_3d<T> const& B) 30 { 31 //distance between the two spheres 32 T d = (A.centre() - B.centre()).length(); 33 T r0 = A.radius(); 34 T r1 = B.radius(); 35 36 //cases 37 //0. if one sphere is completely inside the other one 38 T difR = std::fabs(r0 - r1); 39 if ( d <= difR ) 40 { 41 //return the volume of the smaller sphere 42 T minR = std::min(r0, r1); 43 return (4.0/3.0) * vnl_math::pi * minR * minR * minR; 44 } 45 46 //2. The two spheres intersect... 47 T sumR = r0 + r1; 48 if (d > difR && d < sumR) 49 { 50 //calculate distance to the plane of intersection from the center of circle A 51 T xInt = ( r0*r0 - r1*r1 + d*d ) / (2.0*d) ; 52 53 //calculate height of sperhical caps 54 T h0 = r0 - xInt; 55 T h1 = r1 + xInt - d; 56 57 // volume is sum of tops 58 return (vnl_math::pi/3.0)*( h0*h0 * (3*r0-h0) + h1*h1 * (3*r1-h1) ); 59 } 60 61 //otherwise the two spheres do not intersect 62 return (T) 0.0; 63 } 64 65 66 //--------------------------------------------------------------------------- 67 //Template instantiation 68 #undef BVGL_VOLUME_OF_INTERSECTION_INSTANTIATE 69 #define BVGL_VOLUME_OF_INTERSECTION_INSTANTIATE(T) \ 70 template T bvgl_volume_of_intersection(vgl_sphere_3d<T > const& A, vgl_sphere_3d<T > const& B) 71 72 #endif // bvgl_volume_of_intersection_hxx_ 73