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