gunu_global_interpolate.cpp Source File

00001 /* LIBGUL - Geometry Utility Library
00002  * Copyright (C) 1998-1999 Norbert Irmer
00003  *
00004  * This library is free software; you can redistribute it and/or
00005  * modify it under the terms of the GNU Library General Public
00006  * License as published by the Free Software Foundation; either
00007  * version 2 of the License, or (at your option) any later version.
00008  *
00009  * This library is distributed in the hope that it will be useful,
00010  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012  * Library General Public License for more details.
00013  *
00014  * You should have received a copy of the GNU Library General Public
00015  * License along with this library; if not, write to the
00016  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
00017  * Boston, MA 02111-1307, USA.
00018  */
00019 
00020 #include "stdafx.h"
00021 
00022 #include <iostream>
00023 
00024 #include "gul_types.h"
00025 #include "gunu_parametrize.h"
00026 #include "guma_lineq.h"
00027 #include "gunu_basics.h"
00028 #include "gul_io.h"
00029 
00030 namespace gunu {
00031 
00032 using gul::Ptr;
00033 using gul::point;
00034 using gul::point1;
00035 using gul::point2;
00036 using guma::BandDecomposition;
00037 using guma::BandBackSubstitution;
00038 using gul::dump_matrix;
00039 using gul::dump_surf;
00040 
00041 template< class T >
00042 void CIDecompose( int nQ, const Ptr<T>& Uq, int p, const Ptr<T>& U, 
00043                   Ptr< Ptr<T> >& A, Ptr< Ptr<T> >& L, Ptr<int>& index )
00044 {
00045   Ptr<T> B;
00046   int i,j,span,offset,n,sign;
00047   T u;
00048 
00049   n = nQ-1;
00050 
00051   A[0][p-1] = (T)1;
00052   for( i = p; i < p+p-1; i++ ) A[0][i] = (T)0; 
00053   A[n][p-1] = (T)1;
00054   for( i = 0; i < p-1; i++ ) A[n][i] = (T)0; 
00055 
00056   span = p;
00057   offset = 0;
00058 
00059   for( i = 1; i < n; i++ )
00060   {
00061     u = Uq[i];
00062 
00063     if( (u >= U[span+1]) && (span != n) )
00064       span++;
00065     else
00066       offset--;  // if in the same knotspan as the last point 
00067 
00068     /* cout << "offset = " << offset << "\n"; */
00069 
00070     for( j = 0; j < p-1+offset; j++ ) A[i][j] = (T)0;
00071     for( j = p+p+offset; j < p+p-1; j++ ) A[i][j] = (T)0;
00072     B.use_pointer( &A[i][p-1+offset], p+1 );
00073     CalcBasisFunctions( u, span, p, U, B );
00074   }
00075   BandDecomposition( n+1, p-1, p-1, A, L, &sign, index );
00076 }
00077 
00078 template< class T >
00079 void CIBackSubst(
00080            int nQ, const Ptr< point<T> >& Q, int p, Ptr< Ptr<T> >& A,
00081            const Ptr< Ptr<T> >& L, const Ptr<int>& index, Ptr< point<T> >& P )
00082 {
00083   Ptr<T> b;
00084   int n = nQ-1, i;
00085 
00086   b.reserve_pool(n+1);
00087 
00088   // calculate the x coordinates of the controlpoints
00089   for( i = 0; i <= n; i++ ) b[i] = Q[i].x;
00090   BandBackSubstitution( n+1, p-1, p-1, A, L, index, b ); 
00091   for( i = 0; i <= n; i++ ) P[i].x = b[i];
00092 
00093   // calculate the y coordinates of the controlpoints
00094   for( i = 0; i <= n; i++ ) b[i] = Q[i].y;
00095   BandBackSubstitution( n+1, p-1, p-1, A, L, index, b ); 
00096   for( i = 0; i <= n; i++ ) P[i].y = b[i];
00097 
00098   // calculate the z coordinates of the controlpoints
00099   for( i = 0; i <= n; i++ ) b[i] = Q[i].z;
00100   BandBackSubstitution( n+1, p-1, p-1, A, L, index, b ); 
00101   for( i = 0; i <= n; i++ ) P[i].z = b[i];
00102 }
00103 
00104 template< class T >
00105 void CIBackSubst(
00106            int nQ, const Ptr< point2<T> >& Q, int p, const Ptr< Ptr<T> >& A,
00107            const Ptr< Ptr<T> >& L, const Ptr<int>& index, Ptr< point2<T> >& P )
00108 {
00109   Ptr<T> b;
00110   int n = nQ-1, i;
00111 
00112   b.reserve_pool(n+1);
00113 
00114   // calculate the x coordinates of the controlpoints
00115   for( i = 0; i <= n; i++ ) b[i] = Q[i].x;
00116   BandBackSubstitution( n+1, p-1, p-1, A, L, index, b ); 
00117   for( i = 0; i <= n; i++ ) P[i].x = b[i];
00118 
00119   // calculate the y coordinates of the controlpoints
00120   for( i = 0; i <= n; i++ ) b[i] = Q[i].y;
00121   BandBackSubstitution( n+1, p-1, p-1, A, L, index, b ); 
00122   for( i = 0; i <= n; i++ ) P[i].y = b[i];
00123 }
00124 
00125 template< class T >
00126 void CIBackSubstToCol(
00127            int nQ, const Ptr< point<T> >& Q, int p, const Ptr< Ptr<T> >& A,
00128            const Ptr< Ptr<T> >& L, const Ptr<int>& index, 
00129            int iCol, Ptr< Ptr< point<T> > >& P )
00130 {
00131   Ptr<T> b;
00132   int n = nQ-1, i;
00133 
00134   b.reserve_pool(n+1);
00135 
00136   // calculate the x coordinates of the controlpoints
00137   for( i = 0; i <= n; i++ ) b[i] = Q[i].x;
00138   BandBackSubstitution( n+1, p-1, p-1, A, L, index, b ); 
00139   for( i = 0; i <= n; i++ ) P[i][iCol].x = b[i];
00140 
00141   // calculate the y coordinates of the controlpoints
00142   for( i = 0; i <= n; i++ ) b[i] = Q[i].y;
00143   BandBackSubstitution( n+1, p-1, p-1, A, L, index, b ); 
00144   for( i = 0; i <= n; i++ ) P[i][iCol].y = b[i];
00145 
00146   // calculate the z coordinates of the controlpoints
00147   for( i = 0; i <= n; i++ ) b[i] = Q[i].z;
00148   BandBackSubstitution( n+1, p-1, p-1, A, L, index, b ); 
00149   for( i = 0; i <= n; i++ ) P[i][iCol].z = b[i];
00150 }
00151 
00152 template< class T >
00153 void CIBackSubstToCol(
00154            int nQ, const Ptr< point1<T> >& Q, int p, const Ptr< Ptr<T> >& A,
00155            const Ptr< Ptr<T> >& L, const Ptr<int>& index, 
00156            int iCol, Ptr< Ptr< point1<T> > >& P )
00157 {
00158   Ptr<T> b;
00159   int n = nQ-1, i;
00160 
00161   b.reserve_pool(n+1);
00162 
00163   // calculate the x coordinates of the controlpoints
00164   for( i = 0; i <= n; i++ ) b[i] = Q[i].x;
00165   BandBackSubstitution( n+1, p-1, p-1, A, L, index, b ); 
00166   for( i = 0; i <= n; i++ ) P[i][iCol].x = b[i];
00167 }
00168 
00169 template< class T, class EP >
00170 GULAPI bool GlobalCurveInterpolation( 
00171                                int nQ, const Ptr<EP>& Q, int p,
00172                                Ptr<T>& U, Ptr<EP>& P )                              
00173 {
00174   Ptr<T> Uq,d;
00175   Ptr< Ptr<T> > A,L;  
00176   Ptr<int> index;
00177   T dsum;
00178   int i;
00179 
00180   d.reserve_pool(nQ);
00181   Uq.reserve_pool(nQ);
00182   index.reserve_pool(nQ);
00183 
00184   dsum = ChordLength( nQ, Q, d );
00185   if( dsum == (T)0 ) return false;
00186 
00187   ChordLengthParameters( nQ, dsum, d, Uq );
00188   KnotVectorByAveraging( nQ, Uq, p, U );
00189 
00190   A.reserve_pool(nQ);
00191   L.reserve_pool(nQ);
00192   for( i = 0; i < nQ; i++ )
00193   {
00194     A[i].reserve_pool(p+p-1);
00195     L[i].reserve_pool(p+p-1);
00196   }
00197   index.reserve_pool(nQ);
00198 
00199   CIDecompose(nQ, Uq, p, U, A, L, index);
00200   CIBackSubst<T>(nQ, Q, p, A, L, index, P);
00201 
00202   return true;
00203 }
00204 // template instantiation
00205 template GULAPI bool GlobalCurveInterpolation( 
00206                   int nQ, const Ptr< point<float> >& Q,
00207                   int p, Ptr<float>& U, Ptr< point<float> >& P );
00208 template GULAPI bool GlobalCurveInterpolation( 
00209                   int nQ, const Ptr< point<double> >& Q,
00210                   int p, Ptr<double>& U, Ptr< point<double> >& P );
00211 
00212 template GULAPI bool GlobalCurveInterpolation( 
00213                           int nQ, const Ptr< point2<float> >& Q,
00214                           int p, Ptr<float>& U, Ptr< point2<float> >& P );
00215 template GULAPI bool GlobalCurveInterpolation( 
00216                           int nQ, const Ptr< point2<double> >& Q,
00217                           int p, Ptr<double>& U, Ptr< point2<double> >& P );
00218 
00219 template< class T, class EP >
00220 GULAPI bool GlobalSurfaceInterpolation( 
00221            int nRows, int nCols, const Ptr< Ptr<EP> >& Q,
00222            int pu, int pv, Ptr<T>& U, Ptr<T>& V, Ptr< Ptr<EP> > &P )
00223 {
00224   Ptr<T> Uq, Vq;
00225   int nu,nv,i;
00226   Ptr< Ptr<T> > A,L;  
00227   Ptr<int> index;
00228   Ptr< Ptr<EP> > R;
00229   nu = nCols-1;
00230   nv = nRows-1;
00231 
00232   Uq.reserve_pool(nCols);
00233   Vq.reserve_pool(nRows);
00234 
00235   if( !ChordLengthMeshParams( nRows, nCols, Q, Uq, Vq ) ) return false; 
00236 
00237   KnotVectorByAveraging( nu+1, Uq, pu, U );
00238   KnotVectorByAveraging( nv+1, Vq, pv, V );
00239 
00240   R.reserve_pool(nu+1);
00241   for( i = 0; i <= nu; i++ )
00242     R[i].reserve_pool(nv+1);
00243 
00244   A.reserve_pool(nu+1);
00245   L.reserve_pool(nu+1);
00246   for( i = 0; i <= nu; i++ )
00247   {
00248     A[i].reserve_pool(pu+pu-1);
00249     L[i].reserve_pool(pu+pu-1);
00250   }
00251   index.reserve_pool(nu+1);
00252  
00253   CIDecompose(nu+1, Uq, pu, U, A, L, index);
00254 
00255   for( i = 0; i <= nv; i++ )
00256   {
00257     CIBackSubstToCol(nu+1, Q[i], pu, A, L, index, i, R);
00258   }
00259 
00260   A.reserve_pool(nv+1);
00261   L.reserve_pool(nv+1);
00262   for( i = 0; i <= nv; i++ )
00263   {
00264     A[i].reserve_pool(pv+pv-1);
00265     L[i].reserve_pool(pv+pv-1);
00266   }
00267   index.reserve_pool(nv+1);
00268 
00269   CIDecompose(nv+1, Vq, pv, V, A, L, index);
00270   for( i = 0; i <= nu; i++ )
00271   {
00272     CIBackSubstToCol(nv+1, R[i], pv, A, L, index, i, P);
00273   }
00274 
00275   return true;
00276 }
00277 // template instantiation
00278 template GULAPI bool GlobalSurfaceInterpolation< float,point<float> >( 
00279            int nRows, int nCols, const Ptr< Ptr< point<float> > >& Q,
00280            int pu, int pv, Ptr<float>& U, Ptr<float>& V,
00281            Ptr< Ptr< point<float> > >& P );
00282 template GULAPI bool GlobalSurfaceInterpolation< double,point<double> >( 
00283            int nRows, int nCols, const Ptr< Ptr< point<double> > >& Q,
00284            int pu, int pv, Ptr<double>& U, Ptr<double>& V,
00285            Ptr< Ptr< point<double> > >& P );
00286 
00287 /*
00288 template bool GlobalSurfaceInterpolation< float,point1<float> >( 
00289            int nRows, int nCols, const Ptr< Ptr< point1<float> > >& Q,
00290            int pu, int pv, Ptr<float>& U, Ptr<float>& V,
00291            Ptr< Ptr< point1<float> > >& P );
00292 template bool GlobalSurfaceInterpolation< double,point1<double> >( 
00293            int nRows, int nCols, const Ptr< Ptr< point1<double> > >& Q,
00294            int pu, int pv, Ptr<double>& U, Ptr<double>& V,
00295            Ptr< Ptr< point1<double> > >& P );
00296 */
00297 }
00298