WFMath 0.3.12
ball_funcs.h
00001 // ball_funcs.h (n-dimensional ball implementation)
00002 //
00003 //  The WorldForge Project
00004 //  Copyright (C) 2000, 2001  The WorldForge Project
00005 //
00006 //  This program is free software; you can redistribute it and/or modify
00007 //  it under the terms of the GNU General Public License as published by
00008 //  the Free Software Foundation; either version 2 of the License, or
00009 //  (at your option) any later version.
00010 //
00011 //  This program is distributed in the hope that it will be useful,
00012 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014 //  GNU General Public License for more details.
00015 //
00016 //  You should have received a copy of the GNU General Public License
00017 //  along with this program; if not, write to the Free Software
00018 //  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00019 //
00020 //  For information about WorldForge and its authors, please contact
00021 //  the Worldforge Web Site at http://www.worldforge.org.
00022 //
00023 
00024 // Author: Ron Steinke
00025 
00026 #ifndef WFMATH_BALL_FUNCS_H
00027 #define WFMATH_BALL_FUNCS_H
00028 
00029 #include <wfmath/ball.h>
00030 
00031 #include <wfmath/axisbox.h>
00032 #include <wfmath/miniball.h>
00033 
00034 #include <cassert>
00035 
00036 namespace WFMath {
00037 
00038 template<int dim>
00039 inline bool Ball<dim>::isEqualTo(const Ball<dim>& b, double epsilon) const
00040 {
00041   return Equal(m_center, b.m_center, epsilon)
00042       && Equal(m_radius, b.m_radius, epsilon);
00043 }
00044 
00045 template<int dim>
00046 AxisBox<dim> Ball<dim>::boundingBox() const
00047 {
00048   Point<dim> p_low, p_high;
00049 
00050   for(int i = 0; i < dim; ++i) {
00051     p_low[i] = m_center[i] - m_radius;
00052     p_high[i] = m_center[i] + m_radius;
00053   }
00054 
00055   bool valid = m_center.isValid();
00056 
00057   p_low.setValid(valid);
00058   p_high.setValid(valid);
00059 
00060   return AxisBox<dim>(p_low, p_high, true);
00061 }
00062 
00063 template<int dim, template<class, class> class container>
00064 Ball<dim> BoundingSphere(const container<Point<dim>, std::allocator<Point<dim> > >& c)
00065 {
00066   _miniball::Miniball<dim> m;
00067   _miniball::Wrapped_array<dim> w;
00068 
00069   typename container<Point<dim>, std::allocator<Point<dim> > >::const_iterator i, end = c.end();
00070   bool valid = true;
00071 
00072   for(i = c.begin(); i != end; ++i) {
00073     valid = valid && i->isValid();
00074     for(int j = 0; j < dim; ++j)
00075       w[j] = (*i)[j];
00076     m.check_in(w);
00077   }
00078 
00079   m.build();
00080 
00081 #ifndef NDEBUG
00082   double dummy;
00083 #endif
00084   assert("Check that bounding sphere is good to library accuracy" &&
00085          m.accuracy(dummy) < WFMATH_EPSILON);
00086 
00087   w = m.center();
00088   Point<dim> center;
00089 
00090   for(int j = 0; j < dim; ++j)
00091     center[j] = w[j];
00092 
00093   center.setValid(valid);
00094 
00095   return Ball<dim>(center, std::sqrt(m.squared_radius()));
00096 }
00097 
00098 template<int dim, template<class, class> class container>
00099 Ball<dim> BoundingSphereSloppy(const container<Point<dim>, std::allocator<Point<dim> > >& c)
00100 {
00101   // This is based on the algorithm given by Jack Ritter
00102   // in Volume 2, Number 4 of Ray Tracing News
00103   // <http://www.acm.org/tog/resources/RTNews/html/rtnews7b.html>
00104 
00105   typename container<Point<dim>, std::allocator<Point<dim> > >::const_iterator i = c.begin(),
00106                                                 end = c.end();
00107   if (i == end) {
00108     return Ball<dim>();
00109   }
00110 
00111   CoordType min[dim], max[dim];
00112   typename container<Point<dim>, std::allocator<Point<dim> > >::const_iterator min_p[dim], max_p[dim];
00113   bool valid = i->isValid();
00114 
00115   for(int j = 0; j < dim; ++j) {
00116     min[j] = max[j] = (*i)[j];
00117     min_p[j] = max_p[j] = i;
00118   }
00119 
00120   while(++i != end) {
00121     valid = valid && i->isValid();
00122     for(int j = 0; j < dim; ++j) {
00123       if(min[j] > (*i)[j]) {
00124         min[j] = (*i)[j];
00125         min_p[j] = i;
00126       }
00127       if(max[j] < (*i)[j]) {
00128         max[j] = (*i)[j];
00129         max_p[j] = i;
00130       }
00131     }
00132   }
00133 
00134   CoordType span = -1;
00135   int direction = -1;
00136 
00137   for(int j = 0; j < dim; ++j) {
00138     CoordType new_span = max[j] - min[j];
00139     if(new_span > span) {
00140       span = new_span;
00141       direction = j;
00142     }
00143   }
00144 
00145   assert("Have a direction of maximum size" && direction != -1);
00146 
00147   Point<dim> center = Midpoint(*(min_p[direction]), *(max_p[direction]));
00148   CoordType dist = SloppyDistance(*(min_p[direction]), center);
00149 
00150   for(i = c.begin(); i != end; ++i) {
00151     if(i == min_p[direction] || i == max_p[direction])
00152       continue; // We already have these
00153 
00154     CoordType new_dist = SloppyDistance(*i, center);
00155 
00156     if(new_dist > dist) {
00157       CoordType delta_dist = (new_dist - dist) / 2;
00158       // Even though new_dist may be too large, delta_dist / new_dist
00159       // always gives enough of a shift to include the new point.
00160       center += (*i - center) * delta_dist / new_dist;
00161       dist += delta_dist;
00162       assert("Shifted ball contains new point" &&
00163              SquaredDistance(*i, center) <= dist * dist);
00164     }
00165   }
00166 
00167   center.setValid(valid);
00168 
00169   return Ball<2>(center, dist);
00170 }
00171 
00172 // These two are here, instead of defined in the class, to
00173 // avoid include order problems
00174 
00175 template<int dim>
00176 inline Ball<dim> Point<dim>::boundingSphere() const
00177 {
00178   return Ball<dim>(*this, 0);
00179 }
00180 
00181 template<int dim>
00182 inline Ball<dim> Point<dim>::boundingSphereSloppy() const
00183 {
00184   return Ball<dim>(*this, 0);
00185 }
00186 
00187 } // namespace WFMath
00188 
00189 #endif  // WFMATH_BALL_FUNCS_H