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