#include "creaVtkHeartAngles.h"
#include <math.h>
#include <iostream>
+#include "vtkDoubleArray.h"
+#include "vtkStructuredPoints.h"
+#include "vtkPointData.h"
+#include "vtkDataArray.h"
creaVtkHeartAngles::creaVtkHeartAngles()
{
double alpha;
//Calcular el ángulo que forman las rectas, sabiendo sus vectores directores.
- double cosAlpha = ( ( (x-a)*(xp-a) )+( (y-b)*(yp-b) )+( (z-c)*(zp-c) ) )/sqrt((( (pow(x,2.0)-pow(a,2.0)) + (pow(y,2.0)-pow(b,2.0)) + (pow(z,2.0)-pow(c,2.0)) )*( (pow(xp,2.0)-pow(a,2.0)) + (pow(yp,2.0)-pow(b,2.0)) + (pow(zp,2.0)-pow(c,2.0)) ) ));
+ //double cosAlpha = ( ( (x-a)*(xp-a) )+( (y-b)*(yp-b) )+( (z-c)*(zp-c) ) )/sqrt((( (pow(x,2.0)-pow(a,2.0)) + (pow(y,2.0)-pow(b,2.0)) + (pow(z,2.0)-pow(c,2.0)) )*( (pow(xp,2.0)-pow(a,2.0)) + (pow(yp,2.0)-pow(b,2.0)) + (pow(zp,2.0)-pow(c,2.0)) ) ));
+ double cosAlpha = ( ( x*xp )+( y*yp )+( z*zp ) )/(sqrt( pow(x,2.0) + pow(y,2.0) + pow(z,2.0) )*sqrt ( pow(xp,2.0) + pow(yp,2.0) + pow(zp,2.0) ) );
alpha = (180 * acos(fabs(cosAlpha)))/M_PI;
- //std::cout << "Angulo Alpha: " << alpha << std::endl;
+ /*if((alpha<=90)&&(alpha>=0)){
+ std::cout << " " << a << " " << b <<" " << c <<" " << x <<" " << y <<" " << z <<" " << xp <<" " << yp <<" " << zp;
+ std::cout << " cosAlpha: " << cosAlpha;
+ std::cout << " acos(fabs(cosalpha)) " << acos(fabs(cosAlpha));
+ std::cout << " Angulo Alpha: " << alpha << std::endl ; }*/
double ent = floor(alpha);
double al = alpha - ent;
double min = al * 60;
//std::cout << "minutos: " << min << std::endl;
- std::cout << "Grados: "<< ent <<" minutos: " << floor(min) << std::endl;
+ //std::cout << "Alpha grados: "<< ent <<" minutos: " << floor(min) << " cos: " << cosAlpha << std::endl;
return alpha;
}
double beta;
//Calcular el ángulo que forman las rectas, sabiendo sus vectores directores.
- double cosBeta = ( ( (x-a)*(xp-a) )+( (y-b)*(yp-b) )+( (z-c)*(zp-c) ) )/sqrt((( (pow(x,2.0)-pow(a,2.0)) + (pow(y,2.0)-pow(b,2.0)) + (pow(z,2.0)-pow(c,2.0)) )*( (pow(xp,2.0)-pow(a,2.0)) + (pow(yp,2.0)-pow(b,2.0)) + (pow(zp,2.0)-pow(c,2.0)) ) ));
-
+ //double cosBeta = ( ( (x-a)*(xp-a) )+( (y-b)*(yp-b) )+( (z-c)*(zp-c) ) )/sqrt((( (pow(x,2.0)-pow(a,2.0)) + (pow(y,2.0)-pow(b,2.0)) + (pow(z,2.0)-pow(c,2.0)) )*( (pow(xp,2.0)-pow(a,2.0)) + (pow(yp,2.0)-pow(b,2.0)) + (pow(zp,2.0)-pow(c,2.0)) ) ));
+ double cosBeta = ( ( x*xp )+( y*yp )+( z*zp ) )/(sqrt( pow(x,2.0) + pow(y,2.0) + pow(z,2.0) )*sqrt ( pow(xp,2.0) + pow(yp,2.0) + pow(zp,2.0) ) );
beta = (180 * acos(fabs(cosBeta)))/M_PI;
//std::cout << "Angulo Beta: " << beta << std::endl;
double ent = floor(beta);
double al = beta - ent;
double min = al * 60;
//std::cout << "minutos: " << min << std::endl;
- std::cout << "Grados: "<< ent <<" minutos: " << floor(min) << std::endl;
+ //std::cout << "Beta grados: "<< ent <<" minutos: " << floor(min) << " cos: " << cosBeta << std::endl;
return beta;
}
-double * creaVtkHeartAngles::vectorProjection (double plX, double plY, double plZ, double nX, double nY, double nZ, double vX, double vY, double vZ)
+double * creaVtkHeartAngles::vectorProjection (double nX, double nY, double nZ, double vX, double vY, double vZ)
{
- //punto del plano
- double plx=plX, ply=plY, plz=plZ;
+
//normalPlano
double nx=nX, ny=nY, nz=nZ;
//vector V
double vx=vX, vy=vY, vz=vZ;
- //distancia del punto al plano
- double dist;
//vector proyectado
- double vpx, vpy, vpz;
double proj [3];
- dist = ply - vy;
- //std::cout << dist << std
- vpx = vx - (dist * nx);
- vpy = vy - (dist * ny);
- vpz = vz - (dist * nz);
-
- proj[0]=vpx;
- proj[1]=vpy;
- proj[2]=vpz;
+ proj[0]= vX - (nX*vX);
+ proj[1]= vY - (nY*vY);
+ proj[2]= vZ - (nZ*vZ);
return proj;
}
return puntoCorte;
}
-void creaVtkHeartAngles::calculateImages (vtkImageData* image, double Px, double Py, double Pz, double Vx, double Vy, double Vz, double Nx, double Ny, double Nz, double p2x, double p2y, double p2z)
+void creaVtkHeartAngles::calculateImages (vtkImageData* image, double Px, double Py, double Pz, double Nx, double Ny, double Nz, double p2x, double p2y, double p2z)
{
- alphaImage = image;
+std::cout << "CFT creaVtkHeartAngles::calculateImages Start"<<std::endl;
+
+ alphaImage = vtkImageData::New();
alphaImage->SetExtent( image->GetExtent() );
alphaImage->SetOrigin( image->GetOrigin() );
alphaImage->SetSpacing( image->GetSpacing() );
alphaImage->SetNumberOfScalarComponents( image->GetNumberOfScalarComponents() );
alphaImage->AllocateScalars();
- betaImage = image;
+ int ext[6];
+ image->GetWholeExtent(ext);
+ int dim[3];
+ dim[0]=ext[1]-ext[0]+1;
+ dim[1]=ext[3]-ext[2]+1;
+ dim[2]=ext[5]-ext[4]+1;
+
+ std::cout<<"dim0 "<<dim[0]<<" dim1 "<<dim[1]<<" dim2 "<<dim[2]<<std::endl;
+
+ betaImage = vtkImageData::New();
betaImage->SetExtent( image->GetExtent() );
+ betaImage->SetWholeExtent( image->GetWholeExtent() );
betaImage->SetOrigin( image->GetOrigin() );
betaImage->SetSpacing( image->GetSpacing() );
betaImage->SetScalarTypeToUnsignedChar();
double a;
double b;
- for( int i = 0 ; i < image->GetExtent()[1] ; i++ )
+ int numberOfPoints = dim[0]*dim[1]*dim[2];
+
+ vtkDoubleArray *array;
+ array = (vtkDoubleArray*) image->GetPointData()->GetVectors();
+ //for(int i = 0 ; i < 8600 ; i++) {
+ // std::cout<<" arreglo #"<< i << " "<< array->GetValue(i)<<" tuple3 "<< array->GetTuple3(i)[0]<<" "<< array->GetTuple3(i)[1]<<" "<<array->GetTuple3(i)[2]<< std::endl;
+ //}
+
+ int i, j, k;
+ int numTuple=0;
+
+
+ for( k = 0 ; k < dim[2] ; k++ )
{
- for( int j = 0 ; j < image->GetExtent()[3] ; j++ )
+ for( j = 0 ; j < dim[1] ; j++ )
{
- for( int k = 0 ; k < image->GetExtent()[5] ; k++ )
+ for( i = 0 ; i < dim[0] ; i++ )
{
- unsigned char* pPix = (unsigned char *)image->GetScalarPointer( i , j , k );
- int p1 = pPix[0];
- int p2 = pPix[1];
- int p3 = pPix[2];
-
+ //unsigned char* pPix = (unsigned char *)image->GetScalarPointer( i , j , k );
+
+ double p1 = array->GetTuple3(numTuple)[0];
+ double p2 = array->GetTuple3(numTuple)[1];
+ double p3 = array->GetTuple3(numTuple)[2];
+ numTuple++;
+
+ //std::cout << " p1: " << p1 << " p2: " << p2 <<" p3: " << p3 << std::endl;
+ //std::cout << " Px: " << Px << " Py: " << Py <<" Pz: " << Pz << std::endl;
+ //std::cout << " p2x: " << p2x << " p2y: " << p2y <<" p2z: " << p2z << std::endl;
+
double* vecProj;
- vecProj = vectorProjection (Px, Py, Pz,
- Nx, Ny, Nz,
+ vecProj = vectorProjection (Nx, Ny, Nz,
p1, p2, p3);
+ //std::cout << " vecProj[0]: " << vecProj[0] << " vecProj[1]: " << vecProj[1] <<" vecProj[2]: " << vecProj[2] << std::endl;
a = alpha (Px, Py, Pz,
p1, p2, p3,
vecProj[0], vecProj[1], vecProj[2]);
p2x, p2y, p2z,
vecProj[0], vecProj[1], vecProj[2]);
- unsigned char *zPtr1 = (unsigned char *) alphaImage->GetScalarPointer( i , j , k );
- *zPtr1 = a;
- unsigned char *zPtr2 = (unsigned char *) betaImage->GetScalarPointer( i , j , k );
- *zPtr2 = b;
+ //if(a>=90&&a<=360) { std::cout << " numTuple: " << numTuple << " a: " << a << " b: " << b << std::endl; }
+ if(p1==0&&p2==0&&p3==0){
+ unsigned char *zPtr1 = (unsigned char *) alphaImage->GetScalarPointer( i , j , k );
+ *zPtr1 = (unsigned char)255;
+ unsigned char *zPtr2 = (unsigned char *) betaImage->GetScalarPointer( i , j , k );
+ *zPtr2 = (unsigned char)255;
+ }else{
+ unsigned char *zPtr1 = (unsigned char *) alphaImage->GetScalarPointer( i , j , k );
+ *zPtr1 = (unsigned char)a;
+ unsigned char *zPtr2 = (unsigned char *) betaImage->GetScalarPointer( i , j , k );
+ *zPtr2 = (unsigned char)b;
}
- }
- }
+ } //for i
+ } // for j
+ } // for k
+
+
+std::cout << "CFT creaVtkHeartAngles::calculateImages End"<<std::endl;
}
git://git.creatis.insa-lyon.fr / creaVtk.git / blobdiff