2209 creaVtk Feature New Normal Correction of Heart Angles Box 2013-12-13 16:50

[creaVtk.git] / lib / creaVtk / creaVtkHeartAngles.cpp

diff --git a/lib/creaVtk/creaVtkHeartAngles.cpp b/lib/creaVtk/creaVtkHeartAngles.cpp
index cceaa643963d809df8312c2a1a8e1cc48986a0de..c70472e7df1760e919205b5c30e117d424fd4d1d 100644 (file)
--- a/lib/creaVtk/creaVtkHeartAngles.cpp
+++ b/lib/creaVtk/creaVtkHeartAngles.cpp
@@ -28,6 +28,10 @@
 #include "creaVtkHeartAngles.h"
 #include <math.h>
 #include <iostream>
+#include "vtkDoubleArray.h"
+#include "vtkStructuredPoints.h"
+#include "vtkPointData.h"
+#include "vtkDataArray.h"
 
 creaVtkHeartAngles::creaVtkHeartAngles()
 {
@@ -63,7 +67,7 @@ double creaVtkHeartAngles::alpha (double P0a, double P0b,double P0c, double vx,
        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;
 }
 
@@ -83,34 +87,24 @@ double creaVtkHeartAngles::beta (double P0a, double P0b,double P0c, double P3x,
        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;
 }
@@ -171,9 +165,11 @@ double *creaVtkHeartAngles::intersectionPlaneLine(double plx1, double ply1, doub
        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() );
@@ -181,8 +177,18 @@ void creaVtkHeartAngles::calculateImages (vtkImageData* image, double Px, double
   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();
@@ -194,21 +200,39 @@ void creaVtkHeartAngles::calculateImages (vtkImageData* image, double Px, double
        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( i = 0 ; i < dim[0] ; i++ )
   {
-    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( k = 0 ; k < dim[2] ; k++ )
                {
-                       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]);
@@ -216,14 +240,17 @@ void creaVtkHeartAngles::calculateImages (vtkImageData* image, double Px, double
                        b = beta  (Px, Py, Pz, 
                                                           p2x, p2y, p2z, 
                                                           vecProj[0], vecProj[1], vecProj[2]);
-                       
+                       //std::cout << " a: " << a << " b: " << b << std::endl;
                        unsigned char *zPtr1 = (unsigned char *) alphaImage->GetScalarPointer( i , j , k );
-                       *zPtr1 = a;
+                       *zPtr1 = (unsigned char)a;
                        unsigned char *zPtr2 = (unsigned char *) betaImage->GetScalarPointer( i , j , k );
-                       *zPtr2 = b; 
-                       }
-       } 
-  }
+                       *zPtr2 = (unsigned char)b; 
+                       } //for k
+       } // for j
+  } // for i
+
+
+std::cout << "CFT creaVtkHeartAngles::calculateImages End"<<std::endl;
 }