2 # ---------------------------------------------------------------------
4 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
6 # Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
7 # Previous Authors : Laurent Guigues, Jean-Pierre Roux
8 # CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
10 # This software is governed by the CeCILL-B license under French law and
11 # abiding by the rules of distribution of free software. You can use,
12 # modify and/ or redistribute the software under the terms of the CeCILL-B
13 # license as circulated by CEA, CNRS and INRIA at the following URL
14 # http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
15 # or in the file LICENSE.txt.
17 # As a counterpart to the access to the source code and rights to copy,
18 # modify and redistribute granted by the license, users are provided only
19 # with a limited warranty and the software's author, the holder of the
20 # economic rights, and the successive licensors have only limited
23 # The fact that you are presently reading this means that you have had
24 # knowledge of the CeCILL-B license and that you accept its terms.
25 # ------------------------------------------------------------------------
28 #include "creaVtkHeartAngles.h"
32 creaVtkHeartAngles::creaVtkHeartAngles()
36 creaVtkHeartAngles::~creaVtkHeartAngles()
40 vtkImageData* creaVtkHeartAngles::getAlphaImage ()
45 vtkImageData* creaVtkHeartAngles::getBetaImage ()
50 double creaVtkHeartAngles::alpha (double P0a, double P0b,double P0c, double vx, double vy, double vz, double vxp, double vyp, double vzp)
52 double a=P0a,b=P0b,c=P0c;
53 double x=vx,y=vy,z=vz;
54 double xp=vxp,yp=vyp,zp=vzp;
57 //Calcular el ángulo que forman las rectas, sabiendo sus vectores directores.
58 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)) ) ));
60 alpha = (180 * acos(fabs(cosAlpha)))/M_PI;
61 //std::cout << "Angulo Alpha: " << alpha << std::endl;
62 double ent = floor(alpha);
63 double al = alpha - ent;
65 //std::cout << "minutos: " << min << std::endl;
66 std::cout << "Grados: "<< ent <<" minutos: " << floor(min) << std::endl;
70 double creaVtkHeartAngles::beta (double P0a, double P0b,double P0c, double P3x, double P3y, double P3z, double P4x, double P4y, double P4z)
72 double a=P0a,b=P0b,c=P0c;
73 double x=P3x,y=P3y,z=P3z;
74 double xp=P4x,yp=P4y,zp=P4z;
77 //Calcular el ángulo que forman las rectas, sabiendo sus vectores directores.
78 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)) ) ));
80 beta = (180 * acos(fabs(cosBeta)))/M_PI;
81 //std::cout << "Angulo Beta: " << beta << std::endl;
82 double ent = floor(beta);
83 double al = beta - ent;
85 //std::cout << "minutos: " << min << std::endl;
86 std::cout << "Grados: "<< ent <<" minutos: " << floor(min) << std::endl;
90 double * creaVtkHeartAngles::vectorProjection (double plX, double plY, double plZ, double nX, double nY, double nZ, double vX, double vY, double vZ)
93 double plx=plX, ply=plY, plz=plZ;
95 double nx=nX, ny=nY, nz=nZ;
97 double vx=vX, vy=vY, vz=vZ;
99 //distancia del punto al plano
102 double vpx, vpy, vpz;
106 //std::cout << dist << std
107 vpx = vx - (dist * nx);
108 vpy = vy - (dist * ny);
109 vpz = vz - (dist * nz);
118 double *creaVtkHeartAngles::intersectionPlaneLine(double plx1, double ply1, double plz1, double plx2, double ply2, double plz2, double plx3, double ply3, double plz3, double px1, double py1, double pz1,double px2, double py2, double pz2)
120 double puntoCorte[3];
122 double Plx1 = plx1, Ply1 = ply1, Plz1 = plz1, Plx2 = plx2, Ply2 = ply2, Plz2 = plz2, Plx3 = plx3, Ply3 = ply3, Plz3 = plz3;
123 double Px1 = px1, Py1 = py1, Pz1 = pz1, Px2 = px2, Py2 = py2, Pz2 = pz2;
127 double aX = ((Ply2-Ply1)*(Plz3-Plz1)-(Ply3-Ply1)*(Plz2-Plz1));
128 double a = (Plx1)*((Ply2-Ply1)*(Plz3-Plz1)-(Ply3-Ply1)*(Plz2-Plz1));
130 double bX = ((Plx2-Plx1)*(Plz3-Plz1)-(Plx3-Plx1)*(Plz2-Plz1));
131 double b = (Ply1)*((Plx2-Plx1)*(Plz3-Plz1)-(Plx3-Plx1)*(Plz2-Plz1));
133 double cX = ((Plx2-Plx1)*(Ply3-Ply1)-(Plx3-Plx1)*(Ply2-Ply1));
134 double c = (Plz1)*((Plx2-Plx1)*(Ply3-Ply1)-(Plx3-Plx1)*(Ply2-Ply1));
136 double r = a + b + c;
137 std::cout << "ecuacion: " << aX << "," << a << "," << bX << "," << b << "," << cX << "," << c << std::endl;
138 std::cout << "r: " << r << std::endl;
140 //Calculo del vector director
142 double vdirx = Px2 - Px1;
143 double vdiry = Py2 - Py1;
144 double vdirz = Pz2 - Pz1;
146 std::cout << "vdirx: " << vdirx << " vdiry: " << vdiry << " vdirz: " << vdirz << std::endl;
148 //Se igualan las formulas del plano y la recta
149 //x=Px1+vdirx, y=Py1+vdiry, z=Pz1+vdirz
151 double eRes = aX*Px1 + bX*Py1 + cX*Pz1;
152 double eResT = aX*vdirx + bX*vdiry + cX*vdirz;
153 std::cout << "eRes: " << eRes << " eResT: " << eResT << std::endl;
154 //Calculo del punto solucion
156 std::cout << "Soluciones infinitas... " << std::endl;
159 double t = (((-1)*eRes)-r)/eResT;
161 double Px = Px1 + (t * vdirx);
162 double Py = Py1 + (t * vdiry);
163 double Pz = Pz1 + (t * vdirz);
165 std::cout << "px: " << Px << " py: " << Py << " pz: " << Pz << std::endl;
174 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)
177 alphaImage->SetExtent( image->GetExtent() );
178 alphaImage->SetOrigin( image->GetOrigin() );
179 alphaImage->SetSpacing( image->GetSpacing() );
180 alphaImage->SetScalarTypeToUnsignedChar();
181 alphaImage->SetNumberOfScalarComponents( image->GetNumberOfScalarComponents() );
182 alphaImage->AllocateScalars();
185 betaImage->SetExtent( image->GetExtent() );
186 betaImage->SetOrigin( image->GetOrigin() );
187 betaImage->SetSpacing( image->GetSpacing() );
188 betaImage->SetScalarTypeToUnsignedChar();
189 betaImage->SetNumberOfScalarComponents( image->GetNumberOfScalarComponents() );
190 betaImage->AllocateScalars();
193 pPix = (unsigned char*)image->GetScalarPointer();
197 for( int i = 0 ; i < image->GetExtent()[1] ; i++ )
199 for( int j = 0 ; j < image->GetExtent()[3] ; j++ )
201 for( int k = 0 ; k < image->GetExtent()[5] ; k++ )
203 unsigned char* pPix = (unsigned char *)image->GetScalarPointer( i , j , k );
209 vecProj = vectorProjection (Px, Py, Pz,
212 a = alpha (Px, Py, Pz,
214 vecProj[0], vecProj[1], vecProj[2]);
216 b = beta (Px, Py, Pz,
218 vecProj[0], vecProj[1], vecProj[2]);
220 unsigned char *zPtr1 = (unsigned char *) alphaImage->GetScalarPointer( i , j , k );
222 unsigned char *zPtr2 = (unsigned char *) betaImage->GetScalarPointer( i , j , k );