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         _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;
 122
 123
 124         //CreaMaracasVisu Class which evaluates whether a point is inside the contour or not
 125         ContourExtractData *extract = new ContourExtractData(false);
 126
 127         //Contour list just with the contour created with the points
 128         std::vector<manualBaseModel*> list;
 129         list.push_back(model);
 130
 131         //Y-value asigned arbitrary
 132         extract->SetSizeImageY(500);
 133
 134         //Initialization
 135         extract->SetLstManualContourModel(list);
 136         extract->InitLstContoursLinesYPoints();
 137
 138         cout<<"Processing..."<<endl;
 139
 140         for (int i=0; i<dimX; i++)
 141         {
 142                 if(i%10==0)
 143                         cout<<"Print "<<i<<" of "<<dimX<<endl;
 144
 145                 for (int j=0; j<dimY; j++)
 146                 {
 147                         for (int k=0; k<dimZ; k++)
 148                         {
 149                                 unsigned short *pImage=(unsigned short*)_inImage->GetScalarPointer(i,j,k);
 150                                 //unsigned short *pResult=(unsigned short*)_outImage->GetScalarPointer(i,j,k);
 151                                 double in[4],out[4];
 152                                 in[0]=i;
 153                                 in[1]=j;
 154                                 in[2]=k;
 155                                 in[3]=1;
 156
 157                                 //Transform every point in the imageData to the space resultant of the mathematical transform
 158                                 _transform->MultiplyPoint(in,out);
 159
 160                                 //All the imageData minus the minimum to translate to the positive octant
 161                                 out[1]=out[1]- minY;
 162
 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))
 165                                 {
 166
 167                                         if(_cutInsideOutside==0)//CutInside
 168                                         {
 169                                                 *pImage=0;
 170                                         }
 171
 172                                 }
 173                                 else
 174                                 {
 175                                         if(_cutInsideOutside==1)//CutOutside
 176                                         {
 177                                                 *pImage=0;
 178                                         }
 179                                 }
 180
 181                         } // for k
 182                 } // for j
 183         }// for i
 184
 185         _inImage->Modified();
 186         _inImage->Update();
 187
 188 }
 189
 190 //-------------------------------------------------------------------------
 191 void CutModelPolygon::transformContourPoints(std::vector<double> *vectorOutX,std::vector<double> *vectorOutY,std::vector<double> *vectorOutZ)
 192 {
 193         int numPoints = _points->GetNumberOfPoints();
 194
 195         for(int t=0;t<numPoints;t++)
 196         {
 197                 double point[3],in[4],out[4];
 198                 _points->GetPoint(t,point);
 199                 in[0]=point[0];
 200                 in[1]=point[1];
 201                 in[2]=point[2];
 202                 in[3]=1;
 203
 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
 206
 207                 double _spc[3];
 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]);
 212
 213                 _transform->MultiplyPoint(in,out);
 214                 vectorOutX->push_back(out[0]);
 215                 vectorOutY->push_back(out[1]);
 216                 vectorOutZ->push_back(out[2]);
 217         }
 218
 219 }
 220 //-------------------------------------------------------------------------
 221 void CutModelPolygon::updateTransform(double* v1, double* v2, double* v3)
 222 {
 223
 224         // ImageData Centroid
 225         int i, numPoints=_points->GetNumberOfPoints();
 226         double centerX=0,centerY=0,centerZ=0;
 227         for(i =0; i<numPoints;i++)
 228         {
 229                 double point[3];
 230                 _points->GetPoint(i,point);
 231                 centerX+=point[0];
 232                 centerY+=point[1];
 233                 centerZ+=point[2];
 234         }
 235
 236         centerX=centerX/numPoints;
 237         centerY=centerY/numPoints;
 238         centerZ=centerZ/numPoints;
 239
 240
 241         vtkMatrix4x4 *matrix = vtkMatrix4x4::New();
 242
 243         //First Vector...
 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);
 248
 249
 250         //Second Vector...
 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);
 255
 256         //Third Vector...
 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);
 261
 262         //Fourth Vector...
 263         matrix->SetElement(0,3,centerX);
 264         matrix->SetElement(1,3,centerY);
 265         matrix->SetElement(2,3,centerZ);
 266         matrix->SetElement(3,3,1);
 267
 268         _transform->SetMatrix(matrix);
 269         _transform->Update();
 270
 271         _transform->Inverse();
 272         _transform->Update();
 273         matrix->Delete();
 274 }
 275 //-------------------------------------------------------------------------
 276 void CutModelPolygon::calculateOrthogonalVectors(double* v1, double* v2)
 277 {
 278         double p1[3],p2[3];
 279         _points->GetPoint(0,p1);
 280         _points->GetPoint(1,p2);
 281
 282
 283         v1[0]=p1[0]-p2[0];
 284         v1[1]=p1[1]-p2[1];
 285         v1[2]=p1[2]-p2[2];
 286
 287
 288         double magV1 = sqrt( v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2]);
 289         v1[0]=v1[0]/magV1;
 290         v1[1]=v1[1]/magV1;
 291         v1[2]=v1[2]/magV1;
 292
 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];
 296
 297         double magV2 = sqrt( v2[0]*v2[0] + v2[1]*v2[1] + v2[2]*v2[2]);
 298         v2[0]=v2[0]/magV2;
 299         v2[1]=v2[1]/magV2;
 300         v2[2]=v2[2]/magV2;
 301
 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;
 306 }
 307
 308 //-------------------------------------------------------------------------
 309 void CutModelPolygon::initializeOutputImage()
 310 {
 311         //int ext[6];
 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;
 316
 317         //if (_outImage != NULL ) { _outImage->Delete(); }
 318         //_outImage=vtkImageData::New();
 319         //_outImage->SetDimensions(dimX,dimY,dimZ);
 320         //_outImage->SetScalarTypeToUnsignedShort();
 321         //_outImage->AllocateScalars();
 322 }
 323
 324 //-------------------------------------------------------------------------
 325
 326 manualBaseModel* CutModelPolygon::InitializeContourModel(std::vector<double> pointsX, std::vector<double> pointsY,  std::vector<double> pointsZ)
 327 {
 328         manualBaseModel *model = new manualContourModelPolygon();
 329
 330         int numPoints = pointsX.size();
 331         for(int i = 0; i<numPoints;i++)
 332         {
 333                 model->AddPoint(pointsX[i],pointsY[i],pointsZ[i]);
 334         }
 335
 336         return model;
 337 }
 338
 339
 340
 341 ////////////////////////////////////////////
 342 // Getters and setters
 343 ////////////////////////////////////////////
 344
 345 vtkImageData* CutModelPolygon::getInImage()
 346 {
 347         return _inImage;
 348 }
 349
 350 vtkImageData* CutModelPolygon::getOutImage()
 351 {
 352         return _inImage;
 353 }
 354
 355 vtkPoints* CutModelPolygon::getPoints()
 356 {
 357         return _points;
 358 }
 359
 360 double* CutModelPolygon::getDirection()
 361 {
 362         return _direction;
 363 }
 364
 365 void CutModelPolygon::setInImage(vtkImageData* pImage)
 366 {
 367         _inImage=pImage;
 368 }
 369
 370 void CutModelPolygon::setOutImage(vtkImageData* pImage)
 371 {
 372         _inImage=pImage;
 373 }
 374
 375 void CutModelPolygon::setPoints(vtkPoints *pPoints)
 376 {
 377         _points=pPoints;
 378 }
 379
 380 void CutModelPolygon::setDirection(double *pDirection)
 381 {
 382         _direction=pDirection;
 383 }