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