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)) ) ));
63 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) ) );
65 alpha = (180 * acos(fabs(cosAlpha)))/M_PI;
66 /*if((alpha<=90)&&(alpha>=0)){
67 std::cout << " " << a << " " << b <<" " << c <<" " << x <<" " << y <<" " << z <<" " << xp <<" " << yp <<" " << zp;
68 std::cout << " cosAlpha: " << cosAlpha;
69 std::cout << " acos(fabs(cosalpha)) " << acos(fabs(cosAlpha));
70 std::cout << " Angulo Alpha: " << alpha << std::endl ; }*/
71 double ent = floor(alpha);
72 double al = alpha - ent;
74 //std::cout << "minutos: " << min << std::endl;
75 //std::cout << "Alpha grados: "<< ent <<" minutos: " << floor(min) << " cos: " << cosAlpha << std::endl;
79 double creaVtkHeartAngles::beta (double P0a, double P0b,double P0c, double P3x, double P3y, double P3z, double P4x, double P4y, double P4z)
81 double a=P0a,b=P0b,c=P0c;
82 double x=P3x,y=P3y,z=P3z;
83 double xp=P4x,yp=P4y,zp=P4z;
86 //Calcular el ángulo que forman las rectas, sabiendo sus vectores directores.
87 //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)) ) ));
88 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) ) );
89 beta = (180 * acos(fabs(cosBeta)))/M_PI;
90 //std::cout << "Angulo Beta: " << beta << std::endl;
91 double ent = floor(beta);
92 double al = beta - ent;
94 //std::cout << "minutos: " << min << std::endl;
95 //std::cout << "Beta grados: "<< ent <<" minutos: " << floor(min) << " cos: " << cosBeta << std::endl;
99 double * creaVtkHeartAngles::vectorProjection (double nX, double nY, double nZ, double vX, double vY, double vZ)
103 double nx=nX, ny=nY, nz=nZ;
105 double vx=vX, vy=vY, vz=vZ;
110 proj[0]= vX - (nX*vX);
111 proj[1]= vY - (nY*vY);
112 proj[2]= vZ - (nZ*vZ);
117 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)
119 double puntoCorte[3];
121 double Plx1 = plx1, Ply1 = ply1, Plz1 = plz1, Plx2 = plx2, Ply2 = ply2, Plz2 = plz2, Plx3 = plx3, Ply3 = ply3, Plz3 = plz3;
122 double Px1 = px1, Py1 = py1, Pz1 = pz1, Px2 = px2, Py2 = py2, Pz2 = pz2;
126 double aX = ((Ply2-Ply1)*(Plz3-Plz1)-(Ply3-Ply1)*(Plz2-Plz1));
127 double a = (Plx1)*((Ply2-Ply1)*(Plz3-Plz1)-(Ply3-Ply1)*(Plz2-Plz1));
129 double bX = ((Plx2-Plx1)*(Plz3-Plz1)-(Plx3-Plx1)*(Plz2-Plz1));
130 double b = (Ply1)*((Plx2-Plx1)*(Plz3-Plz1)-(Plx3-Plx1)*(Plz2-Plz1));
132 double cX = ((Plx2-Plx1)*(Ply3-Ply1)-(Plx3-Plx1)*(Ply2-Ply1));
133 double c = (Plz1)*((Plx2-Plx1)*(Ply3-Ply1)-(Plx3-Plx1)*(Ply2-Ply1));
135 double r = a + b + c;
136 std::cout << "ecuacion: " << aX << "," << a << "," << bX << "," << b << "," << cX << "," << c << std::endl;
137 std::cout << "r: " << r << std::endl;
139 //Calculo del vector director
141 double vdirx = Px2 - Px1;
142 double vdiry = Py2 - Py1;
143 double vdirz = Pz2 - Pz1;
145 std::cout << "vdirx: " << vdirx << " vdiry: " << vdiry << " vdirz: " << vdirz << std::endl;
147 //Se igualan las formulas del plano y la recta
148 //x=Px1+vdirx, y=Py1+vdiry, z=Pz1+vdirz
150 double eRes = aX*Px1 + bX*Py1 + cX*Pz1;
151 double eResT = aX*vdirx + bX*vdiry + cX*vdirz;
152 std::cout << "eRes: " << eRes << " eResT: " << eResT << std::endl;
153 //Calculo del punto solucion
155 std::cout << "Soluciones infinitas... " << std::endl;
158 double t = (((-1)*eRes)-r)/eResT;
160 double Px = Px1 + (t * vdirx);
161 double Py = Py1 + (t * vdiry);
162 double Pz = Pz1 + (t * vdirz);
164 std::cout << "px: " << Px << " py: " << Py << " pz: " << Pz << std::endl;
173 void creaVtkHeartAngles::calculateImages (vtkImageData* image, double Px, double Py, double Pz, double Nx, double Ny, double Nz, double p2x, double p2y, double p2z)
175 std::cout << "CFT creaVtkHeartAngles::calculateImages Start"<<std::endl;
177 alphaImage = vtkImageData::New();
178 alphaImage->SetExtent( image->GetExtent() );
179 alphaImage->SetOrigin( image->GetOrigin() );
180 alphaImage->SetSpacing( image->GetSpacing() );
181 alphaImage->SetScalarTypeToUnsignedChar();
182 alphaImage->SetNumberOfScalarComponents( image->GetNumberOfScalarComponents() );
183 alphaImage->AllocateScalars();
186 image->GetWholeExtent(ext);
188 dim[0]=ext[1]-ext[0]+1;
189 dim[1]=ext[3]-ext[2]+1;
190 dim[2]=ext[5]-ext[4]+1;
192 std::cout<<"dim0 "<<dim[0]<<" dim1 "<<dim[1]<<" dim2 "<<dim[2]<<std::endl;
194 betaImage = vtkImageData::New();
195 betaImage->SetExtent( image->GetExtent() );
196 betaImage->SetWholeExtent( image->GetWholeExtent() );
197 betaImage->SetOrigin( image->GetOrigin() );
198 betaImage->SetSpacing( image->GetSpacing() );
199 betaImage->SetScalarTypeToUnsignedChar();
200 betaImage->SetNumberOfScalarComponents( image->GetNumberOfScalarComponents() );
201 betaImage->AllocateScalars();
204 pPix = (unsigned char*)image->GetScalarPointer();
208 int numberOfPoints = dim[0]*dim[1]*dim[2];
210 vtkDoubleArray *array;
211 array = (vtkDoubleArray*) image->GetPointData()->GetVectors();
212 //for(int i = 0 ; i < 8600 ; i++) {
213 // std::cout<<" arreglo #"<< i << " "<< array->GetValue(i)<<" tuple3 "<< array->GetTuple3(i)[0]<<" "<< array->GetTuple3(i)[1]<<" "<<array->GetTuple3(i)[2]<< std::endl;
220 for( k = 0 ; k < dim[2] ; k++ )
222 for( j = 0 ; j < dim[1] ; j++ )
224 for( i = 0 ; i < dim[0] ; i++ )
226 //unsigned char* pPix = (unsigned char *)image->GetScalarPointer( i , j , k );
228 double p1 = array->GetTuple3(numTuple)[0];
229 double p2 = array->GetTuple3(numTuple)[1];
230 double p3 = array->GetTuple3(numTuple)[2];
233 //std::cout << " p1: " << p1 << " p2: " << p2 <<" p3: " << p3 << std::endl;
234 //std::cout << " Px: " << Px << " Py: " << Py <<" Pz: " << Pz << std::endl;
235 //std::cout << " p2x: " << p2x << " p2y: " << p2y <<" p2z: " << p2z << std::endl;
238 vecProj = vectorProjection (Nx, Ny, Nz, p1, p2, p3);
239 //std::cout << " vecProj[0]: " << vecProj[0] << " vecProj[1]: " << vecProj[1] <<" vecProj[2]: " << vecProj[2] << std::endl;
240 a = alpha (Px, Py, Pz,
242 vecProj[0], vecProj[1], vecProj[2]);
244 b = beta (Px, Py, Pz,
246 vecProj[0], vecProj[1], vecProj[2]);
248 //if(a>=90&&a<=360) { std::cout << " numTuple: " << numTuple << " a: " << a << " b: " << b << std::endl; }
249 if(p1==0&&p2==0&&p3==0){
250 unsigned char *zPtr1 = (unsigned char *) alphaImage->GetScalarPointer( i , j , k );
251 *zPtr1 = (unsigned char)255;
252 unsigned char *zPtr2 = (unsigned char *) betaImage->GetScalarPointer( i , j , k );
253 *zPtr2 = (unsigned char)255;
255 unsigned char *zPtr1 = (unsigned char *) alphaImage->GetScalarPointer( i , j , k );
256 *zPtr1 = (unsigned char)a;
257 unsigned char *zPtr2 = (unsigned char *) betaImage->GetScalarPointer( i , j , k );
258 *zPtr2 = (unsigned char)b;
264 std::cout << "CFT creaVtkHeartAngles::calculateImages End"<<std::endl;