ETL/trunk/ETL/_bspline.h

   1 /*! ========================================================================
   2 ** Extended Template and Library
   3 ** B-Spline Class Implementation
   4 ** $Id$
   5 **
   6 ** Copyright (c) 2002 Robert B. Quattlebaum Jr.
   7 **
   8 ** This package is free software; you can redistribute it and/or
   9 ** modify it under the terms of the GNU General Public License as
  10 ** published by the Free Software Foundation; either version 2 of
  11 ** the License, or (at your option) any later version.
  12 **
  13 ** This package is distributed in the hope that it will be useful,
  14 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16 ** General Public License for more details.
  17 **
  18 ** === N O T E S ===========================================================
  19 **
  20 ** This is an internal header file, included by other ETL headers.
  21 ** You should not attempt to use it directly.
  22 **
  23 ** ========================================================================= */
  24
  25 /* === S T A R T =========================================================== */
  26
  27 #ifndef __ETL__BSPLINE_H
  28 #define __ETL__BSPLINE_H
  29
  30 /* === H E A D E R S ======================================================= */
  31
  32 #include <vector>
  33 #include <functional>
  34 #include "_curve_func.h"
  35
  36 /* === M A C R O S ========================================================= */
  37
  38 /* === T Y P E D E F S ===================================================== */
  39
  40 /* === C L A S S E S & S T R U C T S ======================================= */
  41
  42 _ETL_BEGIN_NAMESPACE
  43
  44 template <class T, class K=float, class C=affine_combo<T,K>, class D=distance_func<T> >
  45 class bspline : public std::unary_function<K,T>
  46 {
  47 public:
  48         typedef T value_type;
  49         typedef K knot_type;
  50         typedef std::vector<knot_type>  knot_container;
  51         typedef std::vector<value_type> cpoint_container;
  52         typedef typename knot_container::iterator knot_iterator;
  53         typedef typename cpoint_container::iterator cpoint_iterator;
  54
  55         typedef C affine_func_type;
  56         typedef D distance_func_type;
  57
  58 protected:
  59         affine_func_type affine_func;
  60         distance_func_type distance_func;
  61
  62 private:
  63         int m;
  64         knot_container _knots;
  65         cpoint_container _cpoints;
  66         bool _loop;
  67
  68 public:
  69         bspline():m(2),_loop(false) {  }
  70
  71         int get_m()const { return m-1; };
  72         int set_m(int new_m) { m=new_m+1; return m-1; };
  73
  74         bool set_loop(bool x) { _loop=x; reset_knots(); return _loop; }
  75
  76         knot_container & knots() { return _knots; };
  77         cpoint_container & cpoints() { return _cpoints; };
  78
  79         const knot_container & knots()const { return _knots; };
  80         const cpoint_container & cpoints()const { return _cpoints; };
  81
  82         void reset_knots()
  83         {
  84                 int i;
  85
  86                 if(!_loop)
  87                 {
  88                         _knots.clear();
  89                         if(!_cpoints.size())
  90                                 return;
  91                         while(m>(signed)_cpoints.size())
  92                                 m--;
  93                         for(i=0;i<m;i++)
  94                                 *_knots.insert(_knots.end())=0;
  95                         for(i=1;i<(signed)_cpoints.size()-m+1;i++)
  96                                 *_knots.insert(_knots.end())=i;
  97                         for(i=0;i<m;i++)
  98                                 *_knots.insert(_knots.end())=_cpoints.size()-m+1;
  99                 }
 100                 else
 101                 {
 102                         _knots.clear();
 103                         if(!_cpoints.size())
 104                                 return;
 105                         while(m>(signed)_cpoints.size())
 106                                 m--;
 107                         for(i=0;i<=(signed)_cpoints.size()-m+1;i++)
 108                                 *_knots.insert(_knots.end())=i;
 109                 }
 110         }
 111
 112         int calc_curve_segment(knot_type t)const
 113         {
 114                 int k;
 115                 if(t<0)
 116                         t=0;
 117                 if(t>=_knots.back())
 118                         t=_knots.back()-0.0001;
 119                 for(k=0;_knots[k]>t || _knots[k+1]<=t;k++)
 120                   ;
 121
 122                 return k;
 123         }
 124
 125         knot_container get_segment_knots(int i)const
 126         {
 127                 if(i+1<m)
 128                 {
 129                         knot_container ret(_knots.begin(),_knots.begin()+i+m+1);
 130
 131                         return ret;
 132                 }
 133                 if(i+1>=(signed)_knots.size())
 134                 {
 135                         knot_container ret(_knots.begin()+i-m+1,_knots.end());
 136                         return ret;
 137                 }
 138                 return knot_container(_knots.begin()+i-m+1,     _knots.begin()+i+m);
 139         }
 140
 141         cpoint_container get_segment_cpoints(int i)const
 142         {
 143                 if(i+1<m)
 144                 {
 145                         return cpoint_container();
 146                 }
 147                 if(i+1>=(signed)_knots.size())
 148                 {
 149                         return cpoint_container();
 150                 }
 151                 return cpoint_container(_cpoints.begin()+i-m+1, _cpoints.begin()+i+1);
 152         }
 153
 154         cpoint_container calc_shell(knot_type t, int level)const
 155         {
 156                 int
 157                         i=calc_curve_segment(t),
 158                         j,k;
 159
 160                 knot_container u=get_segment_knots(i);
 161
 162                 cpoint_container d=get_segment_cpoints(i);
 163
 164                 if(!d.size())
 165                         return cpoint_container();
 166
 167                 for(j=0;d.size()>1 && j<level;d.pop_back(),j++)
 168                 {
 169                         for(k=0;k<d.size()-1;k++)
 170                         {
 171                                 d[k]=affine_func(d[k],d[k+1],((t-u[j+k+1])/(u[m+k]-u[j+k+1])));
 172                         }
 173                 }
 174                 return d;
 175         }
 176
 177         value_type operator()(knot_type t)const
 178         {
 179                 return get_curve_val(calc_curve_segment(t),t);
 180         }
 181
 182         value_type get_curve_val(int i,knot_type t)const
 183         {
 184                 int
 185                         j,k;
 186
 187                 knot_container u=get_segment_knots(i);
 188
 189                 cpoint_container d=get_segment_cpoints(i);
 190
 191                 if(!d.size())
 192                         return value_type();
 193
 194                 for(j=0;d.size()>1;d.pop_back(),j++)
 195                 {
 196                         for(k=0;k<(signed)d.size()-1;k++)
 197                         {
 198                                 d[k]=affine_func(d[k],d[k+1],((t-u[j+k+1])/(u[m+k]-u[j+k+1])));
 199                         }
 200                 }
 201                 return d.front();
 202         }
 203
 204         cpoint_iterator find_closest_cpoint(const value_type &point, typename distance_func_type::result_type max)
 205         {
 206                 cpoint_iterator i=_cpoints.begin();
 207                 cpoint_iterator ret=i;
 208                 typename distance_func_type::result_type dist=distance_func(point,_cpoints[0]);
 209
 210                 // The distance function returns "cooked" (ie: squared)
 211                 // distances, so we need to cook our max distance for
 212                 // the comparison to work correctly.
 213                 max=distance_func.cook(max);
 214
 215                 for(++i;i<_cpoints.end();i++)
 216                 {
 217                         typename distance_func_type::result_type thisdist=distance_func(point,*i);
 218
 219                         if(thisdist<dist)
 220                         {
 221                                 dist=thisdist;
 222                                 ret=i;
 223                         }
 224                 }
 225                 if(dist<max)
 226                         return ret;
 227                 return _cpoints.end();
 228         }
 229 };
 230
 231 _ETL_END_NAMESPACE
 232
 233 /* -- F U N C T I O N S ----------------------------------------------------- */
 234
 235 /* -- E N D ----------------------------------------------------------------- */
 236
 237 #endif