2217 creaVtk Support New Normal Correction of Heart Angles Box

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

diff --git a/lib/creaVtk/creaVtkHeartAngles.cpp b/lib/creaVtk/creaVtkHeartAngles.cpp
index c70472e7df1760e919205b5c30e117d424fd4d1d..26de410e8144af06499eac2b2591b01fbaba2e21 100644 (file)
--- a/lib/creaVtk/creaVtkHeartAngles.cpp
+++ b/lib/creaVtk/creaVtkHeartAngles.cpp
@@ -59,10 +59,15 @@ double creaVtkHeartAngles::alpha (double P0a, double P0b,double P0c, double vx,
        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;
@@ -79,8 +84,8 @@ double creaVtkHeartAngles::beta (double P0a, double P0b,double P0c, double P3x,
        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);
@@ -195,44 +200,43 @@ std::cout << "CFT creaVtkHeartAngles::calculateImages Start"<<std::endl;
   betaImage->SetNumberOfScalarComponents( image->GetNumberOfScalarComponents() );
   betaImage->AllocateScalars();
 
-       unsigned char* pPix;
-       pPix = (unsigned char*)image->GetScalarPointer();
-       double a;
-       double b;
+  unsigned char* pPix;
+  pPix = (unsigned char*)image->GetScalarPointer();
+  double a;
+  double b;
 
-       int numberOfPoints = dim[0]*dim[1]*dim[2];
+  int numberOfPoints = dim[0]*dim[1]*dim[2];
 
-       vtkDoubleArray *array;
+  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;
+  int i, j, k;
+  int numTuple=0;
 
 
-       for( i = 0 ; i < dim[0] ; i++ )
+  for( k = 0 ; k < dim[2] ; k++ )
   {
     for( j = 0 ; j < dim[1] ; j++ )
+    {
+       for( i = 0 ; i < dim[0] ; i++ )
        {
-                       for( k = 0 ; k < dim[2] ; k++ )
-               {
-                       //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++;
+               //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;                     
+               //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]);
@@ -240,15 +244,22 @@ std::cout << "CFT creaVtkHeartAngles::calculateImages Start"<<std::endl;
                        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
+                       
+                       //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 i
-
+  } // for k
 
 std::cout << "CFT creaVtkHeartAngles::calculateImages End"<<std::endl;
 }