Oppure

Loading
19/02/10 1:03
emarfi
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



#include<iostream>
#include<fstream>
#include<math.h>
using namespace std;
#include<Eigen/Core>
USING_PART_OF_NAMESPACE_EIGEN
#include<Eigen/Geometry>
#include<Eigen/LU>
#include"Punto.h"
#include"Spline2D.h"

Spline::Spline()
{
	n = 10;//numero di punti
	pt = new Punto[n];
}
Spline::~Spline()
{
	delete [] pt;
}
void Spline::GeneraPunti(void)
{
	cout<<"Random punti generati sono:\n";
	for(int i=0; i<n; i++)
	{
		int bigg = int(double(rand()% 10 + (i*10)));
		int rndX = bigg;
		int rndY = int((double(rand())/RAND_MAX)*9);
		pt[i].SetPunto(rndX,rndY);
		cout<<"T"<<i+1<<": ("<<pt[i].x<<", "<<pt[i].y<<")"<<endl;
	}
	
	/*int rndX[6] = {1,3,5,7,9,13};
	int rndY[6] = {2,8,4,6,-2,3};
	for(int i=0; i<n; i++){
		pt[i].SetPunto(rndX[i],rndY[i]);
	}*/
}
void Spline::SplineInterp(void)
{
	
	MatrixXd A;
	A.setZero(n-2,n-2);
	VectorXd Y;
	Y.setZero(n-2);
	
	for(int i=0; i<n-3; i++)
	{
		A(i,i)=4;
		A(i,i+1)=1;
		A(i+1,i)=1;
		Y(i) = 6*(pt[i].y- 2*pt[i+1].y+pt[i+2].y)/(pow((pt[i+1].x-pt[i].x),2));
	}
	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));
	A(n-3,n-3) = 4;
	
	cout<<" A\n"<<A<<endl;
	cout<<" Y\n"<<Y<<endl;
	VectorXd MM;
	MM.setZero(n-2);
	A.lu().solve(Y,&MM);
	
	VectorXd M;
	M.setZero(n);
	M[0]=0;
	for(int i=0; i<n-2; i++)
	{
		M[i+1]=MM[i];
	}
	M[n-1]=0;
	
	
	////////////////////////////////////////////coefficienti a, b, c, d /////////////////////////////////////////////
	VectorXd a;       VectorXd b;		    VectorXd c;         VectorXd d;
	a.setZero(n-1);	  b.setZero(n-1);    	c.setZero(n-1); 	d.setZero(n-1);
	S.setZero(1000);
	X.setZero(1000);

	for(int i=0; i<n-1; i++)////(da i=1,...,n-1.)
	{
		a(i) = (M(i+1)-M(i))/(6.*(pt[i+1].x-pt[i].x));
		b(i) = M(i)/2;
		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.;
		d(i) = pt[i].y;
	}
		
	cout<<"\ta\n"<<a<<"\n\tb\n "<<b<<"\n\tc\n "<<c<<"\n\td\n "<<d<<endl;

	br=0;//counter di tutti i punti

	for(int i=0; i<n-1; i++)
	{	
		br++;//io primo punto
		X[br] = pt[i].x;
		S[br] = pt[i].y;
		cout<<"X="<<X[br]<<" S="<<S[br]<<endl;
		double step = 1;
		if (pt[i].x < pt[i+1].x)
		{	
			for(double x = (pt[i].x+step); x < pt[i+1].x; x+=step)
			{			
				br++;
				X[br] = x;
				S[br] = a(i)*pow((x-pt[i].x), 3) + b(i)*pow((x-pt[i].x), 2) + 
						c(i)*(x-pt[i].x) + d(i);
				cout<<"X="<<X[br]<<" S="<<S[br]<<endl;
			}
			continue;
		}
		
	}
	br++;//l'ultimo punto
	X[br] = pt[n-1].x;
	S[br] = pt[n-1].y;
	cout<<"counter : "<<br<<" "<<endl;



	for(int i=0; i<n; i++)
	{
		cout<<"T"<<i+1<<": ("<<pt[i].x<<", "<<pt[i].y<<")"<<endl;
	}
}









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 19/02/10 23:35
aaa
19/02/10 8:58
TheKaneB
scusa... ma che razza di linguaggio è?

 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?
aaa
19/02/10 10:16
emarfi
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':online.redwoods.cc.ca.us/instruct/darnold/laproj/Fall98/SkyMeg/…
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 19/02/10 10:44
aaa
19/02/10 16:21
emarfi
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 19/02/10 16:22
aaa
13/06/11 13:07
diegodiego
potresti postare il codice corretto?
non ho capito che cosa intendi per spline non uniforme:)
aaa