lib/maracasVisuLib/src/CutModule/kernel/CutModelPolygon.cxx

   1
   2
   3 #include "CutModelPolygon.h"
   4
   5
   6 CutModelPolygon::CutModelPolygon()
   7 {
   8         _cutInsideOutside=0;
   9         _inImage=NULL;
  10         _transform = vtkTransform::New();
  11 }
  12
  13 CutModelPolygon::~CutModelPolygon()
  14 {
  15         if(_transform!=NULL)
  16         {
  17                 _transform->Delete();
  18         }
  19         if(_points!=NULL)
  20         {
  21                 _points->Delete();
  22         }
  23 }
  24
  25 void CutModelPolygon::processOutImage(int cutInsideOutside)
  26 {
  27         _cutInsideOutside=cutInsideOutside;
  28
  29         int numPoints = _points->GetNumberOfPoints();
  30
  31         double v1[3],v2[3];
  32
  33         // Calculate Orthogonal Vectors using two vectors in the plane and its cross product
  34         calculateOrthogonalVectors(v1,v2);
  35
  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;
  40         cout<<endl;
  41
  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);
  45
  46         //Transform Points coordinates
  47         std::vector<double> vectorOutX;
  48         std::vector<double> vectorOutY;
  49         std::vector<double> vectorOutZ;
  50
  51         transformContourPoints(&vectorOutX,&vectorOutY,&vectorOutZ);
  52
  53
  54         /// Printing Points
  55         int num = vectorOutX.size();
  56         for(int t=0;t<num;t++)
  57         {
  58                 cout<<"Final Point"<<t<<"-X:"<<vectorOutX[t]<<" Y:"<<vectorOutY[t]<<" Z:"<<vectorOutZ[t]<<endl;
  59         }
  60
  61         cutInputImage(vectorOutX,vectorOutY,vectorOutZ);
  62
  63 }
  64
  65 //-------------------------------------------------------------------------
  66 void CutModelPolygon::cutInputImage(std::vector<double> vectorOutX,std::vector<double> vectorOutY,std::vector<double> vectorOutZ)
  67 {
  68
  69         // Find the minimum value in Y axis
  70         int i;
  71         double minY=vectorOutY[0];
  72         for(i=1;i<vectorOutY.size();i++)
  73         {
  74                 if(vectorOutY[i]<minY)
  75                 {
  76                         minY=vectorOutY[i];
  77                 }
  78         }
  79
  80         // All the contour points minus the minimum to translate to the positive octant
  81         /// FIX ME !!!
  82         for(i=0;i<vectorOutY.size();i++)
  83         {
  84
  85                 vectorOutY[i]=vectorOutY[i]-minY;
  86         }
  87
  88
  89         // Creates a poligonal contour model with the points
  90         manualBaseModel *model = InitializeContourModel(vectorOutX,vectorOutY,vectorOutZ);
  91
  92         initializeOutputImage();
  93
  94         int ext[6];
  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;
  99
 100
 101         //CreaMaracasVisu Class which evaluates whether a point is inside the contour or not
 102         ContourExtractData *extract = new ContourExtractData(false);
 103
 104         //Contour list just with the contour created with the points
 105         std::vector<manualBaseModel*> list;
 106         list.push_back(model);
 107
 108         //Y-value asigned arbitrary
 109         extract->SetSizeImageY(500);
 110
 111         //Initialization
 112         extract->SetLstManualContourModel(list);
 113         extract->InitLstContoursLinesYPoints();
 114
 115         cout<<"Processing..."<<endl;
 116
 117         for (int i=0; i<dimX; i++)
 118         {
 119                 if(i%10==0)
 120                         cout<<"Print "<<i<<" of "<<dimX<<endl;
 121
 122                 for (int j=0; j<dimY; j++)
 123                 {
 124                         for (int k=0; k<dimZ; k++)
 125                         {
 126                                 unsigned short *pImage=(unsigned short*)_inImage->GetScalarPointer(i,j,k);
 127                                 //unsigned short *pResult=(unsigned short*)_outImage->GetScalarPointer(i,j,k);
 128                                 double in[4],out[4];
 129                                 in[0]=i;
 130                                 in[1]=j;
 131                                 in[2]=k;
 132                                 in[3]=1;
 133
 134                                 //Transform every point in the imageData to the space resultant of the mathematical transform
 135                                 _transform->MultiplyPoint(in,out);
 136
 137                                 //All the imageData minus the minimum to translate to the positive octant
 138                                 out[1]=out[1]- minY;
 139
 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))
 142                                 {
 143
 144                                         if(_cutInsideOutside==0)//CutInside
 145                                         {
 146                                                 *pImage=0;
 147                                         }
 148
 149                                 }
 150                                 else
 151                                 {
 152                                         if(_cutInsideOutside==1)//CutOutside
 153                                         {
 154                                                 *pImage=0;
 155                                         }
 156                                 }
 157
 158                         } // for k
 159                 } // for j
 160         }// for i
 161
 162         _inImage->Modified();
 163         _inImage->Update();
 164
 165 }
 166
 167 //-------------------------------------------------------------------------
 168 void CutModelPolygon::transformContourPoints(std::vector<double> *vectorOutX,std::vector<double> *vectorOutY,std::vector<double> *vectorOutZ)
 169 {
 170         int numPoints = _points->GetNumberOfPoints();
 171
 172         for(int t=0;t<numPoints;t++)
 173         {
 174                 double point[3],in[4],out[4];
 175                 _points->GetPoint(t,point);
 176                 in[0]=point[0];
 177                 in[1]=point[1];
 178                 in[2]=point[2];
 179                 in[3]=1;
 180
 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
 183
 184                 double _spc[3];
 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]);
 189
 190                 _transform->MultiplyPoint(in,out);
 191                 vectorOutX->push_back(out[0]);
 192                 vectorOutY->push_back(out[1]);
 193                 vectorOutZ->push_back(out[2]);
 194         }
 195
 196 }
 197 //-------------------------------------------------------------------------
 198 void CutModelPolygon::updateTransform(double* v1, double* v2, double* v3)
 199 {
 200
 201         // ImageData Centroid
 202         int i, numPoints=_points->GetNumberOfPoints();
 203         double centerX=0,centerY=0,centerZ=0;
 204         for(i =0; i<numPoints;i++)
 205         {
 206                 double point[3];
 207                 _points->GetPoint(i,point);
 208                 centerX+=point[0];
 209                 centerY+=point[1];
 210                 centerZ+=point[2];
 211         }
 212
 213         centerX=centerX/numPoints;
 214         centerY=centerY/numPoints;
 215         centerZ=centerZ/numPoints;
 216
 217
 218         vtkMatrix4x4 *matrix = vtkMatrix4x4::New();
 219
 220         //First Vector...
 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);
 225
 226
 227         //Second Vector...
 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);
 232
 233         //Third Vector...
 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);
 238
 239         //Fourth Vector...
 240         matrix->SetElement(0,3,centerX);
 241         matrix->SetElement(1,3,centerY);
 242         matrix->SetElement(2,3,centerZ);
 243         matrix->SetElement(3,3,1);
 244
 245         _transform->SetMatrix(matrix);
 246         _transform->Update();
 247
 248         _transform->Inverse();
 249         _transform->Update();
 250         matrix->Delete();
 251 }
 252 //-------------------------------------------------------------------------
 253 void CutModelPolygon::calculateOrthogonalVectors(double* v1, double* v2)
 254 {
 255         double p1[3],p2[3];
 256         _points->GetPoint(0,p1);
 257         _points->GetPoint(1,p2);
 258
 259
 260         v1[0]=p1[0]-p2[0];
 261         v1[1]=p1[1]-p2[1];
 262         v1[2]=p1[2]-p2[2];
 263
 264
 265         double magV1 = sqrt( v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2]);
 266         v1[0]=v1[0]/magV1;
 267         v1[1]=v1[1]/magV1;
 268         v1[2]=v1[2]/magV1;
 269
 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];
 273
 274         double magV2 = sqrt( v2[0]*v2[0] + v2[1]*v2[1] + v2[2]*v2[2]);
 275         v2[0]=v2[0]/magV2;
 276         v2[1]=v2[1]/magV2;
 277         v2[2]=v2[2]/magV2;
 278
 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;
 283 }
 284
 285 //-------------------------------------------------------------------------
 286 void CutModelPolygon::initializeOutputImage()
 287 {
 288         //int ext[6];
 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;
 293
 294         //if (_outImage != NULL ) { _outImage->Delete(); }
 295         //_outImage=vtkImageData::New();
 296         //_outImage->SetDimensions(dimX,dimY,dimZ);
 297         //_outImage->SetScalarTypeToUnsignedShort();
 298         //_outImage->AllocateScalars();
 299 }
 300
 301 //-------------------------------------------------------------------------
 302
 303 manualBaseModel* CutModelPolygon::InitializeContourModel(std::vector<double> pointsX, std::vector<double> pointsY,  std::vector<double> pointsZ)
 304 {
 305         manualBaseModel *model = new manualContourModelPolygon();
 306
 307         int numPoints = pointsX.size();
 308         for(int i = 0; i<numPoints;i++)
 309         {
 310                 model->AddPoint(pointsX[i],pointsY[i],pointsZ[i]);
 311         }
 312
 313         return model;
 314 }
 315
 316
 317
 318 ////////////////////////////////////////////
 319 // Getters and setters
 320 ////////////////////////////////////////////
 321
 322 vtkImageData* CutModelPolygon::getInImage()
 323 {
 324         return _inImage;
 325 }
 326
 327 vtkImageData* CutModelPolygon::getOutImage()
 328 {
 329         return _inImage;
 330 }
 331
 332 vtkPoints* CutModelPolygon::getPoints()
 333 {
 334         return _points;
 335 }
 336
 337 double* CutModelPolygon::getDirection()
 338 {
 339         return _direction;
 340 }
 341
 342 void CutModelPolygon::setInImage(vtkImageData* pImage)
 343 {
 344         _inImage=pImage;
 345 }
 346
 347 void CutModelPolygon::setOutImage(vtkImageData* pImage)
 348 {
 349         _inImage=pImage;
 350 }
 351
 352 void CutModelPolygon::setPoints(vtkPoints *pPoints)
 353 {
 354         _points=pPoints;
 355 }
 356
 357 void CutModelPolygon::setDirection(double *pDirection)
 358 {
 359         _direction=pDirection;
 360 }