1 /*# ---------------------------------------------------------------------
3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
5 # Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
6 # Previous Authors : Laurent Guigues, Jean-Pierre Roux
7 # CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
9 # This software is governed by the CeCILL-B license under French law and
10 # abiding by the rules of distribution of free software. You can use,
11 # modify and/ or redistribute the software under the terms of the CeCILL-B
12 # license as circulated by CEA, CNRS and INRIA at the following URL
13 # http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
14 # or in the file LICENSE.txt.
16 # As a counterpart to the access to the source code and rights to copy,
17 # modify and redistribute granted by the license, users are provided only
18 # with a limited warranty and the software's author, the holder of the
19 # economic rights, and the successive licensors have only limited
22 # The fact that you are presently reading this means that you have had
23 # knowledge of the CeCILL-B license and that you accept its terms.
24 # ------------------------------------------------------------------------ */
26 #include "CutModelPolygon.h"
29 CutModelPolygon::CutModelPolygon()
33 _transform = vtkTransform::New();
36 CutModelPolygon::~CutModelPolygon()
48 void CutModelPolygon::processOutImage(int cutInsideOutside)
50 _cutInsideOutside=cutInsideOutside;
52 int numPoints = _points->GetNumberOfPoints();
56 // Calculate Orthogonal Vectors using two vectors in the plane and its cross product
57 calculateOrthogonalVectors(v1,v2);
59 cout<<"PolyCutter::processOutImage"<<endl;
60 cout<<"V1:"<<" X:"<<v1[0]<<" Y:"<<v1[1]<<" Z:"<<v1[2]<<endl;
61 cout<<"V2:"<<" X:"<<v2[0]<<" Y:"<<v2[1]<<" Z:"<<v2[2]<<endl;
62 cout<<"Direction:"<<" X:"<<_direction[0]<<" Y:"<<_direction[1]<<" Z:"<<_direction[2]<<endl;
65 // Update the inverse transform which it's applied to all the points of the contour
66 // and the points of the input image.
67 updateTransform(v1,v2,_direction);
69 //Transform Points coordinates
70 std::vector<double> vectorOutX;
71 std::vector<double> vectorOutY;
72 std::vector<double> vectorOutZ;
74 transformContourPoints(&vectorOutX,&vectorOutY,&vectorOutZ);
78 int num = vectorOutX.size();
79 for(int t=0;t<num;t++)
81 cout<<"Final Point"<<t<<"-X:"<<vectorOutX[t]<<" Y:"<<vectorOutY[t]<<" Z:"<<vectorOutZ[t]<<endl;
84 cutInputImage(vectorOutX,vectorOutY,vectorOutZ);
88 //-------------------------------------------------------------------------
89 void CutModelPolygon::cutInputImage(std::vector<double> vectorOutX,std::vector<double> vectorOutY,std::vector<double> vectorOutZ)
92 // Find the minimum value in Y axis
94 double minY=vectorOutY[0];
95 for(i=1;i<vectorOutY.size();i++)
97 if(vectorOutY[i]<minY)
103 // All the contour points minus the minimum to translate to the positive octant
105 for(i=0;i<vectorOutY.size();i++)
108 vectorOutY[i]=vectorOutY[i]-minY;
112 // Creates a poligonal contour model with the points
113 manualBaseModel *model = InitializeContourModel(vectorOutX,vectorOutY,vectorOutZ);
115 initializeOutputImage();
118 //EED 2017-01-01 Migration VTK7
119 #if VTK_MAJOR_VERSION <= 5
120 _inImage->GetWholeExtent(ext);
122 _inImage->GetExtent(ext);
124 int dimX=ext[1]-ext[0]+1;
125 int dimY=ext[3]-ext[2]+1;
126 int dimZ=ext[5]-ext[4]+1;
129 //CreaMaracasVisu Class which evaluates whether a point is inside the contour or not
130 ContourExtractData *extract = new ContourExtractData(false);
132 //Contour list just with the contour created with the points
133 std::vector<manualBaseModel*> list;
134 list.push_back(model);
136 //Y-value asigned arbitrary
137 extract->SetSizeImageY(500);
140 extract->SetLstManualContourModel(list);
141 extract->InitLstContoursLinesYPoints();
143 cout<<"Processing..."<<endl;
145 for (int i=0; i<dimX; i++)
148 cout<<"Print "<<i<<" of "<<dimX<<endl;
150 for (int j=0; j<dimY; j++)
152 for (int k=0; k<dimZ; k++)
154 unsigned short *pImage=(unsigned short*)_inImage->GetScalarPointer(i,j,k);
155 //unsigned short *pResult=(unsigned short*)_outImage->GetScalarPointer(i,j,k);
162 //Transform every point in the imageData to the space resultant of the mathematical transform
163 _transform->MultiplyPoint(in,out);
165 //All the imageData minus the minimum to translate to the positive octant
168 //Verify if the point is inside the contour or not
169 if ((out[1]>=0 && out[1]<500) && (extract->isInside(out[0],out[1],0)==true))
172 if(_cutInsideOutside==0)//CutInside
180 if(_cutInsideOutside==1)//CutOutside
190 _inImage->Modified();
191 //EED 2017-01-01 Migration VTK7
192 #if VTK_MAJOR_VERSION <= 5
200 //-------------------------------------------------------------------------
201 void CutModelPolygon::transformContourPoints(std::vector<double> *vectorOutX,std::vector<double> *vectorOutY,std::vector<double> *vectorOutZ)
203 int numPoints = _points->GetNumberOfPoints();
205 for(int t=0;t<numPoints;t++)
207 double point[3],in[4],out[4];
208 _points->GetPoint(t,point);
214 //CHG_LYON_JAN29 Multiply Contour Points by the inverse of the spacing. It's necessary to faire
215 //the hole in the correct position and not in other with a displacement
218 _inImage->GetSpacing(_spc);
219 in[0]=in[0]*(1/_spc[0]);
220 in[1]=in[1]*(1/_spc[1]);
221 in[2]=in[2]*(1/_spc[2]);
223 _transform->MultiplyPoint(in,out);
224 vectorOutX->push_back(out[0]);
225 vectorOutY->push_back(out[1]);
226 vectorOutZ->push_back(out[2]);
230 //-------------------------------------------------------------------------
231 void CutModelPolygon::updateTransform(double* v1, double* v2, double* v3)
234 // ImageData Centroid
235 int i, numPoints=_points->GetNumberOfPoints();
236 double centerX=0,centerY=0,centerZ=0;
237 for(i =0; i<numPoints;i++)
240 _points->GetPoint(i,point);
246 centerX=centerX/numPoints;
247 centerY=centerY/numPoints;
248 centerZ=centerZ/numPoints;
251 vtkMatrix4x4 *matrix = vtkMatrix4x4::New();
254 matrix->SetElement(0,0,v1[0]);
255 matrix->SetElement(1,0,v1[1]);
256 matrix->SetElement(2,0,v1[2]);
257 matrix->SetElement(3,0,0);
261 matrix->SetElement(0,1,v2[0]);
262 matrix->SetElement(1,1,v2[1]);
263 matrix->SetElement(2,1,v2[2]);
264 matrix->SetElement(3,1,0);
267 matrix->SetElement(0,2,v3[0]);
268 matrix->SetElement(1,2,v3[1]);
269 matrix->SetElement(2,2,v3[2]);
270 matrix->SetElement(3,2,0);
273 matrix->SetElement(0,3,centerX);
274 matrix->SetElement(1,3,centerY);
275 matrix->SetElement(2,3,centerZ);
276 matrix->SetElement(3,3,1);
278 _transform->SetMatrix(matrix);
279 _transform->Update();
281 _transform->Inverse();
282 _transform->Update();
285 //-------------------------------------------------------------------------
286 void CutModelPolygon::calculateOrthogonalVectors(double* v1, double* v2)
289 _points->GetPoint(0,p1);
290 _points->GetPoint(1,p2);
298 double magV1 = sqrt( v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2]);
303 v2[0] = v1[1]*_direction[2] - v1[2]*_direction[1];
304 v2[1] = -(v1[0]*_direction[2] - v1[2]*_direction[0]);
305 v2[2] = v1[0]*_direction[1] - v1[1]*_direction[0];
307 double magV2 = sqrt( v2[0]*v2[0] + v2[1]*v2[1] + v2[2]*v2[2]);
312 double magDir = sqrt( _direction[0]*_direction[0] + _direction[1]*_direction[1] + _direction[2]*_direction[2]);
313 _direction[0]=_direction[0]/magDir;
314 _direction[1]=_direction[1]/magDir;
315 _direction[2]=_direction[2]/magDir;
318 //-------------------------------------------------------------------------
319 void CutModelPolygon::initializeOutputImage()
322 //_inImage->GetWholeExtent(ext);
323 //int dimX=ext[1]-ext[0]+1;
324 //int dimY=ext[3]-ext[2]+1;
325 //int dimZ=ext[5]-ext[4]+1;
327 //if (_outImage != NULL ) { _outImage->Delete(); }
328 //_outImage=vtkImageData::New();
329 //_outImage->SetDimensions(dimX,dimY,dimZ);
330 //_outImage->SetScalarTypeToUnsignedShort();
331 //_outImage->AllocateScalars();
334 //-------------------------------------------------------------------------
336 manualBaseModel* CutModelPolygon::InitializeContourModel(std::vector<double> pointsX, std::vector<double> pointsY, std::vector<double> pointsZ)
338 manualBaseModel *model = new manualContourModelPolygon();
340 int numPoints = pointsX.size();
341 for(int i = 0; i<numPoints;i++)
343 model->AddPoint(pointsX[i],pointsY[i],pointsZ[i]);
351 ////////////////////////////////////////////
352 // Getters and setters
353 ////////////////////////////////////////////
355 vtkImageData* CutModelPolygon::getInImage()
360 vtkImageData* CutModelPolygon::getOutImage()
365 vtkPoints* CutModelPolygon::getPoints()
370 double* CutModelPolygon::getDirection()
375 void CutModelPolygon::setInImage(vtkImageData* pImage)
380 void CutModelPolygon::setOutImage(vtkImageData* pImage)
385 void CutModelPolygon::setPoints(vtkPoints *pPoints)
390 void CutModelPolygon::setDirection(double *pDirection)
392 _direction=pDirection;