3 #include "CutModelPolygon.h"
6 CutModelPolygon::CutModelPolygon()
10 _transform = vtkTransform::New();
13 CutModelPolygon::~CutModelPolygon()
25 void CutModelPolygon::processOutImage(int cutInsideOutside)
27 _cutInsideOutside=cutInsideOutside;
29 int numPoints = _points->GetNumberOfPoints();
33 // Calculate Orthogonal Vectors using two vectors in the plane and its cross product
34 calculateOrthogonalVectors(v1,v2);
36 cout<<"PolyCutter::processOutImage"<<endl;
37 cout<<"V1:"<<" X:"<<v1[0]<<" Y:"<<v1[1]<<" Z:"<<v1[2]<<endl;
38 cout<<"V2:"<<" X:"<<v2[0]<<" Y:"<<v2[1]<<" Z:"<<v2[2]<<endl;
39 cout<<"Direction:"<<" X:"<<_direction[0]<<" Y:"<<_direction[1]<<" Z:"<<_direction[2]<<endl;
42 // Update the inverse transform which it's applied to all the points of the contour
43 // and the points of the input image.
44 updateTransform(v1,v2,_direction);
46 //Transform Points coordinates
47 std::vector<double> vectorOutX;
48 std::vector<double> vectorOutY;
49 std::vector<double> vectorOutZ;
51 transformContourPoints(&vectorOutX,&vectorOutY,&vectorOutZ);
55 int num = vectorOutX.size();
56 for(int t=0;t<num;t++)
58 cout<<"Final Point"<<t<<"-X:"<<vectorOutX[t]<<" Y:"<<vectorOutY[t]<<" Z:"<<vectorOutZ[t]<<endl;
61 cutInputImage(vectorOutX,vectorOutY,vectorOutZ);
65 //-------------------------------------------------------------------------
66 void CutModelPolygon::cutInputImage(std::vector<double> vectorOutX,std::vector<double> vectorOutY,std::vector<double> vectorOutZ)
69 // Find the minimum value in Y axis
71 double minY=vectorOutY[0];
72 for(i=1;i<vectorOutY.size();i++)
74 if(vectorOutY[i]<minY)
80 // All the contour points minus the minimum to translate to the positive octant
82 for(i=0;i<vectorOutY.size();i++)
85 vectorOutY[i]=vectorOutY[i]-minY;
89 // Creates a poligonal contour model with the points
90 manualBaseModel *model = InitializeContourModel(vectorOutX,vectorOutY,vectorOutZ);
92 initializeOutputImage();
95 _inImage->GetWholeExtent(ext);
96 int dimX=ext[1]-ext[0]+1;
97 int dimY=ext[3]-ext[2]+1;
98 int dimZ=ext[5]-ext[4]+1;
101 //CreaMaracasVisu Class which evaluates whether a point is inside the contour or not
102 ContourExtractData *extract = new ContourExtractData(false);
104 //Contour list just with the contour created with the points
105 std::vector<manualBaseModel*> list;
106 list.push_back(model);
108 //Y-value asigned arbitrary
109 extract->SetSizeImageY(500);
112 extract->SetLstManualContourModel(list);
113 extract->InitLstContoursLinesYPoints();
115 cout<<"Processing..."<<endl;
117 for (int i=0; i<dimX; i++)
120 cout<<"Print "<<i<<" of "<<dimX<<endl;
122 for (int j=0; j<dimY; j++)
124 for (int k=0; k<dimZ; k++)
126 unsigned short *pImage=(unsigned short*)_inImage->GetScalarPointer(i,j,k);
127 //unsigned short *pResult=(unsigned short*)_outImage->GetScalarPointer(i,j,k);
134 //Transform every point in the imageData to the space resultant of the mathematical transform
135 _transform->MultiplyPoint(in,out);
137 //All the imageData minus the minimum to translate to the positive octant
140 //Verify if the point is inside the contour or not
141 if ((out[1]>=0 && out[1]<500) && (extract->isInside(out[0],out[1],0)==true))
144 if(_cutInsideOutside==0)//CutInside
152 if(_cutInsideOutside==1)//CutOutside
162 _inImage->Modified();
167 //-------------------------------------------------------------------------
168 void CutModelPolygon::transformContourPoints(std::vector<double> *vectorOutX,std::vector<double> *vectorOutY,std::vector<double> *vectorOutZ)
170 int numPoints = _points->GetNumberOfPoints();
172 for(int t=0;t<numPoints;t++)
174 double point[3],in[4],out[4];
175 _points->GetPoint(t,point);
181 //CHG_LYON_JAN29 Multiply Contour Points by the inverse of the spacing. It's necessary to faire
182 //the hole in the correct position and not in other with a displacement
185 _inImage->GetSpacing(_spc);
186 in[0]=in[0]*(1/_spc[0]);
187 in[1]=in[1]*(1/_spc[1]);
188 in[2]=in[2]*(1/_spc[2]);
190 _transform->MultiplyPoint(in,out);
191 vectorOutX->push_back(out[0]);
192 vectorOutY->push_back(out[1]);
193 vectorOutZ->push_back(out[2]);
197 //-------------------------------------------------------------------------
198 void CutModelPolygon::updateTransform(double* v1, double* v2, double* v3)
201 // ImageData Centroid
202 int i, numPoints=_points->GetNumberOfPoints();
203 double centerX=0,centerY=0,centerZ=0;
204 for(i =0; i<numPoints;i++)
207 _points->GetPoint(i,point);
213 centerX=centerX/numPoints;
214 centerY=centerY/numPoints;
215 centerZ=centerZ/numPoints;
218 vtkMatrix4x4 *matrix = vtkMatrix4x4::New();
221 matrix->SetElement(0,0,v1[0]);
222 matrix->SetElement(1,0,v1[1]);
223 matrix->SetElement(2,0,v1[2]);
224 matrix->SetElement(3,0,0);
228 matrix->SetElement(0,1,v2[0]);
229 matrix->SetElement(1,1,v2[1]);
230 matrix->SetElement(2,1,v2[2]);
231 matrix->SetElement(3,1,0);
234 matrix->SetElement(0,2,v3[0]);
235 matrix->SetElement(1,2,v3[1]);
236 matrix->SetElement(2,2,v3[2]);
237 matrix->SetElement(3,2,0);
240 matrix->SetElement(0,3,centerX);
241 matrix->SetElement(1,3,centerY);
242 matrix->SetElement(2,3,centerZ);
243 matrix->SetElement(3,3,1);
245 _transform->SetMatrix(matrix);
246 _transform->Update();
248 _transform->Inverse();
249 _transform->Update();
252 //-------------------------------------------------------------------------
253 void CutModelPolygon::calculateOrthogonalVectors(double* v1, double* v2)
256 _points->GetPoint(0,p1);
257 _points->GetPoint(1,p2);
265 double magV1 = sqrt( v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2]);
270 v2[0] = v1[1]*_direction[2] - v1[2]*_direction[1];
271 v2[1] = -(v1[0]*_direction[2] - v1[2]*_direction[0]);
272 v2[2] = v1[0]*_direction[1] - v1[1]*_direction[0];
274 double magV2 = sqrt( v2[0]*v2[0] + v2[1]*v2[1] + v2[2]*v2[2]);
279 double magDir = sqrt( _direction[0]*_direction[0] + _direction[1]*_direction[1] + _direction[2]*_direction[2]);
280 _direction[0]=_direction[0]/magDir;
281 _direction[1]=_direction[1]/magDir;
282 _direction[2]=_direction[2]/magDir;
285 //-------------------------------------------------------------------------
286 void CutModelPolygon::initializeOutputImage()
289 //_inImage->GetWholeExtent(ext);
290 //int dimX=ext[1]-ext[0]+1;
291 //int dimY=ext[3]-ext[2]+1;
292 //int dimZ=ext[5]-ext[4]+1;
294 //if (_outImage != NULL ) { _outImage->Delete(); }
295 //_outImage=vtkImageData::New();
296 //_outImage->SetDimensions(dimX,dimY,dimZ);
297 //_outImage->SetScalarTypeToUnsignedShort();
298 //_outImage->AllocateScalars();
301 //-------------------------------------------------------------------------
303 manualBaseModel* CutModelPolygon::InitializeContourModel(std::vector<double> pointsX, std::vector<double> pointsY, std::vector<double> pointsZ)
305 manualBaseModel *model = new manualContourModelPolygon();
307 int numPoints = pointsX.size();
308 for(int i = 0; i<numPoints;i++)
310 model->AddPoint(pointsX[i],pointsY[i],pointsZ[i]);
318 ////////////////////////////////////////////
319 // Getters and setters
320 ////////////////////////////////////////////
322 vtkImageData* CutModelPolygon::getInImage()
327 vtkImageData* CutModelPolygon::getOutImage()
332 vtkPoints* CutModelPolygon::getPoints()
337 double* CutModelPolygon::getDirection()
342 void CutModelPolygon::setInImage(vtkImageData* pImage)
347 void CutModelPolygon::setOutImage(vtkImageData* pImage)
352 void CutModelPolygon::setPoints(vtkPoints *pPoints)
357 void CutModelPolygon::setDirection(double *pDirection)
359 _direction=pDirection;