Questo sito utilizza cookies, anche di terze parti, per mostrare pubblicità e servizi in linea con il tuo account. Leggi l'informativa sui cookies.
Username: Password: oppure
C/C++ - spline perche' non funziona?
Forum - C/C++ - spline perche' non funziona?

Avatar
emarfi (Normal User)
Newbie


Messaggi: 3
Iscritto: 19/02/2010

Segnala al moderatore
Postato alle 2:03
Venerdì, 19/02/2010
Il problema e' interpolare e poi visualizzare (su un graffo), n punti scielti a caso con
una spline cubica.
ecco cosa ho fatto fin ora:
non so' cosa sia il problema, ho guardato e riguardato e riguardato il codice e non so' perche non funzioni.
i punti in excel formano un graffo, ma le tangenti a sinistra dai punti dati non coincidono con quelle alla destra, come devono avere in una spline cubica.

sarei contentissimo anche per qualche idea, grazie

Codice sorgente - presumibilmente C++

  1. #include<iostream>
  2. #include<fstream>
  3. #include<math.h>
  4. using namespace std;
  5. #include<Eigen/Core>
  6. USING_PART_OF_NAMESPACE_EIGEN
  7. #include<Eigen/Geometry>
  8. #include<Eigen/LU>
  9. #include"Punto.h"
  10. #include"Spline2D.h"
  11.  
  12. Spline::Spline()
  13. {
  14.         n = 10;//numero di punti
  15.         pt = new Punto[n];
  16. }
  17. Spline::~Spline()
  18. {
  19.         delete [] pt;
  20. }
  21. void Spline::GeneraPunti(void)
  22. {
  23.         cout<<"Random punti generati sono:\n";
  24.         for(int i=0; i<n; i++)
  25.         {
  26.                 int bigg = int(double(rand()% 10 + (i*10)));
  27.                 int rndX = bigg;
  28.                 int rndY = int((double(rand())/RAND_MAX)*9);
  29.                 pt[i].SetPunto(rndX,rndY);
  30.                 cout<<"T"<<i+1<<": ("<<pt[i].x<<", "<<pt[i].y<<")"<<endl;
  31.         }
  32.        
  33.         /*int rndX[6] = {1,3,5,7,9,13};
  34.         int rndY[6] = {2,8,4,6,-2,3};
  35.         for(int i=0; i<n; i++){
  36.                 pt[i].SetPunto(rndX[i],rndY[i]);
  37.         }*/
  38. }
  39. void Spline::SplineInterp(void)
  40. {
  41.        
  42.         MatrixXd A;
  43.         A.setZero(n-2,n-2);
  44.         VectorXd Y;
  45.         Y.setZero(n-2);
  46.        
  47.         for(int i=0; i<n-3; i++)
  48.         {
  49.                 A(i,i)=4;
  50.                 A(i,i+1)=1;
  51.                 A(i+1,i)=1;
  52.                 Y(i) = 6*(pt[i].y- 2*pt[i+1].y+pt[i+2].y)/(pow((pt[i+1].x-pt[i].x),2));
  53.         }
  54.         Y(n-3) = 6*(pt[n-3].y- 2*pt[n-2].y+pt[n-1].y)/(pow((pt[n-2].x-pt[n-3].x),2));
  55.         A(n-3,n-3) = 4;
  56.        
  57.         cout<<" A\n"<<A<<endl;
  58.         cout<<" Y\n"<<Y<<endl;
  59.         VectorXd MM;
  60.         MM.setZero(n-2);
  61.         A.lu().solve(Y,&MM);
  62.        
  63.         VectorXd M;
  64.         M.setZero(n);
  65.         M[0]=0;
  66.         for(int i=0; i<n-2; i++)
  67.         {
  68.                 M[i+1]=MM[i];
  69.         }
  70.         M[n-1]=0;
  71.        
  72.        
  73.         ////////////////////////////////////////////coefficienti a, b, c, d /////////////////////////////////////////////
  74.         VectorXd a;       VectorXd b;               VectorXd c;         VectorXd d;
  75.         a.setZero(n-1);   b.setZero(n-1);       c.setZero(n-1);         d.setZero(n-1);
  76.         S.setZero(1000);
  77.         X.setZero(1000);
  78.  
  79.         for(int i=0; i<n-1; i++)////(da i=1,...,n-1.)
  80.         {
  81.                 a(i) = (M(i+1)-M(i))/(6.*(pt[i+1].x-pt[i].x));
  82.                 b(i) = M(i)/2;
  83.                 c(i) = (pt[i+1].y-pt[i].y)/(pt[i+1].x-pt[i].x) - (M(i+1)+2*M(i))*(pt[i+1].x-pt[i].x)/6.;
  84.                 d(i) = pt[i].y;
  85.         }
  86.                
  87.         cout<<"\ta\n"<<a<<"\n\tb\n "<<b<<"\n\tc\n "<<c<<"\n\td\n "<<d<<endl;
  88.  
  89.         br=0;//counter di tutti i punti
  90.  
  91.         for(int i=0; i<n-1; i++)
  92.         {      
  93.                 br++;//io primo punto
  94.                 X[br] = pt[i].x;
  95.                 S[br] = pt[i].y;
  96.                 cout<<"X="<<X[br]<<" S="<<S[br]<<endl;
  97.                 double step = 1;
  98.                 if (pt[i].x < pt[i+1].x)
  99.                 {      
  100.                         for(double x = (pt[i].x+step); x < pt[i+1].x; x+=step)
  101.                         {                      
  102.                                 br++;
  103.                                 X[br] = x;
  104.                                 S[br] = a(i)*pow((x-pt[i].x), 3) + b(i)*pow((x-pt[i].x), 2) +
  105.                                                 c(i)*(x-pt[i].x) + d(i);
  106.                                 cout<<"X="<<X[br]<<" S="<<S[br]<<endl;
  107.                         }
  108.                         continue;
  109.                 }
  110.                
  111.         }
  112.         br++;//l'ultimo punto
  113.         X[br] = pt[n-1].x;
  114.         S[br] = pt[n-1].y;
  115.         cout<<"counter : "<<br<<" "<<endl;
  116.  
  117.  
  118.  
  119.         for(int i=0; i<n; i++)
  120.         {
  121.                 cout<<"T"<<i+1<<": ("<<pt[i].x<<", "<<pt[i].y<<")"<<endl;
  122.         }
  123. }



ecco i punti:
1     2    
2     6,26136    
3     8    
4     5,96591    
5     4    
6     5,375    
7     6    
8     2,03409    
9     -2    
10     -6,47727    
11     -6,04545    
12     -2,34091    
13     3    

il file allegato e' il grafo della spline

EDIT by HeDo: Il codice va racchiuso tra tag code, oltretutto credo sia "grafo" e non "graffo"... vabbe essere sardi mah...

Ultima modifica effettuata da emarfi il 20/02/2010 alle 0:35
PM Quote
Avatar
TheKaneB (Member)
Guru^2


Messaggi: 1787
Iscritto: 26/06/2009

Segnala al moderatore
Postato alle 9:58
Venerdì, 19/02/2010
scusa... ma che razza di linguaggio è?

Codice sorgente - presumibilmente Plain Text

  1. a(i) = (M(i+1)-M(i))/(6.*(pt[i+1].x-pt[i].x));



mi puzza un po' di matlab....
fammi indovinare: copia-incolla allo stato brado?


Software Failure: Guru Meditation
Forum su Informatica, Elettronica, Robotica e Tecnologia: http://www.nonsoloamiga.com
PM Quote
Avatar
emarfi (Normal User)
Newbie


Messaggi: 3
Iscritto: 19/02/2010

Segnala al moderatore
Postato alle 11:16
Venerdì, 19/02/2010
no e' c++, (visual studio 2008) ma con l'uso delle librerie eigen, si possono trovare con google. eigen permette di "operare" in modo semplice con le matrici, vettori...
il codice funziona, ma l'algoritmo non va bene. e' l'implementazione di, vedi qua':http://online.redwoods.cc.ca.us/instruct/darnold/laproj/Fa ...
ilcodice sopra non e' tutto il codice del programma, e' il codice della classe Spline, che poi si chiama nel main.

Ultima modifica effettuata da emarfi il 19/02/2010 alle 11:44
PM Quote
Avatar
emarfi (Normal User)
Newbie


Messaggi: 3
Iscritto: 19/02/2010

Segnala al moderatore
Postato alle 17:21
Venerdì, 19/02/2010
Ho trovato il problema, l' algoritmo funziona per una spline uniforme, cioe' di intervallo (x(i+1)-x(i))costante.
Per una spline non-uniforme ho dovtuo cambiare un po' le matrici.

Ultima modifica effettuata da emarfi il 19/02/2010 alle 17:22
PM Quote
Avatar
diegodiego (Normal User)
Newbie


Messaggi: 1
Iscritto: 13/06/2011

Segnala al moderatore
Postato alle 15:07
Lunedì, 13/06/2011
potresti postare il codice corretto?
non ho capito che cosa intendi per spline non uniforme:)

PM Quote