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"
31 #include "vtkDoubleArray.h"
32 #include "vtkStructuredPoints.h"
33 #include "vtkPointData.h"
34 #include "vtkDataArray.h"
36 creaVtkHeartAngles::creaVtkHeartAngles()
40 creaVtkHeartAngles::~creaVtkHeartAngles()
44 vtkImageData* creaVtkHeartAngles::getAlphaImage ()
49 vtkImageData* creaVtkHeartAngles::getBetaImage ()
54 double creaVtkHeartAngles::alpha (double P0a, double P0b,double P0c, double vx, double vy, double vz, double vxp, double vyp, double vzp)
56 double a=P0a,b=P0b,c=P0c;
57 double x=vx,y=vy,z=vz;
58 double xp=vxp,yp=vyp,zp=vzp;
61 //Calcular el ángulo que forman las rectas, sabiendo sus vectores directores.
62 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)) ) ));
64 alpha = (180 * acos(fabs(cosAlpha)))/M_PI;
65 //std::cout << "Angulo Alpha: " << alpha << std::endl;
66 double ent = floor(alpha);
67 double al = alpha - ent;
69 //std::cout << "minutos: " << min << std::endl;
70 //std::cout << "Alpha grados: "<< ent <<" minutos: " << floor(min) << " cos: " << cosAlpha << std::endl;
74 double creaVtkHeartAngles::beta (double P0a, double P0b,double P0c, double P3x, double P3y, double P3z, double P4x, double P4y, double P4z)
76 double a=P0a,b=P0b,c=P0c;
77 double x=P3x,y=P3y,z=P3z;
78 double xp=P4x,yp=P4y,zp=P4z;
81 //Calcular el ángulo que forman las rectas, sabiendo sus vectores directores.
82 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)) ) ));
84 beta = (180 * acos(fabs(cosBeta)))/M_PI;
85 //std::cout << "Angulo Beta: " << beta << std::endl;
86 double ent = floor(beta);
87 double al = beta - ent;
89 //std::cout << "minutos: " << min << std::endl;
90 //std::cout << "Beta grados: "<< ent <<" minutos: " << floor(min) << " cos: " << cosBeta << std::endl;
94 double * creaVtkHeartAngles::vectorProjection (double nX, double nY, double nZ, double vX, double vY, double vZ)
98 double nx=nX, ny=nY, nz=nZ;
100 double vx=vX, vy=vY, vz=vZ;
105 proj[0]= vX - (nX*vX);
106 proj[1]= vY - (nY*vY);
107 proj[2]= vZ - (nZ*vZ);
112 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)
114 double puntoCorte[3];
116 double Plx1 = plx1, Ply1 = ply1, Plz1 = plz1, Plx2 = plx2, Ply2 = ply2, Plz2 = plz2, Plx3 = plx3, Ply3 = ply3, Plz3 = plz3;
117 double Px1 = px1, Py1 = py1, Pz1 = pz1, Px2 = px2, Py2 = py2, Pz2 = pz2;
121 double aX = ((Ply2-Ply1)*(Plz3-Plz1)-(Ply3-Ply1)*(Plz2-Plz1));
122 double a = (Plx1)*((Ply2-Ply1)*(Plz3-Plz1)-(Ply3-Ply1)*(Plz2-Plz1));
124 double bX = ((Plx2-Plx1)*(Plz3-Plz1)-(Plx3-Plx1)*(Plz2-Plz1));
125 double b = (Ply1)*((Plx2-Plx1)*(Plz3-Plz1)-(Plx3-Plx1)*(Plz2-Plz1));
127 double cX = ((Plx2-Plx1)*(Ply3-Ply1)-(Plx3-Plx1)*(Ply2-Ply1));
128 double c = (Plz1)*((Plx2-Plx1)*(Ply3-Ply1)-(Plx3-Plx1)*(Ply2-Ply1));
130 double r = a + b + c;
131 std::cout << "ecuacion: " << aX << "," << a << "," << bX << "," << b << "," << cX << "," << c << std::endl;
132 std::cout << "r: " << r << std::endl;
134 //Calculo del vector director
136 double vdirx = Px2 - Px1;
137 double vdiry = Py2 - Py1;
138 double vdirz = Pz2 - Pz1;
140 std::cout << "vdirx: " << vdirx << " vdiry: " << vdiry << " vdirz: " << vdirz << std::endl;
142 //Se igualan las formulas del plano y la recta
143 //x=Px1+vdirx, y=Py1+vdiry, z=Pz1+vdirz
145 double eRes = aX*Px1 + bX*Py1 + cX*Pz1;
146 double eResT = aX*vdirx + bX*vdiry + cX*vdirz;
147 std::cout << "eRes: " << eRes << " eResT: " << eResT << std::endl;
148 //Calculo del punto solucion
150 std::cout << "Soluciones infinitas... " << std::endl;
153 double t = (((-1)*eRes)-r)/eResT;
155 double Px = Px1 + (t * vdirx);
156 double Py = Py1 + (t * vdiry);
157 double Pz = Pz1 + (t * vdirz);
159 std::cout << "px: " << Px << " py: " << Py << " pz: " << Pz << std::endl;
168 void creaVtkHeartAngles::calculateImages (vtkImageData* image, double Px, double Py, double Pz, double Nx, double Ny, double Nz, double p2x, double p2y, double p2z)
170 std::cout << "CFT creaVtkHeartAngles::calculateImages Start"<<std::endl;
172 alphaImage = vtkImageData::New();
173 alphaImage->SetExtent( image->GetExtent() );
174 alphaImage->SetOrigin( image->GetOrigin() );
175 alphaImage->SetSpacing( image->GetSpacing() );
176 alphaImage->SetScalarTypeToUnsignedChar();
177 alphaImage->SetNumberOfScalarComponents( image->GetNumberOfScalarComponents() );
178 alphaImage->AllocateScalars();
181 image->GetWholeExtent(ext);
183 dim[0]=ext[1]-ext[0]+1;
184 dim[1]=ext[3]-ext[2]+1;
185 dim[2]=ext[5]-ext[4]+1;
187 std::cout<<"dim0 "<<dim[0]<<" dim1 "<<dim[1]<<" dim2 "<<dim[2]<<std::endl;
189 betaImage = vtkImageData::New();
190 betaImage->SetExtent( image->GetExtent() );
191 betaImage->SetWholeExtent( image->GetWholeExtent() );
192 betaImage->SetOrigin( image->GetOrigin() );
193 betaImage->SetSpacing( image->GetSpacing() );
194 betaImage->SetScalarTypeToUnsignedChar();
195 betaImage->SetNumberOfScalarComponents( image->GetNumberOfScalarComponents() );
196 betaImage->AllocateScalars();
199 pPix = (unsigned char*)image->GetScalarPointer();
203 int numberOfPoints = dim[0]*dim[1]*dim[2];
205 vtkDoubleArray *array;
206 array = (vtkDoubleArray*) image->GetPointData()->GetVectors();
207 //for(int i = 0 ; i < 8600 ; i++) {
208 // std::cout<<" arreglo #"<< i << " "<< array->GetValue(i)<<" tuple3 "<< array->GetTuple3(i)[0]<<" "<< array->GetTuple3(i)[1]<<" "<<array->GetTuple3(i)[2]<< std::endl;
215 for( i = 0 ; i < dim[0] ; i++ )
217 for( j = 0 ; j < dim[1] ; j++ )
219 for( k = 0 ; k < dim[2] ; k++ )
221 //unsigned char* pPix = (unsigned char *)image->GetScalarPointer( i , j , k );
223 double p1 = array->GetTuple3(numTuple)[0];
224 double p2 = array->GetTuple3(numTuple)[1];
225 double p3 = array->GetTuple3(numTuple)[2];
228 //std::cout << " p1: " << p1 << " p2: " << p2 <<" p3: " << p3 << std::endl;
229 //std::cout << " Px: " << Px << " Py: " << Py <<" Pz: " << Pz << std::endl;
230 //std::cout << " p2x: " << p2x << " p2y: " << p2y <<" p2z: " << p2z << std::endl;
233 vecProj = vectorProjection (Nx, Ny, Nz,
235 //std::cout << " vecProj[0]: " << vecProj[0] << " vecProj[1]: " << vecProj[1] <<" vecProj[2]: " << vecProj[2] << std::endl;
236 a = alpha (Px, Py, Pz,
238 vecProj[0], vecProj[1], vecProj[2]);
240 b = beta (Px, Py, Pz,
242 vecProj[0], vecProj[1], vecProj[2]);
243 //std::cout << " a: " << a << " b: " << b << std::endl;
244 unsigned char *zPtr1 = (unsigned char *) alphaImage->GetScalarPointer( i , j , k );
245 *zPtr1 = (unsigned char)a;
246 unsigned char *zPtr2 = (unsigned char *) betaImage->GetScalarPointer( i , j , k );
247 *zPtr2 = (unsigned char)b;
253 std::cout << "CFT creaVtkHeartAngles::calculateImages End"<<std::endl;