#include "creaVtkHeartAngles.h"
#include <math.h>
#include <iostream>
+#include "vtkDoubleArray.h"
+#include "vtkStructuredPoints.h"
+#include "vtkPointData.h"
+#include "vtkDataArray.h"
creaVtkHeartAngles::creaVtkHeartAngles()
{
{
}
+vtkImageData* creaVtkHeartAngles::getAlphaImage ()
+{
+ return alphaImage;
+}
+
+vtkImageData* creaVtkHeartAngles::getBetaImage ()
+{
+ return betaImage;
+};
+
double creaVtkHeartAngles::alpha (double P0a, double P0b,double P0c, double vx, double vy, double vz, double vxp, double vyp, double vzp)
{
double a=P0a,b=P0b,c=P0c;
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;
}
-vtkImageData* creaVtkHeartAngles::calculateAngleAlpha (vtkImageData* image, double Px, double Py, double Pz, double Vx, double Vy, double Vz, double Nx, double Ny, double Nz /*, double* pPlane*/)
+void creaVtkHeartAngles::calculateImages (vtkImageData* image, double Px, double Py, double Pz, double Nx, double Ny, double Nz, double p2x, double p2y, double p2z)
{
- unsigned char* pPix;
- pPix = (unsigned char*)image->GetScalarPointer();
-
- for(int i=0;i<1041;i=i+3){
- int p1 = pPix[i];
- int p2 = pPix[i+1];
- int p3 = pPix[i+2];
-
- double* vecProj;
- vecProj = vectorProjection (Px, Py, Pz, Nx, Ny, Nz, p1, p2, p3);
- double a, b;
-
- //double* cutPoint;
- //cutPoint = intersectionPlaneLine(pPlane[0], pPlane[1], pPlane[2], pPlane[3], pPlane[4], pPlane[5], pPlane[6], pPlane[7], pPlane[8], px1, py1, pz1, px2, py2, pz2);
-
- a = alpha (p1, p2, p3,
- p1, p2, p3, //?? es el vector?
- vecProj[0], vecProj[1], vecProj[2]);
-
- }
-
-
+std::cout << "CFT creaVtkHeartAngles::calculateImages Start"<<std::endl;
+
+ alphaImage = vtkImageData::New();
+ alphaImage->SetExtent( image->GetExtent() );
+ alphaImage->SetOrigin( image->GetOrigin() );
+ alphaImage->SetSpacing( image->GetSpacing() );
+ alphaImage->SetScalarTypeToUnsignedChar();
+ alphaImage->SetNumberOfScalarComponents( image->GetNumberOfScalarComponents() );
+ alphaImage->AllocateScalars();
+
+ 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();
+ betaImage->SetNumberOfScalarComponents( image->GetNumberOfScalarComponents() );
+ betaImage->AllocateScalars();
+
+ unsigned char* pPix;
+ pPix = (unsigned char*)image->GetScalarPointer();
+ double a;
+ double b;
+
+ 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( j = 0 ; j < dim[1] ; j++ )
+ {
+ for( i = 0 ; i < dim[0] ; i++ )
+ {
+ //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 (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]);
+
+ b = beta (Px, Py, Pz,
+ p2x, p2y, p2z,
+ vecProj[0], vecProj[1], vecProj[2]);
+
+ //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;
}
-vtkImageData* creaVtkHeartAngles::calculateAngleBeta (vtkImageData* image, double Px, double Py, double Pz, double Vx, double Vy, double Vz, double Nx, double Ny, double Nz /*, double* pPlane*/)
-{
- unsigned char* pPix;
- pPix = (unsigned char*)image->GetScalarPointer();
-
- for(int i=0;i<1041;i=i+3){
- int p1 = pPix[i];
- int p2 = pPix[i+1];
- int p3 = pPix[i+2];
-
- double* vecProj;
- vecProj = vectorProjection (Px, Py, Pz, Nx, Ny, Nz, p1, p2, p3);
- double a, b;
- //double* cutPoint;
- //cutPoint = intersectionPlaneLine(pPlane[0], pPlane[1], pPlane[2], pPlane[3], pPlane[4], pPlane[5], pPlane[6], pPlane[7], pPlane[8], px1, py1, pz1, px2, py2, pz2);
-
- b = beta (0, 0, 0,
- p1, p2, p3, //?? es el vector?
- vecProj[0], vecProj[1], vecProj[2]);
-
- }
-
-
-}
git://git.creatis.insa-lyon.fr / creaVtk.git / blobdiff