]> Creatis software - creaVtk.git/blobdiff - lib/creaVtk/creaVtkHeartAngles.cpp
2209 creaVtk Feature New Normal Correction of Heart Angles Box 2013-12-13 16:50
[creaVtk.git] / lib / creaVtk / creaVtkHeartAngles.cpp
index 916f0b7190426d90c60115b3165941541a03a609..c70472e7df1760e919205b5c30e117d424fd4d1d 100644 (file)
 #include "creaVtkHeartAngles.h"
 #include <math.h>
 #include <iostream>
+#include "vtkDoubleArray.h"
+#include "vtkStructuredPoints.h"
+#include "vtkPointData.h"
+#include "vtkDataArray.h"
 
 creaVtkHeartAngles::creaVtkHeartAngles()
 {
@@ -37,6 +41,16 @@ 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;
@@ -53,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;
 }
 
@@ -73,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;
 }
@@ -161,55 +165,93 @@ double *creaVtkHeartAngles::intersectionPlaneLine(double plx1, double ply1, doub
        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)
 {
+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;
 
-       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]);
-
-       }
-       
-       
-}
+       int numberOfPoints = dim[0]*dim[1]*dim[2];
 
-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();
+       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;
+       //}
 
-       for(int i=0;i<1041;i=i+3){
-               int p1 = pPix[i];
-               int p2 = pPix[i+1];
-               int p3 = pPix[i+2];
+       int i, j, k;
+       int numTuple=0;
 
-       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);
+       for( i = 0 ; i < dim[0] ; i++ )
+  {
+    for( j = 0 ; j < dim[1] ; j++ )
+       {
+                       for( k = 0 ; k < dim[2] ; k++ )
+               {
+                       //unsigned char* pPix = (unsigned char *)image->GetScalarPointer( i , j , k );
 
-       b = beta (0, 0, 0, 
-                                               p1, p2, p3, //?? es el vector?
-                                               vecProj[0], vecProj[1], vecProj[2]);
+                       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]);
+                       //std::cout << " a: " << a << " b: " << b << std::endl;
+                       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 k
+       } // for j
+  } // for i
+
+
+std::cout << "CFT creaVtkHeartAngles::calculateImages End"<<std::endl;
 }
 
+
+