vcglib: sphere3.h Source File

00001 /****************************************************************************
00002 * VCGLib                                                            o o     *
00003 * Visual and Computer Graphics Library                            o     o   *
00004 *                                                                _   O  _   *
00005 * Copyright(C) 2004                                                \/)\/    *
00006 * Visual Computing Lab                                            /\/|      *
00007 * ISTI - Italian National Research Council                           |      *
00008 *                                                                    \      *
00009 * All rights reserved.                                                      *
00010 *                                                                           *
00011 * This program is free software; you can redistribute it and/or modify      *   
00012 * it under the terms of the GNU General Public License as published by      *
00013 * the Free Software Foundation; either version 2 of the License, or         *
00014 * (at your option) any later version.                                       *
00015 *                                                                           *
00016 * This program is distributed in the hope that it will be useful,           *
00017 * but WITHOUT ANY WARRANTY; without even the implied warranty of            *
00018 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
00019 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt)          *
00020 * for more details.                                                         *
00021 *                                                                           *
00022 ****************************************************************************/
00023 /****************************************************************************
00024   History
00025 
00026 $Log: not supported by cvs2svn $
00027 Revision 1.9  2006/07/06 12:40:34  ganovelli
00028 typdef ..ScalarType  added
00029 
00030 Revision 1.8  2005/02/22 14:18:15  ponchio
00031 assert addded.
00032 
00033 Revision 1.7  2005/02/21 17:03:03  ponchio
00034 Added Tight creation.
00035 
00036 Revision 1.6  2004/12/01 16:06:59  ponchio
00037 Distance
00038 
00039 Revision 1.5  2004/09/29 13:55:33  ponchio
00040 Added Distance shpere - point.
00041 
00042 Revision 1.4  2004/04/02 09:49:01  ponchio
00043 Ehm... a couople of small errors.
00044 
00045 Revision 1.3  2004/04/02 09:44:13  ponchio
00046 Sphere ->Sphere3
00047 
00048 Revision 1.2  2004/03/25 17:25:46  ponchio
00049 #include sbagliato.
00050 
00051 Revision 1.1  2004/03/21 17:51:57  ponchio
00052 First version.
00053 
00054 
00055 ****************************************************************************/
00056 
00057 #ifndef VCG_SPHERE_H
00058 #define VCG_SPHERE_H
00059 
00060 #include <assert.h>
00061 #include <vcg/space/point3.h>
00062 #include <vector>
00063 
00064 namespace vcg {
00065 
00075 template <class T> class Sphere3 {
00076 protected:
00077   Point3<T> _center;
00078   T _radius;
00079  public:
00080         typedef T ScalarType;
00081   Sphere3(): _radius(-1) {}
00082   Sphere3(const Point3<T> &center, T radius): _center(center), _radius(radius) {}
00083   
00084   T &Radius() { return _radius; }
00085   const T &Radius() const { return _radius; }
00086   Point3<T> &Center() { return _center; }
00087   const Point3<T> &Center() const { return _center; }
00088 
00089         bool IsEmpty() const { return _radius < 0; }
00091   bool IsIn(const Point3<T> &p) const;
00092   bool IsIn(const Sphere3<T> &p) const;
00093 
00094   void Add(const Point3<T> &p);
00095   void Add(const Sphere3 &sphere);
00096   void Intersect(const Sphere3 &sphere);
00097         
00098 
00099   int CreateFromBox(int n, const Point3<T> *points);
00100   //makes 36 iterations over the data... but get good results.
00101         int CreateTight(int n, const Point3<T> *points, 
00102                   T threshold = 1.01, T speed = 0.6);
00103         int CreateTight(const std::vector<Point3< T> > & points, 
00104                   T threshold = 1.01, T speed = 0.6);
00105 };
00106 
00107 template <class T> T Distance(const Sphere3<T> &sphere, 
00108                               const Point3<T> &point) {
00109   T dist = Distance(point, sphere.Center()) - sphere.Radius();
00110   if(dist < 0) dist = 0;
00111   return dist;
00112 }
00113 
00114 template <class T> T Distance(const Sphere3<T> &sphere, 
00115                               const Sphere3<T> &s) {
00116   T dist = Distance(s.Center(), sphere.Center()) 
00117                     - sphere.Radius() - s.Radius();
00118   if(dist < 0) dist = 0;
00119   return dist;
00120 }
00121 
00122 typedef Sphere3<float> Sphere3f;
00123 typedef Sphere3<double> Sphere3d;
00124 
00125 template <class T> void Sphere3<T>::Add(const Sphere3<T> &sphere) {
00126   if(IsEmpty()) {
00127                         *this = sphere;
00128                         return;
00129                 }
00130                 Point3<T> dist = sphere.Center() - _center;
00131     float distance = dist.Norm();
00132     float fartest = sphere.Radius() + distance;
00133                 if(fartest <= _radius) 
00134       return;
00135 
00136     float nearest = sphere.Radius() - distance;
00137     if(nearest >= _radius) {
00138       *this = sphere;
00139       return;
00140     }
00141 
00142     if(distance < 0.001*(_radius + sphere.Radius())) {
00143       _radius += distance;
00144       return;
00145     }
00146 
00147     _center += dist * ((fartest - _radius) / (distance * 2));
00148     _radius = (_radius + fartest)/2;
00149 }
00150 
00151 template <class T> void Sphere3<T>::Add(const Point3<T> &p) {
00152         if(IsEmpty()) {
00153                 _center = p;
00154                 _radius = 0;
00155         }
00156         Point3<T> dist = p - _center;
00157         float fartest = dist.Norm();
00158         if(fartest <= _radius) return;
00159         _center += dist * ((fartest - _radius) / (fartest*2));
00160         _radius = (_radius + fartest)/2;
00161 }  
00162 
00163 template <class T> bool Sphere3<T>::IsIn(const Point3<T> &p) const {
00164   if(IsEmpty()) return false;
00165   Point3<T> dist = p - _center;
00166   return dist.SquaredNorm() <= _radius*_radius;
00167 }
00168 
00169 template <class T> bool Sphere3<T>::IsIn(const Sphere3<T> &p) const {
00170   if(IsEmpty()) return false;
00171   Point3<T> dist = p.Center() - _center;
00172   float distance = dist.Norm();
00173   return distance + p.Radius() < _radius;
00174 }
00175 
00176 template <class T> void Sphere3<T>::Intersect(const Sphere3<T> &s) {
00177   float dist = Distance(_center, s.Center());
00178   float r = 0.5 * (_radius + s.Radius() - dist);
00179   if(r < 0) {
00180     _radius = -1;
00181     return;
00182   }
00183   _center = (s.Center()*(_radius - r) + _center*(s.Radius() - r))/dist;
00184   _radius = r;
00185 }
00186 
00187  template <class T> int Sphere3<T>::CreateFromBox(int n,  const Point3<T> *points) {
00188    Point3f max(-1e100, -1e100, -1e100);
00189    Point3f min(1e100, 1e100, 1e100);
00190    for(int i = 0; i < n; i++) {
00191      for(int k = 0; k < 3; k++) {
00192        if(max[k] < points[i][k]) max[k] = points[i][k];
00193        if(min[k] > points[i][k]) min[k] = points[i][k];
00194      }
00195    }
00196    Center() = (min + max)/2;
00197    Radius() = (min - max).Norm()/2;
00198    return 0;
00199  }
00200  template <class T> int Sphere3<T>::CreateTight(int n, const Point3<T> *points, 
00201                                                 T threshold, T speed) { 
00202    //This is quantized gradient descent... really ugly. But simple :P
00203    //TODO step should adapt to terrain...
00204    for(int i = 0; i < n; i++)
00205      Add(points[i]);
00206    Radius() *= 1.0001;
00207 
00208    Point3<T> center;
00209    //Test with 6 directions
00210 
00211    Point3f pert[6];
00212    T step = Radius()/8;
00213    int count = 0;
00214    while(1) {
00215      count++;
00216      T radius = Radius();
00217      pert[0] = Point3f(step, 0, 0);
00218      pert[1] = -pert[0];
00219      pert[2] = Point3f(0, step, 0);
00220      pert[3] = -pert[2];
00221      pert[4] = Point3f(0, 0, step);
00222      pert[5] = -pert[4];
00223 
00224      int best = 6;
00225      T best_radius = Radius()/threshold;
00226 
00227      for(int k = 0; k < 6; k++) {
00228        center = Center() + pert[k];
00229        radius = 0;
00230        for(int i = 0; i < n; i++) {
00231          float r = Distance(center, points[i]);
00232          if(r > radius)
00233            radius = r;
00234        }
00235        if(radius < best_radius) {
00236          best = k;
00237          best_radius = radius;
00238        }
00239      }
00240      if(best != 6) {
00241        Center() = Center() + pert[best];
00242        Radius() = best_radius;
00243      }
00244      step *= speed; 
00245      if(step <= Radius() * (threshold - 1))
00246        break;
00247    }
00248    Radius() *= 1.01;
00249    
00250    //Test we did it correctly.
00251    for(int i = 0; i < n; i++) 
00252      assert(IsIn(points[i]));
00253    
00254    return count;
00255  }
00256 
00257 template <class T> int Sphere3<T>::CreateTight(const std::vector<Point3<T> > & points,
00258                                                                                                                                                                                          T threshold, T speed){
00259                 return (points.empty())? -1 :CreateTight(points.size(),&(*points.begin()),threshold,speed);
00260 }
00261 
00262 } //namespace
00263 
00264 
00265 #endif