common/clitkDicomRTStruct2ImageFilter.cxx

   1 /*=========================================================================
   2   Program:         vv http://www.creatis.insa-lyon.fr/rio/vv
   3   Main authors :   XX XX XX
   4
   5   Authors belongs to:
   6   - University of LYON           http://www.universite-lyon.fr/
   7   - Léon Bérard cancer center    http://www.centreleonberard.fr
   8   - CREATIS CNRS laboratory      http://www.creatis.insa-lyon.fr
   9
  10   This software is distributed WITHOUT ANY WARRANTY; without even
  11   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  12   PURPOSE.  See the copyright notices for more information.
  13
  14   It is distributed under dual licence
  15   - BSD       http://www.opensource.org/licenses/bsd-license.php
  16   - CeCILL-B  http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17
  18   =========================================================================*/
  19
  20 #include <iterator>
  21 #include <algorithm>
  22
  23 // clitk
  24 #include "clitkDicomRTStruct2ImageFilter.h"
  25 #include "clitkImageCommon.h"
  26
  27 // vtk
  28 #include <vtkVersion.h>
  29 #include <vtkPolyDataToImageStencil.h>
  30 #include <vtkSmartPointer.h>
  31 #include <vtkImageStencil.h>
  32 #include <vtkLinearExtrusionFilter.h>
  33 #include <vtkMetaImageWriter.h>
  34
  35
  36 //--------------------------------------------------------------------
  37 clitk::DicomRTStruct2ImageFilter::DicomRTStruct2ImageFilter()
  38 {
  39   mROI = NULL;
  40   mWriteOutput = false;
  41   mCropMask = true;
  42 }
  43 //--------------------------------------------------------------------
  44
  45
  46 //--------------------------------------------------------------------
  47 clitk::DicomRTStruct2ImageFilter::~DicomRTStruct2ImageFilter()
  48 {
  49
  50 }
  51 //--------------------------------------------------------------------
  52
  53
  54 //--------------------------------------------------------------------
  55 bool clitk::DicomRTStruct2ImageFilter::ImageInfoIsSet() const
  56 {
  57   return mSize.size() && mSpacing.size() && mOrigin.size();
  58 }
  59 //--------------------------------------------------------------------
  60
  61
  62 //--------------------------------------------------------------------
  63 void clitk::DicomRTStruct2ImageFilter::SetWriteOutputFlag(bool b)
  64 {
  65   mWriteOutput = b;
  66 }
  67 //--------------------------------------------------------------------
  68
  69
  70 //--------------------------------------------------------------------
  71 void clitk::DicomRTStruct2ImageFilter::SetROI(clitk::DicomRT_ROI * roi)
  72 {
  73   mROI = roi;
  74 }
  75 //--------------------------------------------------------------------
  76
  77
  78 //--------------------------------------------------------------------
  79 void clitk::DicomRTStruct2ImageFilter::SetCropMaskEnabled(bool b)
  80 {
  81   mCropMask = b;
  82 }
  83 //--------------------------------------------------------------------
  84
  85
  86 //--------------------------------------------------------------------
  87 void clitk::DicomRTStruct2ImageFilter::SetOutputImageFilename(std::string s)
  88 {
  89   mOutputFilename = s;
  90   mWriteOutput = true;
  91 }
  92 //--------------------------------------------------------------------
  93
  94
  95 //--------------------------------------------------------------------
  96 void clitk::DicomRTStruct2ImageFilter::SetImage(vvImage::Pointer image)
  97 {
  98   if (image->GetNumberOfDimensions() != 3) {
  99     std::cerr << "Error. Please provide a 3D image." << std::endl;
 100     exit(0);
 101   }
 102   mSpacing.resize(3);
 103   mOrigin.resize(3);
 104   mSize.resize(3);
 105   mDirection.resize(3);
 106   mTransformMatrix = image->GetTransform()[0]->GetMatrix();
 107   for(unsigned int i=0; i<3; i++) {
 108     mSpacing[i] = image->GetSpacing()[i];
 109     mOrigin[i] = image->GetOrigin()[i];
 110     mSize[i] = image->GetSize()[i];
 111     mDirection[i].resize(3);
 112     for(unsigned int j=0; j<3; j++)
 113       mDirection[i][j] = image->GetDirection()[i][j];
 114   }
 115 }
 116 //--------------------------------------------------------------------
 117
 118
 119 //--------------------------------------------------------------------
 120 void clitk::DicomRTStruct2ImageFilter::SetImageFilename(std::string f)
 121 {
 122   itk::ImageIOBase::Pointer header = clitk::readImageHeader(f);
 123   if (header->GetNumberOfDimensions() < 3) {
 124     std::cerr << "Error. Please provide a 3D image instead of " << f << std::endl;
 125     exit(0);
 126   }
 127   if (header->GetNumberOfDimensions() > 3) {
 128     std::cerr << "Warning dimension > 3 are ignored" << std::endl;
 129   }
 130   mSpacing.resize(3);
 131   mOrigin.resize(3);
 132   mSize.resize(3);
 133   mDirection.resize(3);
 134   for(unsigned int i=0; i<3; i++) {
 135     mSpacing[i] = header->GetSpacing(i);
 136     mOrigin[i] = header->GetOrigin(i);
 137     mSize[i] = header->GetDimensions(i);
 138     mDirection[i].resize(3);
 139     for(unsigned int j=0; j<3; j++)
 140       mDirection[i][j] = header->GetDirection(i)[j];
 141   }
 142 }
 143 //--------------------------------------------------------------------
 144
 145
 146 //--------------------------------------------------------------------
 147 void clitk::DicomRTStruct2ImageFilter::SetOutputOrigin(const double* origin)
 148 {
 149   std::copy(origin,origin+3,std::back_inserter(mOrigin));
 150 }
 151 //--------------------------------------------------------------------
 152
 153
 154 //--------------------------------------------------------------------
 155 void clitk::DicomRTStruct2ImageFilter::SetOutputSpacing(const double* spacing)
 156 {
 157   std::copy(spacing,spacing+3,std::back_inserter(mSpacing));
 158 }
 159 //--------------------------------------------------------------------
 160
 161
 162 //--------------------------------------------------------------------
 163 void clitk::DicomRTStruct2ImageFilter::SetOutputSize(const unsigned long* size)
 164 {
 165   std::copy(size,size+3,std::back_inserter(mSize));
 166 }
 167 //--------------------------------------------------------------------
 168
 169
 170 //--------------------------------------------------------------------
 171 void clitk::DicomRTStruct2ImageFilter::Update()
 172 {
 173   if (!mROI) {
 174     std::cerr << "Error. No ROI set, please use SetROI." << std::endl;
 175     exit(0);
 176   }
 177   if (!ImageInfoIsSet()) {
 178     std::cerr << "Error. Please provide image info (spacing/origin) with SetImageFilename" << std::endl;
 179     exit(0);
 180   }
 181
 182   // Get Mesh
 183   vtkPolyData * mesh = mROI->GetMesh();
 184
 185   // Get bounds
 186   double *bounds=mesh->GetBounds();
 187
 188   //Change mOrigin, mSize and mSpacing with respect to the directions
 189   // Spacing is influenced by input direction
 190   std::vector<double> tempSpacing;
 191   tempSpacing.resize(3);
 192   for(int i=0; i<3; i++) {
 193     tempSpacing[i] = 0.0;
 194     for(int j=0; j<3; j++) {
 195       tempSpacing[i] += mDirection[i][j] * mSpacing[j];
 196     }
 197   }
 198   for(int i=0; i<3; i++)
 199     mSpacing[i] = tempSpacing[i];
 200
 201   // Size is influenced by affine transform matrix and input direction
 202   // Size is converted to double, transformed and converted back to size type.
 203   std::vector<double> tempSize;
 204   tempSize.resize(3);
 205   for(int i=0; i<3; i++) {
 206     tempSize[i] = 0.0;
 207     for(int j=0; j<3; j++) {
 208       tempSize[i] += mDirection[i][j] * mSize[j];
 209     }
 210   }
 211   for(int i=0; i<3; i++) {
 212     if (tempSize[i] < 0.0) {
 213       tempSize[i] *= -1;
 214       mOrigin[i] += mSpacing[i]*(tempSize[i] -1);
 215       mSpacing[i] *= -1;
 216     }
 217     mSize[i] = lrint(tempSize[i]);
 218   }
 219
 220   // Compute origin
 221   std::vector<double> origin;
 222   origin.resize(3);
 223   origin[0] = floor((bounds[0]-mOrigin[0])/mSpacing[0]-2)*mSpacing[0]+mOrigin[0];
 224   origin[1] = floor((bounds[2]-mOrigin[1])/mSpacing[1]-2)*mSpacing[1]+mOrigin[1];
 225   origin[2] = floor((bounds[4]-mOrigin[2])/mSpacing[2]-2)*mSpacing[2]+mOrigin[2];
 226
 227   // Compute extend
 228   std::vector<double> extend;
 229   extend.resize(3);
 230   extend[0] = ceil((bounds[1]-origin[0])/mSpacing[0]+4);
 231   extend[1] = ceil((bounds[3]-origin[1])/mSpacing[1]+4);
 232   extend[2] = ceil((bounds[5]-origin[2])/mSpacing[2]+4);
 233   // If no crop, set initial image size/origin
 234   if (!mCropMask) {
 235     for(int i=0; i<3; i++) {
 236       origin[i] = mOrigin[i];
 237       extend[i] = mSize[i]-1;
 238     }
 239   }
 240
 241 std::cout << "origin " << origin[0] << " " << origin[1] << " " << origin[2] << std::endl;
 242 std::cout << "extend " << extend[0] << " " << extend[1] << " " << extend[2] << std::endl;
 243   // Create new output image
 244   mBinaryImage = vtkSmartPointer<vtkImageData>::New();
 245 #if VTK_MAJOR_VERSION <= 5
 246   mBinaryImage->SetScalarTypeToUnsignedChar();
 247 #endif
 248   mBinaryImage->SetOrigin(&origin[0]);
 249   mBinaryImage->SetSpacing(&mSpacing[0]);
 250   mBinaryImage->SetExtent(0, extend[0],
 251                           0, extend[1],
 252                           0, extend[2]);
 253 #if VTK_MAJOR_VERSION <= 5
 254   mBinaryImage->AllocateScalars();
 255 #else
 256   mBinaryImage->AllocateScalars(VTK_UNSIGNED_CHAR, 1);
 257 #endif
 258
 259   memset(mBinaryImage->GetScalarPointer(), 0,
 260          mBinaryImage->GetDimensions()[0]*mBinaryImage->GetDimensions()[1]*mBinaryImage->GetDimensions()[2]*sizeof(unsigned char));
 261
 262   // Extrude
 263   vtkSmartPointer<vtkLinearExtrusionFilter> extrude=vtkSmartPointer<vtkLinearExtrusionFilter>::New();
 264 #if VTK_MAJOR_VERSION <= 5
 265   extrude->SetInput(mesh);
 266 #else
 267   extrude->SetInputData(mesh);
 268 #endif
 269   ///We extrude in the -slice_spacing direction to respect the FOCAL convention (NEEDED !)
 270   extrude->SetVector(0, 0, -mSpacing[2]);
 271
 272   // Binarization
 273   vtkSmartPointer<vtkPolyDataToImageStencil> sts=vtkSmartPointer<vtkPolyDataToImageStencil>::New();
 274   //The following line is extremely important
 275   //http://www.nabble.com/Bug-in-vtkPolyDataToImageStencil--td23368312.html#a23370933
 276   sts->SetTolerance(0);
 277   sts->SetInformationInput(mBinaryImage);
 278 #if VTK_MAJOR_VERSION <= 5
 279   sts->SetInput(extrude->GetOutput());
 280 #else
 281   sts->SetInputConnection(extrude->GetOutputPort(0));
 282 #endif
 283   //sts->SetInput(mesh);
 284
 285   vtkSmartPointer<vtkImageStencil> stencil=vtkSmartPointer<vtkImageStencil>::New();
 286 #if VTK_MAJOR_VERSION <= 5
 287   stencil->SetStencil(sts->GetOutput());
 288 #else
 289   stencil->SetStencilConnection(sts->GetOutputPort(0));
 290 #endif
 291 #if VTK_MAJOR_VERSION <= 5
 292   stencil->SetInput(mBinaryImage);
 293 #else
 294   stencil->SetInputData(mBinaryImage);
 295 #endif
 296   stencil->ReverseStencilOn();
 297   stencil->Update();
 298
 299   /*
 300   vtkSmartPointer<vtkMetaImageWriter> w = vtkSmartPointer<vtkMetaImageWriter>::New();
 301   w->SetInput(stencil->GetOutput());
 302   w->SetFileName("binary2.mhd");
 303   w->Write();
 304   */
 305
 306   mBinaryImage->ShallowCopy(stencil->GetOutput());
 307
 308   if (mWriteOutput) {
 309     typedef itk::Image<unsigned char, 3> ImageType;
 310     typedef itk::VTKImageToImageFilter<ImageType> ConnectorType;
 311     ConnectorType::Pointer connector = ConnectorType::New();
 312     connector->SetInput(GetOutput());
 313     connector->Update();
 314     clitk::writeImage<ImageType>(connector->GetOutput(), mOutputFilename);
 315   }
 316 }
 317 //--------------------------------------------------------------------
 318
 319
 320
 321 //--------------------------------------------------------------------
 322 vtkImageData * clitk::DicomRTStruct2ImageFilter::GetOutput()
 323 {
 324   assert(mBinaryImage);
 325   return mBinaryImage;
 326 }
 327 //--------------------------------------------------------------------
 328