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

   1 /*# ---------------------------------------------------------------------
   2 #
   3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
   4 #                        pour la Sant�)
   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
   8 #
   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.
  15 #
  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
  20 #  liability.
  21 #
  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 # ------------------------------------------------------------------------ */
  25
  26 #include "CutModelPolygon.h"
  27
  28
  29 CutModelPolygon::CutModelPolygon()
  30 {
  31         _cutInsideOutside=0;
  32         _inImage=NULL;
  33         _transform = vtkTransform::New();
  34 }
  35
  36 CutModelPolygon::~CutModelPolygon()
  37 {
  38         if(_transform!=NULL)
  39         {
  40                 _transform->Delete();
  41         }
  42         if(_points!=NULL)
  43         {
  44                 _points->Delete();
  45         }
  46 }
  47
  48 void CutModelPolygon::processOutImage(int cutInsideOutside)
  49 {
  50         _cutInsideOutside=cutInsideOutside;
  51
  52         int numPoints = _points->GetNumberOfPoints();
  53
  54         double v1[3],v2[3];
  55
  56         // Calculate Orthogonal Vectors using two vectors in the plane and its cross product
  57         calculateOrthogonalVectors(v1,v2);
  58
  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;
  63         cout<<endl;
  64
  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);
  68
  69         //Transform Points coordinates
  70         std::vector<double> vectorOutX;
  71         std::vector<double> vectorOutY;
  72         std::vector<double> vectorOutZ;
  73
  74         transformContourPoints(&vectorOutX,&vectorOutY,&vectorOutZ);
  75
  76
  77         /// Printing Points
  78         int num = vectorOutX.size();
  79         for(int t=0;t<num;t++)
  80         {
  81                 cout<<"Final Point"<<t<<"-X:"<<vectorOutX[t]<<" Y:"<<vectorOutY[t]<<" Z:"<<vectorOutZ[t]<<endl;
  82         }
  83
  84         cutInputImage(vectorOutX,vectorOutY,vectorOutZ);
  85
  86 }
  87
  88 //-------------------------------------------------------------------------
  89 void CutModelPolygon::cutInputImage(std::vector<double> vectorOutX,std::vector<double> vectorOutY,std::vector<double> vectorOutZ)
  90 {
  91
  92         // Find the minimum value in Y axis
  93         int i;
  94         double minY=vectorOutY[0];
  95         for(i=1;i<vectorOutY.size();i++)
  96         {
  97                 if(vectorOutY[i]<minY)
  98                 {
  99                         minY=vectorOutY[i];
 100                 }
 101         }
 102
 103         // All the contour points minus the minimum to translate to the positive octant
 104         /// FIX ME !!!
 105         for(i=0;i<vectorOutY.size();i++)
 106         {
 107
 108                 vectorOutY[i]=vectorOutY[i]-minY;
 109         }
 110
 111
 112         // Creates a poligonal contour model with the points
 113         manualBaseModel *model = InitializeContourModel(vectorOutX,vectorOutY,vectorOutZ);
 114
 115         initializeOutputImage();
 116
 117         int ext[6];
 118 //EED 2017-01-01 Migration VTK7
 119 #if VTK_MAJOR_VERSION <= 5
 120         _inImage->GetWholeExtent(ext);
 121 #else
 122         _inImage->GetExtent(ext);
 123 #endif
 124         int dimX=ext[1]-ext[0]+1;
 125         int dimY=ext[3]-ext[2]+1;
 126         int dimZ=ext[5]-ext[4]+1;
 127
 128
 129         //CreaMaracasVisu Class which evaluates whether a point is inside the contour or not
 130         ContourExtractData *extract = new ContourExtractData(false);
 131
 132         //Contour list just with the contour created with the points
 133         std::vector<manualBaseModel*> list;
 134         list.push_back(model);
 135
 136         //Y-value asigned arbitrary
 137         extract->SetSizeImageY(500);
 138
 139         //Initialization
 140         extract->SetLstManualContourModel(list);
 141         extract->InitLstContoursLinesYPoints();
 142
 143         cout<<"Processing..."<<endl;
 144
 145         for (int i=0; i<dimX; i++)
 146         {
 147                 if(i%10==0)
 148                         cout<<"Print "<<i<<" of "<<dimX<<endl;
 149
 150                 for (int j=0; j<dimY; j++)
 151                 {
 152                         for (int k=0; k<dimZ; k++)
 153                         {
 154                                 unsigned short *pImage=(unsigned short*)_inImage->GetScalarPointer(i,j,k);
 155                                 //unsigned short *pResult=(unsigned short*)_outImage->GetScalarPointer(i,j,k);
 156                                 double in[4],out[4];
 157                                 in[0]=i;
 158                                 in[1]=j;
 159                                 in[2]=k;
 160                                 in[3]=1;
 161
 162                                 //Transform every point in the imageData to the space resultant of the mathematical transform
 163                                 _transform->MultiplyPoint(in,out);
 164
 165                                 //All the imageData minus the minimum to translate to the positive octant
 166                                 out[1]=out[1]- minY;
 167
 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))
 170                                 {
 171
 172                                         if(_cutInsideOutside==0)//CutInside
 173                                         {
 174                                                 *pImage=0;
 175                                         }
 176
 177                                 }
 178                                 else
 179                                 {
 180                                         if(_cutInsideOutside==1)//CutOutside
 181                                         {
 182                                                 *pImage=0;
 183                                         }
 184                                 }
 185
 186                         } // for k
 187                 } // for j
 188         }// for i
 189
 190         _inImage->Modified();
 191 //EED 2017-01-01 Migration VTK7
 192 #if VTK_MAJOR_VERSION <= 5
 193         _inImage->Update();
 194 #else
 195         // ..
 196 #endif
 197
 198 }
 199
 200 //-------------------------------------------------------------------------
 201 void CutModelPolygon::transformContourPoints(std::vector<double> *vectorOutX,std::vector<double> *vectorOutY,std::vector<double> *vectorOutZ)
 202 {
 203         int numPoints = _points->GetNumberOfPoints();
 204
 205         for(int t=0;t<numPoints;t++)
 206         {
 207                 double point[3],in[4],out[4];
 208                 _points->GetPoint(t,point);
 209                 in[0]=point[0];
 210                 in[1]=point[1];
 211                 in[2]=point[2];
 212                 in[3]=1;
 213
 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
 216
 217                 double _spc[3];
 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]);
 222
 223                 _transform->MultiplyPoint(in,out);
 224                 vectorOutX->push_back(out[0]);
 225                 vectorOutY->push_back(out[1]);
 226                 vectorOutZ->push_back(out[2]);
 227         }
 228
 229 }
 230 //-------------------------------------------------------------------------
 231 void CutModelPolygon::updateTransform(double* v1, double* v2, double* v3)
 232 {
 233
 234         // ImageData Centroid
 235         int i, numPoints=_points->GetNumberOfPoints();
 236         double centerX=0,centerY=0,centerZ=0;
 237         for(i =0; i<numPoints;i++)
 238         {
 239                 double point[3];
 240                 _points->GetPoint(i,point);
 241                 centerX+=point[0];
 242                 centerY+=point[1];
 243                 centerZ+=point[2];
 244         }
 245
 246         centerX=centerX/numPoints;
 247         centerY=centerY/numPoints;
 248         centerZ=centerZ/numPoints;
 249
 250
 251         vtkMatrix4x4 *matrix = vtkMatrix4x4::New();
 252
 253         //First Vector...
 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);
 258
 259
 260         //Second Vector...
 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);
 265
 266         //Third Vector...
 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);
 271
 272         //Fourth Vector...
 273         matrix->SetElement(0,3,centerX);
 274         matrix->SetElement(1,3,centerY);
 275         matrix->SetElement(2,3,centerZ);
 276         matrix->SetElement(3,3,1);
 277
 278         _transform->SetMatrix(matrix);
 279         _transform->Update();
 280
 281         _transform->Inverse();
 282         _transform->Update();
 283         matrix->Delete();
 284 }
 285 //-------------------------------------------------------------------------
 286 void CutModelPolygon::calculateOrthogonalVectors(double* v1, double* v2)
 287 {
 288         double p1[3],p2[3];
 289         _points->GetPoint(0,p1);
 290         _points->GetPoint(1,p2);
 291
 292
 293         v1[0]=p1[0]-p2[0];
 294         v1[1]=p1[1]-p2[1];
 295         v1[2]=p1[2]-p2[2];
 296
 297
 298         double magV1 = sqrt( v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2]);
 299         v1[0]=v1[0]/magV1;
 300         v1[1]=v1[1]/magV1;
 301         v1[2]=v1[2]/magV1;
 302
 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];
 306
 307         double magV2 = sqrt( v2[0]*v2[0] + v2[1]*v2[1] + v2[2]*v2[2]);
 308         v2[0]=v2[0]/magV2;
 309         v2[1]=v2[1]/magV2;
 310         v2[2]=v2[2]/magV2;
 311
 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;
 316 }
 317
 318 //-------------------------------------------------------------------------
 319 void CutModelPolygon::initializeOutputImage()
 320 {
 321         //int ext[6];
 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;
 326
 327         //if (_outImage != NULL ) { _outImage->Delete(); }
 328         //_outImage=vtkImageData::New();
 329         //_outImage->SetDimensions(dimX,dimY,dimZ);
 330         //_outImage->SetScalarTypeToUnsignedShort();
 331         //_outImage->AllocateScalars();
 332 }
 333
 334 //-------------------------------------------------------------------------
 335
 336 manualBaseModel* CutModelPolygon::InitializeContourModel(std::vector<double> pointsX, std::vector<double> pointsY,  std::vector<double> pointsZ)
 337 {
 338         manualBaseModel *model = new manualContourModelPolygon();
 339
 340         int numPoints = pointsX.size();
 341         for(int i = 0; i<numPoints;i++)
 342         {
 343                 model->AddPoint(pointsX[i],pointsY[i],pointsZ[i]);
 344         }
 345
 346         return model;
 347 }
 348
 349
 350
 351 ////////////////////////////////////////////
 352 // Getters and setters
 353 ////////////////////////////////////////////
 354
 355 vtkImageData* CutModelPolygon::getInImage()
 356 {
 357         return _inImage;
 358 }
 359
 360 vtkImageData* CutModelPolygon::getOutImage()
 361 {
 362         return _inImage;
 363 }
 364
 365 vtkPoints* CutModelPolygon::getPoints()
 366 {
 367         return _points;
 368 }
 369
 370 double* CutModelPolygon::getDirection()
 371 {
 372         return _direction;
 373 }
 374
 375 void CutModelPolygon::setInImage(vtkImageData* pImage)
 376 {
 377         _inImage=pImage;
 378 }
 379
 380 void CutModelPolygon::setOutImage(vtkImageData* pImage)
 381 {
 382         _inImage=pImage;
 383 }
 384
 385 void CutModelPolygon::setPoints(vtkPoints *pPoints)
 386 {
 387         _points=pPoints;
 388 }
 389
 390 void CutModelPolygon::setDirection(double *pDirection)
 391 {
 392         _direction=pDirection;
 393 }