00001 /******************************************************************************
00002 
00003  Copyright 2008 Departamento de Realidad Virtual
00004  y Unidad de Cómputo Aplicado DGSGA, UNAM.
00005 
00006 
00007  This file is part of RBF++.
00008 
00009  RBF++ is free software: you can redistribute it and/or modify
00010  it under the terms of the GNU General Public License as published by
00011  the Free Software Foundation, either version 3 of the License, or
00012  (at your option) any later version.
00013 
00014  RBF++ is distributed in the hope that it will be useful,
00015  but WITHOUT ANY WARRANTY; without even the implied warranty of
00016  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00017  GNU General Public License for more details.
00018 
00019  You should have received a copy of the GNU General Public License
00020  along with RBF++. If not, see <http://www.gnu.org/licenses/>.
00021 
00022 
00023 *******************************************************************************/
00024 
00025 /* 
00026  Author: Daniel Cervantes Cabrera.
00027  Supervision: Pedro Gonzalez Casanova.
00028  Project: RBF++
00029  Institution: DGSCA, UNAM.
00030  Date: 3/05/08
00031  Description: Thin spline interpolation example.
00032  */
00033 
00034 
00035 #include <LA/Point2dT.h>
00036 #include <LA/Point3dT.h>
00037 #include <LA/Matrix.h>
00038 #include <RBF/GrammMatrix.h>
00039 #include <RBF/TPSInterp.h>
00040 
00041 #include <vector>
00042 #include <stdlib.h>
00043 #include <time.h>
00044 #include <math.h>
00045 
00046 
00047 float f(float x, float y){
00048   if (0 <= x &&  x <= 1 && 0 <= y && y <= 1)
00049     return 0.75 * exp(-((9*x - 2)*(9*x-2) + (9*x - 2)*(9*x - 2))/4.0) 
00050       + 0.75 * exp(-((9*x - 2)*(9*x-2)/49 - (9*y - 2)*(9*y - 2))/10.0) 
00051       + 0.5 * exp(-((9*x - 7)*(9*x-7) + (9*y - 3)*(9*y - 3))/4.0)
00052       - 0.2 * exp(-(9*x - 4)*(9*x-4) - (9*y - 7)*(9*y - 0)); 
00053   else
00054     return 0;  
00055 }
00056 
00057 
00058 float f2(float x, float y)
00059 {
00060   if (0 <= x &&  x <= 1 && 0 <= y && y <= 1)
00061     return ( 5.0 / 4.0 + cos(5.4 * y) ) / 
00062       ( 6 + 6 * pow(( 3 * x - 1 ), 2) );  
00063   else
00064     return 0;
00065 
00066 }
00067 
00068 int main(int argc, char * argv[]){
00069 
00070   
00071   if(argc < 2){
00072     std::cout << "Usar ./rbfinterp na" << std::endl 
00073               << "na = numeros aleatorios " 
00074               <<  std::endl;
00075     
00076     return 0;
00077   }
00078 
00079 
00080   int N = atoi(argv[1]);
00081 
00082   
00083 
00084   RBF::TPSInterp<double> tpsinterpf;  
00085   std::vector<double> d;
00086   std::vector<LA::Point3dd> points;       
00087   LA::Point2dd p2d;
00088 
00089   points.resize(N);
00090 
00091 
00092   srand(time(0));
00093   for(int i=0; i< N; i++){
00094     points[i].X() =  float(rand())/RAND_MAX;
00095     points[i].Y() =  float(rand())/RAND_MAX;
00096     points[i].Z() = f(points[i].X(),points[i].Y());
00097       
00098   }
00099 
00100   tpsinterpf.build(points);
00101   
00102   N = 50;
00103 
00104   
00105   FILE *  fpt = fopen("interpolation.txt","w");
00106   
00107   for(int i=0; i < N; i++)
00108     for(int j=0; j < N; j++){
00109       p2d[0] = float(i)/N;
00110       p2d[1] = float(j)/N;
00111       fprintf(fpt,"%4.6f %4.6f %4.6f\n",p2d[0],p2d[1],tpsinterpf(p2d));         
00112     }         
00113   
00114   fclose(fpt);
00115 }