common/clitkXml2DicomRTStructFilter.txx

   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 // clitk
  21 #include "clitkXml2DicomRTStructFilter.h"
  22
  23 // xml parser
  24 #include "../utilities/pugixml/pugixml.hpp"
  25
  26 // itk
  27 #include "itkImage.h"
  28 #include "itkImageFileReader.h"
  29
  30 // vtk
  31 #include <vtkSmartPointer.h>
  32 #include <vtkMetaImageWriter.h>
  33 #include <vtkImageData.h>
  34 #include <vtkPolygon.h>
  35 #include <vtkAppendPolyData.h>
  36 #include <vtkPointData.h>
  37
  38 // FIXME to remove
  39 #include "vtkPolyDataMapper.h"
  40 #include "vtkPolyDataMapper2D.h"
  41 #include "vtkRenderWindowInteractor.h"
  42 #include "vtkPolyDataReader.h"
  43 #include "vtkRenderWindow.h"
  44 #include "vtkRenderer.h"
  45 #include "vtkCamera.h"
  46 #include "vtkProperty.h"
  47 #include "vtkProperty2D.h"
  48 #include <vtksys/SystemTools.hxx>
  49
  50 // gdcm
  51 #include <vtkRTStructSetProperties.h>
  52 #include <vtkGDCMPolyDataReader.h>
  53 #include <vtkGDCMPolyDataWriter.h>
  54
  55 //--------------------------------------------------------------------
  56 template<class PixelType>
  57 clitk::Xml2DicomRTStructFilter<PixelType>::Xml2DicomRTStructFilter()
  58 {
  59   m_ImageMHD = "";
  60   m_StructureSetFilename = "";
  61   m_DicomFolder = "";
  62   m_OutputFilename = "default-output.dcm";
  63 }
  64 //--------------------------------------------------------------------
  65
  66
  67 //--------------------------------------------------------------------
  68 template<class PixelType>
  69 clitk::Xml2DicomRTStructFilter<PixelType>::~Xml2DicomRTStructFilter()
  70 {
  71 }
  72 //--------------------------------------------------------------------
  73
  74
  75 //--------------------------------------------------------------------
  76 template<class PixelType>
  77 void clitk::Xml2DicomRTStructFilter<PixelType>::Update()
  78 {
  79   // Check this is a RT-Struct
  80   gdcm::Reader gdcm_reader;
  81   gdcm_reader.SetFileName(m_StructureSetFilename.c_str());
  82   if (!gdcm_reader.Read()) {
  83     clitkExceptionMacro("Error could not open the file '" << m_StructureSetFilename << std::endl);
  84   }
  85   gdcm::MediaStorage ms;
  86   ms.SetFromFile(gdcm_reader.GetFile());
  87   if (ms != gdcm::MediaStorage::RTStructureSetStorage) {
  88     clitkExceptionMacro("File '" << m_StructureSetFilename << "' is not a DICOM-RT-Struct file." << std::endl);
  89   }
  90
  91   // Read rt struct
  92   vtkSmartPointer<vtkGDCMPolyDataReader> reader = vtkGDCMPolyDataReader::New();
  93   reader->SetFileName(m_StructureSetFilename.c_str());
  94   reader->Update();
  95
  96   //Check id use pixel is set and if imageMHD is present in such case
  97   typedef itk::Image<double, 3> Input3DType;
  98   typedef itk::ImageFileReader<Input3DType> InputReader3DType;
  99   InputReader3DType::Pointer mhdReader = InputReader3DType::New();
 100   if (m_UsePixel) {
 101     if (m_ImageMHD != "") {
 102       mhdReader->SetFileName(m_ImageMHD);
 103       mhdReader->Update();
 104     } else {
 105       std::cout << "Set MHD image (option -j) when pixel is set" << std::endl;
 106       return;
 107     }
 108   }
 109
 110   // Get properties
 111   vtkRTStructSetProperties * p = reader->GetRTStructSetProperties();
 112   if (GetVerboseFlag()) {
 113     std::cout << "Number of structures in the dicom-rt-struct : "
 114               << p->GetNumberOfStructureSetROIs() << std::endl;
 115   }
 116
 117   std::map<std::string, vtkSmartPointer<vtkAppendPolyData> > mapName2Data;
 118
 119   pugi::xml_document doc;
 120   pugi::xml_parse_result result = doc.load_file(m_InputFilename.c_str());
 121   if (!result) {
 122     std::cout << "ERROR: no result" << std::endl;
 123     return;
 124   }
 125
 126   //Take the main dict and look for point_mm in each struct in each slice
 127   pugi::xml_node mainDict = doc.child("plist").child("dict").child("array");
 128   for (pugi::xml_node_iterator mainDictIt = mainDict.begin(); mainDictIt != mainDict.end(); ++mainDictIt) { //Look for all slice (one slice per dict)
 129     if (!std::strcmp(mainDictIt->name(), "dict")) {
 130       pugi::xml_node_iterator sliceDescriptionIt = mainDictIt->begin();
 131       while (sliceDescriptionIt != mainDictIt->end() && std::abs(std::strcmp(sliceDescriptionIt->child_value(), "ImageIndex"))) //Look for the ImageIndex for the current slice
 132         ++sliceDescriptionIt;
 133       int sliceNb = 0;
 134       if (sliceDescriptionIt != mainDictIt->end()) { //take the following node to have the values
 135         ++sliceDescriptionIt;
 136         sliceNb = std::atoi(sliceDescriptionIt->child_value());
 137       }
 138       for (pugi::xml_node_iterator sliceIt = mainDictIt->child("array").begin(); sliceIt != mainDictIt->child("array").end(); ++sliceIt) { //Look for all struct in the current slice (one struct per dict)
 139         if (!std::strcmp(sliceIt->name(), "dict")) {
 140           pugi::xml_node_iterator structureIt = sliceIt->begin();
 141           while (structureIt != sliceIt->end() && std::abs(std::strcmp(structureIt->child_value(), "Name"))) //Look for the name for the current struct
 142             ++structureIt;
 143           if (structureIt != sliceIt->end()) { //take the following node to have the values
 144             ++structureIt;
 145             std::string name(structureIt->child_value());
 146             if (m_UsePixel) {
 147               while (structureIt != sliceIt->end() && std::abs(std::strcmp(structureIt->child_value(), "Point_px"))) //Look for the Point_px for the current struct
 148                 ++structureIt;
 149             } else {
 150               while (structureIt != sliceIt->end() && std::abs(std::strcmp(structureIt->child_value(), "Point_mm"))) //Look for the Point_mm for the current struct
 151                 ++structureIt;
 152             }
 153             if (structureIt != sliceIt->end()) { //take the following node to have the values
 154               ++structureIt;
 155               // Insert all points for 1 struct into vtkPoints
 156               vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
 157               for (pugi::xml_node_iterator pointIt = structureIt->begin(); pointIt != structureIt->end(); ++pointIt) { //Look for all points in mm in the current struct (one point per string)
 158                 if (!std::strcmp(pointIt->name(), "string")) {
 159                   // Split the string to save the 2 or 3 values (remove the useless char) as double
 160                   char* copy = (char*)pointIt->child_value();
 161                   std::vector<char*> v;
 162                   char* chars_array = strtok(copy, ", ");
 163                   while(chars_array) {
 164                       v.push_back(chars_array);
 165                       chars_array = strtok(NULL, ", ");
 166                   }
 167                   double x, y, z;
 168                   v[0] = v[0] + 1;
 169                   if (m_UsePixel) {
 170                     v[1][strlen(v[1])-1] = '\0';
 171                     x = std::atof(v[0]);
 172                     y = std::atof(v[1]);
 173                     z = sliceNb;
 174                     Input3DType::PointType position;
 175                     Input3DType::IndexType index;
 176                     index[0] = x;
 177                     index[1] = y;
 178                     index[2] = z;
 179                     mhdReader->GetOutput()->TransformIndexToPhysicalPoint(index, position);
 180                     x = position[0];
 181                     y = position[1];
 182                     z = position[2];
 183                   } else {
 184                     v[2][strlen(v[2])-1] = '\0';
 185                     x = std::atof(v[0]);
 186                     y = std::atof(v[1]);
 187                     z = std::atof(v[2]);
 188                   }
 189                   points->InsertNextPoint(x, y, z);
 190                 }
 191               }
 192
 193               //Create the polygon attached to the points
 194               vtkSmartPointer<vtkPolygon> polygon = vtkSmartPointer<vtkPolygon>::New();
 195               polygon->GetPointIds()->SetNumberOfIds(points->GetNumberOfPoints()); //make a quad
 196               for (int pointIdx=0; pointIdx<points->GetNumberOfPoints(); ++pointIdx)
 197                 polygon->GetPointIds()->SetId(pointIdx, pointIdx);
 198               // Add the polygon to a list of polygons
 199               vtkSmartPointer<vtkCellArray> polygons = vtkSmartPointer<vtkCellArray>::New();
 200               polygons->InsertNextCell(polygon);
 201               // Create a PolyData
 202               vtkSmartPointer<vtkPolyData> polygonPolyData = vtkSmartPointer<vtkPolyData>::New();
 203               polygonPolyData->SetPoints(points);
 204               polygonPolyData->SetPolys(polygons);
 205
 206               //Append the polyData into the map at the correct stuct name entry
 207               std::map<std::string, vtkSmartPointer<vtkAppendPolyData> >::const_iterator it = mapName2Data.find(name);
 208               if (it == mapName2Data.end())
 209                 mapName2Data[name] = vtkSmartPointer<vtkAppendPolyData>::New();
 210 #if VTK_MAJOR_VERSION <= 5
 211               mapName2Data[name]->AddInput(polygonPolyData);
 212 #else
 213               mapName2Data[name]->AddInputData(polygonPolyData);
 214 #endif
 215             }
 216           }
 217         }
 218       }
 219     }
 220   }
 221
 222   for (std::map<std::string, vtkSmartPointer<vtkAppendPolyData> >::iterator it = mapName2Data.begin(); it != mapName2Data.end(); ++it)
 223     it->second->Update();
 224
 225   // number of contours
 226   int numMasks = mapName2Data.size();
 227
 228   // Init writer
 229   vtkGDCMPolyDataWriter * writer = vtkGDCMPolyDataWriter::New();
 230   writer->SetNumberOfInputPorts(numMasks);
 231   writer->SetFileName(m_OutputFilename.c_str());
 232   writer->SetMedicalImageProperties(reader->GetMedicalImageProperties());
 233
 234   // List of already present rois
 235   vtkStringArray* roiNames = vtkStringArray::New();
 236   vtkStringArray* roiAlgorithms = vtkStringArray::New();
 237   vtkStringArray* roiTypes = vtkStringArray::New();
 238   roiNames->SetNumberOfValues(numMasks);
 239   roiAlgorithms->SetNumberOfValues(numMasks);
 240   roiTypes->SetNumberOfValues(numMasks);
 241
 242   // Add new structures
 243   std::map<std::string, vtkSmartPointer<vtkAppendPolyData> >::iterator it = mapName2Data.begin();
 244   for (unsigned int i = 0; i < numMasks; ++i) {
 245 #if VTK_MAJOR_VERSION <= 5
 246     writer->SetInput(i, it->second->GetOutput());
 247 #else
 248     writer->SetInputData(i, it->second->GetOutput());
 249 #endif
 250     roiNames->InsertValue(i, it->first);
 251     roiAlgorithms->InsertValue(i, "CLITK_CREATED");
 252     roiTypes->InsertValue(i, "coucou"); //Roi type instead of coucou
 253     ++it;
 254   }
 255
 256   /*
 257   //  Visu DEBUG
 258   vtkPolyDataMapper *cubeMapper = vtkPolyDataMapper::New();
 259   cubeMapper->SetInputData( mapName2Data.begin()->second->GetOutput() );
 260   double range[2];
 261   mapName2Data.begin()->second->GetOutput()->GetPointData()->GetScalars()->GetRange(range);
 262   std::cout << numMasks << std::endl;
 263   cubeMapper->SetScalarRange(range);
 264   std::cout << range[0] << " " << range[1] << std::endl;
 265   vtkActor *cubeActor = vtkActor::New();
 266   cubeActor->SetMapper(cubeMapper);
 267   vtkProperty * property = cubeActor->GetProperty();
 268   property->SetRepresentationToWireframe();
 269
 270   vtkRenderer *renderer = vtkRenderer::New();
 271   vtkRenderWindow *renWin = vtkRenderWindow::New();
 272   renWin->AddRenderer(renderer);
 273
 274   vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
 275   iren->SetRenderWindow(renWin);
 276
 277   renderer->AddActor(cubeActor);
 278   renderer->ResetCamera();
 279   renderer->SetBackground(1,1,1);
 280
 281   renWin->SetSize(300,300);
 282
 283   renWin->Render();
 284   iren->Start();
 285   */
 286   // End visu
 287
 288   // Write (need to read dicom for slice id)
 289   vtkRTStructSetProperties* theProperties = vtkRTStructSetProperties::New();
 290   writer->SetRTStructSetProperties(theProperties);
 291   if (GetVerboseFlag()) {
 292     std::cout << "Looking for dicom info, study instance "
 293               << p->GetStudyInstanceUID() << std::endl;
 294   }
 295   writer->InitializeRTStructSet(m_DicomFolder,
 296                                 reader->GetRTStructSetProperties()->GetStructureSetLabel(),
 297                                 reader->GetRTStructSetProperties()->GetStructureSetName(),
 298                                 roiNames, roiAlgorithms, roiTypes);
 299   writer->Write();
 300   reader->Delete();
 301   roiNames->Delete();
 302   roiTypes->Delete();
 303   roiAlgorithms->Delete();
 304   writer->Delete();
 305 }
 306 //--------------------------------------------------------------------
 307