UnifBsplines.h

00001 //===========================================================================
00002 // SINTEF Multilevel B-spline Approximation library - version 1.1
00003 //
00004 // Copyright (C) 2000-2005 SINTEF ICT, Applied Mathematics, Norway.
00005 //
00006 // This program is free software; you can redistribute it and/or          
00007 // modify it under the terms of the GNU General Public License            
00008 // as published by the Free Software Foundation version 2 of the License. 
00009 //
00010 // This program is distributed in the hope that it will be useful,        
00011 // but WITHOUT ANY WARRANTY; without even the implied warranty of         
00012 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          
00013 // GNU General Public License for more details.                           
00014 //
00015 // You should have received a copy of the GNU General Public License      
00016 // along with this program; if not, write to the Free Software            
00017 // Foundation, Inc.,                                                      
00018 // 59 Temple Place - Suite 330,                                           
00019 // Boston, MA  02111-1307, USA.                                           
00020 //
00021 // Contact information: e-mail: tor.dokken@sintef.no                      
00022 // SINTEF ICT, Department of Applied Mathematics,                         
00023 // P.O. Box 124 Blindern,                                                 
00024 // 0314 Oslo, Norway.                                                     
00025 //
00026 // Other licenses are also available for this software, notably licenses
00027 // for:
00028 // - Building commercial software.                                        
00029 // - Building software whose source code you wish to keep private.        
00030 //===========================================================================
00031 #ifndef _UNIFBSPLINES_H_
00032 #define _UNIFBSPLINES_H_
00033 
00034 #include <GenMatrix.h>
00035 
00036 #include <cmath>
00037 using namespace std;
00038 
00039 //#define MBA_CUBIC_C1 1 
00040 
00041 #ifdef  MBA_CUBIC_C1
00042 // Basis functions
00043 inline double B_0(double t) {double tmp = 1.0-t; return tmp*tmp*tmp/2.0;}
00044 inline double B_1(double t) {return (2.5*t*t*t - 4.5*t*t + 1.5*t + 0.5);}
00045 inline double B_2(double t) {return (-2.5*t*t*t + 3.0*t*t);}
00046 inline double B_3(double t) {return t*t*t/2.0;}
00047 
00048 // And their derivatives
00049 inline double dB_0(double t) {return 1.5*(1.0-t)*(t-1.0);}
00050 inline double dB_1(double t) {return (7.5*t*t - 9.0*t + 1.5);}
00051 inline double dB_2(double t) {return (-7.5*t*t + 6.0*t);}
00052 inline double dB_3(double t) {return 1.5*t*t;}
00053 
00054 // PRE EVALUATTED tensor products of B-splines, in a grid index
00055 // -------------------------------------------------------------
00056 
00057 // Pre-evaluated tensors equivalent to B(k,s)*B(l,t)
00058 static double tensor_BB[2][2] = {
00059   {1./4., 1./4.},
00060   {1./4., 1./4.}
00061 };
00062 
00063 // ???? only temp. for compilation
00064 // equivalent to dB(k,s)*B(l,t).
00065 static double tensor_dBB[2][2] = {
00066   {-3./4., -3./4.},
00067   { 3./4.,  3./4.}
00068 };
00069 
00070 // The transpose of the matrix above, for the y-components,
00071 // equivalent to B(k,s)*dB(l,t)
00072 static double tensor_BdB[2][2] = {
00073   {-3./4., 3./4.},
00074   {-3./4., 3./4.}
00075 };
00076 
00077 #else
00078 // The uniform (cardinal), univariate cubic C2 basis functions
00079 // 0 <= t < 1;
00080 inline double B_0(double t) {double tmp = 1.0-t; return tmp*tmp*tmp/6.0;}
00081 
00082 //static inline double B_1(double t) {return (3.0*t*t*t - 6.0*t*t+4.0)/6.0;}
00083 static double div46 = 4.0/6.0;
00084 inline double B_1(double t) {return (0.5*t*t*t - t*t + div46);}
00085 
00086 //static inline double B_2(double t) {return (-3.0*t*t*t + 3.0*t*t + 3.0*t + 1.0)/6.0;}
00087 static double div16 = 1.0/6.0;
00088 inline double B_2(double t) {return (-0.5*t*t*t + 0.5*t*t + 0.5*t + div16);}
00089 
00090 inline double B_3(double t) {return t*t*t/6.0;}
00091 
00092 // FIRST DERIVATIVES to the uniform cubic C2 basis functions
00093 inline double dB_0(double t) {return 0.5*(1.0-t)*(t-1.0);}
00094 inline double dB_1(double t) {return (1.5*t*t - 2.0*t);}
00095 inline double dB_2(double t) {return (-(1.5*t*t) + t + 0.5);}
00096 inline double dB_3(double t) {return t*t/2.0;}
00097 
00098 // Lagt til, Vegard
00099 // SECOND DERIVATIVES to the uniform cubic C2 basis functions
00100 inline double ddB_0(double t) {return (-t+1.0);}
00101 inline double ddB_1(double t) {return (3.0*t-2.0);}
00102 inline double ddB_2(double t) {return (-3.0*t+1.0);}
00103 inline double ddB_3(double t) {return t;}
00104 
00105 // PRE EVALUATTED tensor products of B-splines, in a grid index
00106 // -------------------------------------------------------------
00107 
00108 // Pre-evaluated tensors equivalent to B(k,s)*B(l,t)
00109 static double tensor_BB[3][3] = {
00110   {1./36., 1./9., 1./36.},
00111   {1./9., 4./9., 1./9.},
00112   {1./36., 1./9., 1./36.}
00113 };
00114 
00115 
00116 // Pre-evaluated tensor products of B-splines, in a grid index
00117 // Mixed derivative an positional tensor for the x-component.
00118 // equivalent to dB(k,s)*B(l,t).
00119 // (The second row with zeros is not used in evaluations.)
00120 static double tensor_dBB[3][3] = {
00121   {-1./12., -1./3., -1./12.},
00122   {     0.,     0.,      0.},
00123   { 1./12.,  1./3.,  1./12.}
00124 };
00125 
00126 // The transpose of the matrix above, for the y-components,
00127 // equivalent to B(k,s)*dB(l,t)
00128 static double tensor_BdB[3][3] = {
00129   {-1./12., 0., 1./12.},
00130   { -1./3., 0., 1./3.},
00131   {-1./12., 0., 1./12.}
00132 };
00133 
00134 #endif
00135 
00136 
00137 // Generic interface for evaluation to all univariate basis functions
00138 inline double B(int k, double t) {
00139   if (k == 0)
00140     return B_0(t);
00141   else if (k == 1)
00142     return B_1(t);
00143   else if (k == 2)
00144     return B_2(t);
00145   else
00146     return B_3(t);
00147 }
00148 
00149 
00150 // Generic interface to first derivatives
00151 inline double dB(int k, double t) {
00152   if (k == 0)
00153     return dB_0(t);
00154   else if (k == 1)
00155     return dB_1(t);
00156   else if (k == 2)
00157     return dB_2(t);
00158   else
00159     return dB_3(t);
00160 }
00161 
00162 // Lagt til, Vegard
00163 // Generic interface to second derivatives
00164 inline double ddB(int k, double t) {
00165   if (k == 0)
00166     return ddB_0(t);
00167   else if (k == 1)
00168     return ddB_1(t);
00169   else if (k == 2)
00170     return ddB_2(t);
00171   else
00172     return ddB_3(t);
00173 }
00174 
00175 // Evaluation of a bivariate tensor product basis function
00176 inline double w(int k, int l, double s, double t) {
00177   return B(k,s)*B(l,t);
00178 }
00179 
00180 // s,t \in [0,1) (but special on gridlines m and n)
00181 // i,j \in [-1, ???
00182 inline void ijst(int m, int n, double uc, double vc, int& i, int& j, double& s, double& t) {
00183         //int i = std::min((int)uc - 1, m-2);
00184         //int j = std::min((int)vc - 1, n-2);
00185 
00186 #ifdef  MBA_CUBIC_C1
00187         i = 2*((int)uc) - 1;
00188         j = 2*((int)vc) - 1;
00189 #else
00190         i = (int)uc - 1;
00191         j = (int)vc - 1;
00192 #endif
00193 
00194 #ifdef WIN32
00195   s = uc - floor(uc);
00196         t = vc - floor(vc);
00197 #else
00198   s = uc - std::floor(uc);
00199   t = vc - std::floor(vc);
00200 #endif
00201         
00202         // adjust for x or y on gridlines m and n (since impl. has 0 <= x <= m and 0 <= y <= n 
00203 #ifdef  MBA_CUBIC_C1
00204         if (i == 2*m-1) {
00205                 i-=2;
00206                 s = 1;
00207         }
00208         if (j == 2*n-1) {
00209                 j-=2;
00210                 t = 1;
00211         }
00212 #else
00213         if (i == m-1) {
00214                 i--;
00215                 s = 1;
00216         }
00217         if (j == n-1) {
00218                 j--;
00219                 t = 1;
00220         }
00221 #endif
00222 }
00223 
00224 inline void WKLandSum2(double s, double t, double w_kl[4][4], double& sum_w_ab2) {
00225   int k,l;
00226   sum_w_ab2 = 0.0; 
00227   for (k = 0; k <= 3; k++) {
00228     for (l = 0; l <=3; l++) {
00229       double tmp = w(k, l, s, t);
00230       w_kl[k][l] = tmp;
00231       sum_w_ab2 += (tmp*tmp); 
00232     }
00233   }
00234 }
00235 
00236 // for check, should give unity
00237 // inline double sumWKL(double s, double t) {
00238 //  int k,l;
00239 //  double sum_w_ab = 0.0; 
00240 //  for (k = 0; k <= 3; k++) {
00241 //    for (l = 0; l <=3; l++) {
00242 //      sum_w_ab += w(k, l, s, t);
00243 //    }
00244 //  }
00245 //  return sum_w_ab;
00246 //}
00247 
00248 // -------------------------------------------------------------------------
00249 // Refinement (The Oslo algorithm)
00250 // Generic implementation for functional and parametric case
00251 // There are 4 different configurations of vertices in the refined lattice:
00252 
00253 // For the uniform bicubic C2 case
00254 
00255 // 1) On an existing vertex
00256 template <class Type>
00257 inline Type phi_2i_2j(GenMatrix<Type>& mat, int i, int j) {
00258   return 1.0/64.0*( mat(i-1,j-1) + mat(i-1,j+1) + mat(i+1,j-1) + mat(i+1,j+1)
00259     +6.0*(mat(i-1,j) + mat(i,j-1) + mat(i,j+1) + mat(i+1,j)) + 36.0*mat(i,j) );
00260 }
00261 
00262 // 2) On an existing vertical line
00263 template <class Type>
00264 inline Type phi_2i_2jPluss1(GenMatrix<Type>& mat, int i, int j) {
00265   return 1.0/16.0*( mat(i-1,j) + mat(i-1,j+1) + mat(i+1,j) + mat(i+1,j+1)
00266     +6.0*(mat(i,j) + mat(i,j+1)) );
00267 }
00268 
00269 // 3) On an existing horizontal line
00270 template <class Type>
00271 inline Type phi_2iPlus1_2j(GenMatrix<Type>& mat, int i, int j) {
00272   return 1.0/16.0*( mat(i,j-1) + mat(i,j+1) + mat(i+1,j-1) + mat(i+1,j+1)
00273     +6.0*(mat(i,j) + mat(i+1,j)) );
00274 }
00275 
00276 // 4) a new crossing in both directions
00277 template <class Type>
00278 inline Type phi_2iPlus1_2jPlus1(GenMatrix<Type>& mat, int i, int j) {
00279   return 1.0/4.0*( mat(i,j) + mat(i,j+1) + mat(i+1,j) + mat(i+1,j+1) );
00280 }
00281 
00282 
00283 // The Oslo Algorithm For the uniform bicubic C1 case 
00284 // There are 16 configurations:
00285 // ----------------------------
00286 //
00287 //      |       :
00288 //      |       :
00289 //      |       :
00290 //.....9|13...15:4...........
00291 //     7|11    3:16
00292 //      |       :
00293 //      |       :
00294 //_____5|2____12:14__________
00295 //     1|6     8:10
00296 //      |       :
00297 //      |       :
00298 //      |       :
00299 
00300 // !!!! int i, int j denotes (2i) and (2j) in the coefficient matrix PHI !!!
00301 // 1) SW of an existing grid point
00302 template <class Type>
00303 inline Type phi_4iMin1_4jMin1(GenMatrix<Type>& mat, int i, int j) {
00304   return 1.0/16.0*( 9.0*mat(i-1,j-1) + 3.0*(mat(i-1,j) + mat(i,j-1)) + mat(i,j) );
00305 }
00306 
00307 // 2) NE of an existing grid point
00308 template <class Type>
00309 inline Type phi_4i_4j(GenMatrix<Type>& mat, int i, int j) {
00310   return 1.0/16.0*( mat(i-1,j-1) + 3.0*(mat(i-1,j) + mat(i,j-1)) + 9.0*mat(i,j) );
00311 }
00312 
00313 // 3) SW of a new grid point
00314 template <class Type>
00315 inline Type phi_4iPlus1_4jPlus1(GenMatrix<Type>& mat, int i, int j) {
00316   return 1.0/64.0*( mat(i-1,j-1) + 5.0*(mat(i-1,j) + mat(i,j-1)) + 25.0*mat(i,j) + 2.0*(mat(i-1,j+1) + mat(i+1,j-1)) + 10.0*(mat(i+1,j) + mat(i,j+1)) + 4.0*mat(i+1,j+1) );
00317 }
00318 
00319 // 4) NE of a new grid point
00320 template <class Type>
00321 inline Type phi_4iPlus2_4jPlus2(GenMatrix<Type>& mat, int i, int j) {
00322   return 1.0/64.0*( 4.0*mat(i,j) + 5.0*mat(i+2,j+1) + 2*(mat(i+2,j) + mat(i,j+2)) + mat(i+2,j+2) + 25*mat(i+1,j+1) + 10.0*(mat(i+1,j) + mat(i,j+1)) + 5.0*mat(i+1,j+2) );
00323 }
00324 
00325 // 5) NW of an existing grid point
00326 template <class Type>
00327 inline Type phi_4iMin1_4j(GenMatrix<Type>& mat, int i, int j) {
00328   return 1.0/16.0*( 3.0*(mat(i-1,j-1) + mat(i,j)) + 9.0*mat(i-1,j) + mat(i,j-1) );
00329 }
00330 
00331 // 6) SE of an existing grid point
00332 template <class Type>
00333 inline Type phi_4i_4jMin1(GenMatrix<Type>& mat, int i, int j) {
00334   return 1.0/16.0*( 3.0*(mat(i-1,j-1) + mat(i,j)) + mat(i-1,j) + 9.0*mat(i,j-1) );
00335 }
00336 
00337 // 7) SW of a new grid point on an existing vertical line
00338 template <class Type>
00339 inline Type phi_4iMin1_4jPlus1(GenMatrix<Type>& mat, int i, int j) {
00340   return 1.0/32.0*( 3.0*mat(i-1,j-1) + 15.0*mat(i-1,j) + mat(i,j-1) + 5.0*mat(i,j) + 2.0*mat(i,j+1) + 6.0*mat(i-1,j+1) );
00341 }
00342 
00343 // 8) SW of a new grid point on an existing horizontal line
00344 template <class Type>
00345 inline Type phi_4iPlus1_4jMin1(GenMatrix<Type>& mat, int i, int j) {
00346   return 1.0/32.0*( 3.0*mat(i-1,j-1) + mat(i-1,j) + 15*mat(i,j-1) + 5.0*mat(i,j) + 2.0*mat(i+1,j) + 6.0*mat(i+1,j-1) );
00347 }
00348 
00349 // 9) NW of a new grid point on an existing vertical line
00350 template <class Type>
00351 inline Type phi_4iMin1_4jPlus2(GenMatrix<Type>& mat, int i, int j) {
00352   return 1.0/32.0*( 6.0*mat(i-1,j) + 2.0*mat(i,j) + 15*mat(i-1,j+1) + 5.0*mat(i,j+1) + mat(i,j+2) + 3.0*mat(i-1,j+2) );
00353 }
00354 
00355 // 10) SE of a new grid point on an existing horizontal line
00356 template <class Type>
00357 inline Type phi_4iPlus2_4jMin1(GenMatrix<Type>& mat, int i, int j) {
00358   return 1.0/32.0*( 6.0*mat(i,j-1) + 2.0*mat(i,j) + mat(i+2,j) + 15.0*mat(i+1,j-1) + 5.0*mat(i+1,j) + 3.0*mat(i+2,j-1) );
00359 }
00360 
00361 // 11) SE of a new grid point on an existing vertical line
00362 template <class Type>
00363 inline Type phi_4i_4jPlus1(GenMatrix<Type>& mat, int i, int j) {
00364   return 1.0/32.0*( mat(i-1,j-1) + 5.0*mat(i-1,j) + 3.0*mat(i,j-1) + 15.0*mat(i,j) + 2.0*mat(i-1,j+1) + 6.0*mat(i,j+1) );
00365 }
00366 
00367 // 12) NW of a new grid point on an existing horizontal line
00368 template <class Type>
00369 inline Type phi_4iPlus1_4j(GenMatrix<Type>& mat, int i, int j) {
00370   return 1.0/32.0*( mat(i-1,j-1) + 3.0*mat(i-1,j) + 5.0*mat(i,j-1) + 15.0*mat(i,j) + 2.0*mat(i+1,j-1) + 6.0*mat(i+1,j) );
00371 }
00372 
00373 // 13) NE of a new grid point on an existing vertical line
00374 template <class Type>
00375 inline Type phi_4i_4jPlus2(GenMatrix<Type>& mat, int i, int j) {
00376   return 1.0/32.0*( 2.0*mat(i-1,j) + 6.0*mat(i,j) + 5.0*mat(i-1,j+1) + 15.0*mat(i,j+1) + mat(i-1,j+2) + 3.0*mat(i,j+2) );
00377 }
00378 
00379 // 14) NE of a new grid point on an existing horizontal line
00380 template <class Type>
00381 inline Type phi_4iPlus2_4j(GenMatrix<Type>& mat, int i, int j) {
00382   return 1.0/32.0*( 2.0*mat(i,j-1) + 6.0*mat(i,j) + 5.0*mat(i+1,j-1) + 15.0*mat(i+1,j) + mat(i+2,j-1) + 3.0*mat(i+2,j) );
00383 }
00384 
00385 // 15) NW of a new grid point
00386 template <class Type>
00387 inline Type phi_4iPlus1_4jPlus2(GenMatrix<Type>& mat, int i, int j) {
00388   return 1.0/64.0*( 2.0*(mat(i-1,j) + mat(i+1,j+2)) + 10.0*(mat(i,j) + mat(i+1,j+1)) + 5.0*(mat(i-1,j+1) + mat(i,j+2)) + 25.0*mat(i,j+1) + 4.0*mat(i+1,j) + mat(i-1,j+2) );
00389 }
00390 
00391 // 16) SE of a new grid point
00392 template <class Type>
00393 inline Type phi_4iPlus2_4jPlus1(GenMatrix<Type>& mat, int i, int j) {
00394   return 1.0/64.0*( 2.0*(mat(i,j-1) + mat(i+2,j+1)) + 10.0*(mat(i,j) + mat(i+1,j+1)) + 5.0*(mat(i+1,j-1) + mat(i+2,j)) + 4.0*mat(i,j+1) + 25.0*mat(i+1,j) + mat(i+2,j-1) );
00395 }
00396 
00397 #endif // _UNIFBSPLINES_H_

Generated on Wed Nov 28 12:27:19 2007 for MBA by  doxygen 1.5.1