vv/vvMeshReader.cxx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================**/
  18 #include <algorithm>
  19
  20 #include <QApplication>
  21
  22 #include <gdcmFile.h>
  23 #if GDCM_MAJOR_VERSION == 2
  24   #include <gdcmReader.h>
  25   #include <gdcmTag.h>
  26   #include <gdcmAttribute.h>
  27 #else
  28   #include <gdcm.h>
  29   #include <gdcmSQItem.h>
  30 #endif
  31
  32 #include <vtkSmartPointer.h>
  33 #include <vtkAppendPolyData.h>
  34 #include <vtkCellArray.h>
  35 #include <vtkMetaImageWriter.h>
  36 #include <vtkMetaImageReader.h>
  37 #include <vtkPolyDataWriter.h>
  38 #include <vtkImageData.h>
  39 #include <vtkImageCast.h>
  40 #include <vtkImageGaussianSmooth.h>
  41
  42 #include "clitkCommon.h"
  43 #include "vvMeshReader.h"
  44 #include "vvProgressDialog.h"
  45 #include <clitkCommon.h>
  46
  47 vvMeshReader::vvMeshReader() :
  48   vtk_mode(false)
  49 {}
  50
  51 void vvMeshReader::Update()
  52 {
  53   //Show a progress bar only when opening a DC-struct (ie. multiple contours)
  54   vvProgressDialog progress("Opening " + filename,(!vtk_mode) && (selected_contours.size()>1));
  55   this->start();
  56   while (this->isRunning()) {
  57     progress.SetProgress(result.size(),selected_contours.size());
  58     this->wait(50);
  59     qApp->processEvents();
  60   }
  61 }
  62
  63 void vvMeshReader::run()
  64 {
  65   ///Verify the important stuff has been set
  66   assert(filename != "");
  67   if (vtk_mode) { //Read vtkPolyData
  68     vvMesh::Pointer m=vvMesh::New();
  69     m->ReadFromVTK(filename.c_str());
  70     if (vf) m->propagateContour(vf);
  71     m->ComputeMasks(image->GetVTKImages()[0],false); //don't extrude the contour
  72     result.push_back(m);
  73   } else { //Read a Dicom-struct file
  74     assert(selected_contours.size() > 0);
  75     assert(image);
  76     std::vector<vvMesh::Pointer> contour_stacks=readSelectedContours();
  77     for (std::vector<vvMesh::Pointer>::iterator i=contour_stacks.begin();
  78          i!=contour_stacks.end(); i++) {
  79       (*i)->ComputeMasks(image->GetVTKImages()[0],true); //Remesh the contour
  80       (*i)->ComputeMeshes();
  81       if (vf) (*i)->propagateContour(vf);
  82       result.push_back(*i);
  83     }
  84   }
  85 }
  86
  87 std::vector<std::pair<int,std::string> > vvMeshReader::GetROINames()
  88 {
  89   assert(filename!="");
  90   std::vector<std::pair<int, std::string> > roi_names;
  91 #if GDCM_MAJOR_VERSION == 2
  92   // duplicate code from  clitk::DicomRT_StructureSet::Read
  93   gdcm::Reader reader;
  94   reader.SetFileName( filename.c_str() );
  95   reader.Read();
  96
  97   const gdcm::DataSet &ds = reader.GetFile().GetDataSet();
  98
  99   gdcm::Tag tssroisq(0x3006,0x0020);
 100   const gdcm::DataElement &ssroisq = ds.GetDataElement( tssroisq );
 101   gdcm::SmartPointer<gdcm::SequenceOfItems> roi_seq = ssroisq.GetValueAsSQ();
 102   assert(roi_seq); // TODO error message
 103   for(unsigned int ridx = 0; ridx < roi_seq->GetNumberOfItems(); ++ridx)
 104     {
 105     gdcm::Item & item = roi_seq->GetItem( ridx + 1); // Item starts at 1
 106     const gdcm::DataSet& nestedds = item.GetNestedDataSet();
 107     if( nestedds.FindDataElement( gdcm::Tag(0x3006,0x22) ) )
 108       {
 109       gdcm::Attribute<0x3006,0x26> roiname;
 110       roiname.SetFromDataSet( nestedds );
 111       std::string name = roiname.GetValue();      // 0x3006,0x0026 = [ROI Name]
 112       gdcm::Attribute<0x3006,0x0022> roinumber;
 113       roinumber.SetFromDataSet( nestedds );
 114       int nb = roinumber.GetValue();  // 0x3006,0x0022 = [ROI Number]
 115
 116       roi_names.push_back(make_pair(nb,name));
 117       }
 118     }
 119 #else
 120   gdcm::File reader;
 121   reader.SetFileName(filename.c_str());
 122   reader.SetMaxSizeLoadEntry(16384);
 123   reader.Load();
 124
 125   gdcm::SeqEntry * roi_info=reader.GetSeqEntry(0x3006,0x0020);
 126   assert(roi_info);
 127   // DD("ici");
 128   //int n=0;
 129   for (gdcm::SQItem* i=roi_info->GetFirstSQItem(); i!=0; i=roi_info->GetNextSQItem())
 130     if (i->GetEntryValue(0x3006,0x0022)!= gdcm::GDCM_UNFOUND)
 131       roi_names.push_back(make_pair(atoi(i->GetEntryValue(0x3006,0x0022).c_str()),i->GetEntryValue(0x3006,0x0026)));
 132 #endif
 133   return roi_names;
 134 }
 135
 136 std::vector<vvMesh::Pointer> vvMeshReader::readSelectedContours()
 137 {
 138   std::vector<vvMesh::Pointer> result;
 139 #if GDCM_MAJOR_VERSION == 2
 140 #else
 141   gdcm::File reader;
 142   reader.SetFileName(filename.c_str());
 143   reader.SetMaxSizeLoadEntry(16384);
 144   reader.Load();
 145
 146   gdcm::SeqEntry * rois=reader.GetSeqEntry(0x3006,0x0039);
 147   ///We need to iterate both on the contours themselves, and on the contour info
 148   gdcm::SeqEntry * roi_info=reader.GetSeqEntry(0x3006,0x0020);
 149   gdcm::SQItem* k=roi_info->GetFirstSQItem();
 150   for(gdcm::SQItem* i=rois->GetFirstSQItem(); i!=0; i=rois->GetNextSQItem()) { //loop over ROIS
 151     assert(k!=0);
 152     vtkSmartPointer<vtkAppendPolyData> append=vtkSmartPointer<vtkAppendPolyData>::New();
 153     std::istringstream ss(i->GetEntryValue(0x3006,0x0084));
 154     int roi_number;
 155     ss >> roi_number;
 156     if (std::find(selected_contours.begin(),selected_contours.end(),roi_number) != selected_contours.end()) { //Only read selected ROIs
 157       vvMesh::Pointer current_roi=vvMesh::New();
 158       std::vector<double> rgb=clitk::parse_string<double>(i->GetEntryValue(0x3006,0x002a),'\\');
 159       assert(rgb.size()==3);
 160       current_roi->r=rgb[0]/255;
 161       current_roi->g=rgb[1]/255;
 162       current_roi->b=rgb[2]/255;
 163       current_roi->structure_name=k->GetEntryValue(0x3006,0x0026);
 164       gdcm::SeqEntry * contours=i->GetSeqEntry(0x3006,0x0040);
 165       double z0=-1; //Used to determine spacing between slices, assumed to be constant
 166       for(gdcm::SQItem* j=contours->GetFirstSQItem(); j!=0; j=contours->GetNextSQItem()) { //loop over 2D contours
 167         std::string contour_type=j->GetEntryValue(0x3006,0x0042);
 168         if (contour_type=="CLOSED_PLANAR ") {
 169           int point_number=clitk::parse_value<int>(j->GetEntryValue(0x3006,0x0046));
 170           std::vector<float> points=clitk::parse_string<float>(j->GetEntryValue(0x3006,0x0050),'\\');
 171           assert(points.size() == static_cast<unsigned int>(point_number)*3);
 172           if (z0 == -1) //First contour
 173             z0=points[2];
 174           //2nd contour, spacing not yet set. Need to be sure we are on a different slice,
 175           //sometimes there is more than one closed contour per slice
 176           else if (current_roi->GetSpacing()==-1 && points[2] != z0 )
 177             current_roi->SetSpacing(points[2]-z0);
 178           vtkPolyData * contour=vtkPolyData::New();
 179           contour->Allocate(); //for cell structures
 180           contour->SetPoints(vtkPoints::New());
 181           vtkIdType ids[2];
 182           for (unsigned int idx=0; idx<points.size(); idx+=3) {
 183             contour->GetPoints()->InsertNextPoint(points[idx],points[idx+1],points[idx+2]);
 184             ids[0]=idx/3;
 185             ids[1]=(ids[0]+1)%point_number; //0-1,1-2,...,n-1-0
 186             contour->GetLines()->InsertNextCell(2,ids);
 187           }
 188           append->AddInput(contour);
 189         } else if (contour_type == "POINT ")
 190           ; // silently ignore POINT type since we don't need them at the moment
 191         else
 192           std::cerr << "Warning: contour type " << contour_type << " not handled!" << std::endl;
 193       }
 194       append->Update();
 195       current_roi->AddMesh(append->GetOutput());
 196       result.push_back(current_roi);
 197     } else {
 198       //std::cerr << "Warning: ignoring ROI #" << roi_number << std::endl;
 199     }
 200     k=roi_info->GetNextSQItem(); //increment the second loop variable
 201   }
 202 #endif
 203   return result;
 204 }
 205