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 _inImage->GetWholeExtent(ext);
119 int dimX=ext[1]-ext[0]+1;
120 int dimY=ext[3]-ext[2]+1;
121 int dimZ=ext[5]-ext[4]+1;
124 //CreaMaracasVisu Class which evaluates whether a point is inside the contour or not
125 ContourExtractData *extract = new ContourExtractData(false);
127 //Contour list just with the contour created with the points
128 std::vector<manualBaseModel*> list;
129 list.push_back(model);
131 //Y-value asigned arbitrary
132 extract->SetSizeImageY(500);
135 extract->SetLstManualContourModel(list);
136 extract->InitLstContoursLinesYPoints();
138 cout<<"Processing..."<<endl;
140 for (int i=0; i<dimX; i++)
143 cout<<"Print "<<i<<" of "<<dimX<<endl;
145 for (int j=0; j<dimY; j++)
147 for (int k=0; k<dimZ; k++)
149 unsigned short *pImage=(unsigned short*)_inImage->GetScalarPointer(i,j,k);
150 //unsigned short *pResult=(unsigned short*)_outImage->GetScalarPointer(i,j,k);
157 //Transform every point in the imageData to the space resultant of the mathematical transform
158 _transform->MultiplyPoint(in,out);
160 //All the imageData minus the minimum to translate to the positive octant
163 //Verify if the point is inside the contour or not
164 if ((out[1]>=0 && out[1]<500) && (extract->isInside(out[0],out[1],0)==true))
167 if(_cutInsideOutside==0)//CutInside
175 if(_cutInsideOutside==1)//CutOutside
185 _inImage->Modified();
190 //-------------------------------------------------------------------------
191 void CutModelPolygon::transformContourPoints(std::vector<double> *vectorOutX,std::vector<double> *vectorOutY,std::vector<double> *vectorOutZ)
193 int numPoints = _points->GetNumberOfPoints();
195 for(int t=0;t<numPoints;t++)
197 double point[3],in[4],out[4];
198 _points->GetPoint(t,point);
204 //CHG_LYON_JAN29 Multiply Contour Points by the inverse of the spacing. It's necessary to faire
205 //the hole in the correct position and not in other with a displacement
208 _inImage->GetSpacing(_spc);
209 in[0]=in[0]*(1/_spc[0]);
210 in[1]=in[1]*(1/_spc[1]);
211 in[2]=in[2]*(1/_spc[2]);
213 _transform->MultiplyPoint(in,out);
214 vectorOutX->push_back(out[0]);
215 vectorOutY->push_back(out[1]);
216 vectorOutZ->push_back(out[2]);
220 //-------------------------------------------------------------------------
221 void CutModelPolygon::updateTransform(double* v1, double* v2, double* v3)
224 // ImageData Centroid
225 int i, numPoints=_points->GetNumberOfPoints();
226 double centerX=0,centerY=0,centerZ=0;
227 for(i =0; i<numPoints;i++)
230 _points->GetPoint(i,point);
236 centerX=centerX/numPoints;
237 centerY=centerY/numPoints;
238 centerZ=centerZ/numPoints;
241 vtkMatrix4x4 *matrix = vtkMatrix4x4::New();
244 matrix->SetElement(0,0,v1[0]);
245 matrix->SetElement(1,0,v1[1]);
246 matrix->SetElement(2,0,v1[2]);
247 matrix->SetElement(3,0,0);
251 matrix->SetElement(0,1,v2[0]);
252 matrix->SetElement(1,1,v2[1]);
253 matrix->SetElement(2,1,v2[2]);
254 matrix->SetElement(3,1,0);
257 matrix->SetElement(0,2,v3[0]);
258 matrix->SetElement(1,2,v3[1]);
259 matrix->SetElement(2,2,v3[2]);
260 matrix->SetElement(3,2,0);
263 matrix->SetElement(0,3,centerX);
264 matrix->SetElement(1,3,centerY);
265 matrix->SetElement(2,3,centerZ);
266 matrix->SetElement(3,3,1);
268 _transform->SetMatrix(matrix);
269 _transform->Update();
271 _transform->Inverse();
272 _transform->Update();
275 //-------------------------------------------------------------------------
276 void CutModelPolygon::calculateOrthogonalVectors(double* v1, double* v2)
279 _points->GetPoint(0,p1);
280 _points->GetPoint(1,p2);
288 double magV1 = sqrt( v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2]);
293 v2[0] = v1[1]*_direction[2] - v1[2]*_direction[1];
294 v2[1] = -(v1[0]*_direction[2] - v1[2]*_direction[0]);
295 v2[2] = v1[0]*_direction[1] - v1[1]*_direction[0];
297 double magV2 = sqrt( v2[0]*v2[0] + v2[1]*v2[1] + v2[2]*v2[2]);
302 double magDir = sqrt( _direction[0]*_direction[0] + _direction[1]*_direction[1] + _direction[2]*_direction[2]);
303 _direction[0]=_direction[0]/magDir;
304 _direction[1]=_direction[1]/magDir;
305 _direction[2]=_direction[2]/magDir;
308 //-------------------------------------------------------------------------
309 void CutModelPolygon::initializeOutputImage()
312 //_inImage->GetWholeExtent(ext);
313 //int dimX=ext[1]-ext[0]+1;
314 //int dimY=ext[3]-ext[2]+1;
315 //int dimZ=ext[5]-ext[4]+1;
317 //if (_outImage != NULL ) { _outImage->Delete(); }
318 //_outImage=vtkImageData::New();
319 //_outImage->SetDimensions(dimX,dimY,dimZ);
320 //_outImage->SetScalarTypeToUnsignedShort();
321 //_outImage->AllocateScalars();
324 //-------------------------------------------------------------------------
326 manualBaseModel* CutModelPolygon::InitializeContourModel(std::vector<double> pointsX, std::vector<double> pointsY, std::vector<double> pointsZ)
328 manualBaseModel *model = new manualContourModelPolygon();
330 int numPoints = pointsX.size();
331 for(int i = 0; i<numPoints;i++)
333 model->AddPoint(pointsX[i],pointsY[i],pointsZ[i]);
341 ////////////////////////////////////////////
342 // Getters and setters
343 ////////////////////////////////////////////
345 vtkImageData* CutModelPolygon::getInImage()
350 vtkImageData* CutModelPolygon::getOutImage()
355 vtkPoints* CutModelPolygon::getPoints()
360 double* CutModelPolygon::getDirection()
365 void CutModelPolygon::setInImage(vtkImageData* pImage)
370 void CutModelPolygon::setOutImage(vtkImageData* pImage)
375 void CutModelPolygon::setPoints(vtkPoints *pPoints)
380 void CutModelPolygon::setDirection(double *pDirection)
382 _direction=pDirection;