X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=vv%2FvvMeshReader.cxx;h=b04163853d202f14f1a435f89b50c9444f38a036;hb=HEAD;hp=980e17d157346cc42e6f6e3e151a94b80a7a5277;hpb=0b7c9b1e1215634b02cbd38d4e4ba101d6111ba8;p=clitk.git

diff --git a/vv/vvMeshReader.cxx b/vv/vvMeshReader.cxx
index 980e17d..b041638 100644
--- a/vv/vvMeshReader.cxx
+++ b/vv/vvMeshReader.cxx
@@ -1,9 +1,9 @@
 /*=========================================================================
   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
 
-  Authors belong to: 
+  Authors belong to:
   - University of LYON              http://www.universite-lyon.fr/
-  - Léon Bérard cancer center       http://oncora1.lyon.fnclcc.fr
+  - Léon Bérard cancer center       http://www.centreleonberard.fr
   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
 
   This software is distributed WITHOUT ANY WARRANTY; without even
@@ -14,14 +14,27 @@
 
   - BSD        See included LICENSE.txt file
   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
-======================================================================-====*/
+  ===========================================================================**/
+
+// std
 #include <algorithm>
 
+// qt
 #include <QApplication>
 
+// gdcm 
+#include <gdcmFile.h>
+#if GDCM_MAJOR_VERSION >= 2
+#include <gdcmReader.h>
+#include <gdcmTag.h>
+#include <gdcmAttribute.h>
+#else
 #include <gdcm.h>
 #include <gdcmSQItem.h>
+#endif
 
+// vtk
+#include <vtkVersion.h>
 #include <vtkSmartPointer.h>
 #include <vtkAppendPolyData.h>
 #include <vtkCellArray.h>
@@ -32,165 +45,336 @@
 #include <vtkImageCast.h>
 #include <vtkImageGaussianSmooth.h>
 
+// clitk
+#include "clitkCommon.h"
 #include "clitkCommon.h"
+#include "clitkDicomRT_StructureSet.h"
 #include "vvMeshReader.h"
 #include "vvProgressDialog.h"
 
+//------------------------------------------------------------------------------
 vvMeshReader::vvMeshReader() :
-    vtk_mode(false)
+vtk_mode(false)
 {}
+//------------------------------------------------------------------------------
+
 
+//------------------------------------------------------------------------------
 void vvMeshReader::Update()
 {
-    //Show a progress bar only when opening a DC-struct (ie. multiple contours)
-    vvProgressDialog progress("Opening " + filename,(!vtk_mode) && (selected_contours.size()>1));
-    this->start();
-    while (this->isRunning())
-    {
-        progress.SetProgress(result.size(),selected_contours.size());
-        this->wait(50);
-        qApp->processEvents();
-    }
+  //Show a progress bar only when opening a DC-struct (ie. multiple contours)
+  vvProgressDialog progress("Opening " + filename,(!vtk_mode) && (selected_contours.size()>1));
+  this->start();
+  while (this->isRunning()) {
+    progress.SetProgress(result.size(),selected_contours.size());
+    this->wait(50);
+    qApp->processEvents();
+  }
 }
+//------------------------------------------------------------------------------
+
 
+//------------------------------------------------------------------------------
 void vvMeshReader::run()
 {
-    ///Verify the important stuff has been set
-    assert(filename != "");
-    if (vtk_mode) //Read vtkPolyData
-    {
-        vvMesh::Pointer m=vvMesh::New();
-        m->ReadFromVTK(filename.c_str());
-        if (vf) m->propagateContour(vf);
-        m->ComputeMasks(image->GetVTKImages()[0],true);
-        result.push_back(m);
-    }
-    else //Read a Dicom-struct file
-    {
-        assert(selected_contours.size() > 0);
-        assert(image);
-        std::vector<vvMesh::Pointer> contour_stacks=readSelectedContours();
-        for (std::vector<vvMesh::Pointer>::iterator i=contour_stacks.begin();
-                i!=contour_stacks.end();i++)
-        {
-            (*i)->ComputeMasks(image->GetVTKImages()[0],true); //Remesh the contour
-            (*i)->ComputeMeshes();
-            if (vf) (*i)->propagateContour(vf);
-            result.push_back(*i);
-        }
+  ///Verify the important stuff has been set
+  assert(filename != "");
+  if (vtk_mode) { //Read vtkPolyData
+    vvMesh::Pointer m=vvMesh::New();
+    m->ReadFromVTK(filename.c_str());
+    if (vf) m->propagateContour(vf);
+    m->ComputeMasks(image->GetVTKImages()[0],false); //don't extrude the contour
+    result.push_back(m);
+  } else { //Read a Dicom-struct file
+    assert(selected_contours.size() > 0);
+    assert(image);
+    std::vector<vvMesh::Pointer> contour_stacks=readSelectedContours();
+    for (std::vector<vvMesh::Pointer>::iterator i=contour_stacks.begin();
+         i!=contour_stacks.end(); i++) {
+      (*i)->ComputeMasks(image->GetVTKImages()[0],true); //Remesh the contour
+      (*i)->ComputeMeshes();
+      if (vf) (*i)->propagateContour(vf);
+      result.push_back(*i);
     }
+  }
 }
+//------------------------------------------------------------------------------
 
-template<class ElementType>
-ElementType parse_value(std::string str)
-{
-    std::istringstream parser(str);
-    ElementType value;
-    parser >> value;
-    assert(!parser.fail());
-    return value;
-}
 
-template<class ElementType>
-std::vector<ElementType> parse_string(std::string str,char delim)
+//------------------------------------------------------------------------------
+std::vector<std::pair<int,std::string> > vvMeshReader::GetROINames()
 {
-    std::istringstream ss(str);
-    std::string token;
-    std::vector<ElementType> result;
-    while (getline(ss,token,delim))
-    {
-        result.push_back(parse_value<ElementType>(token));
+  assert(filename!="");
+  std::vector<std::pair<int, std::string> > roi_names;
+
+#if CLITK_USE_SYSTEM_GDCM == 1
+
+  // Read RT-struct data
+  vtkSmartPointer<vtkGDCMPolyDataReader> areader = vtkGDCMPolyDataReader::New();
+  areader->SetFileName(filename.c_str());
+  areader->Update();
+
+  // get info on roi names
+  vtkRTStructSetProperties * p = areader->GetRTStructSetProperties();
+  int n = p->GetNumberOfStructureSetROIs();
+  
+  for(unsigned int i=0; i<n; i++) {
+    std::string name = p->GetStructureSetROIName(i);
+    int nb = p->GetStructureSetROINumber(i);
+    roi_names.push_back(make_pair(nb,name));
+  }
+
+#else
+#if GDCM_MAJOR_VERSION >= 2
+
+    // duplicate code from  clitk::DicomRT_StructureSet::Read
+    gdcm::Reader * reader = new gdcm::Reader;
+    reader->SetFileName( filename.c_str() );
+    reader->Read();
+
+    const gdcm::DataSet &ds = reader->GetFile().GetDataSet();
+
+    // Check file type
+    //Verify if the file is a RT-Structure-Set dicom file
+    gdcm::File * mFile = &(reader->GetFile());
+    gdcm::MediaStorage ms;
+    ms.SetFromFile(*mFile);
+    if( ms != gdcm::MediaStorage::RTStructureSetStorage )
+      {
+        std::cerr << "Error. the file " << filename
+                  << " is not a Dicom Struct ? (must have a SOP Class UID [0008|0016] = 1.2.840.10008.5.1.4.1.1.481.3 ==> [RT Structure Set Storage])"
+                  << std::endl;
+        exit(0);
+      }
+
+    gdcm::Attribute<0x8,0x60> modality;
+    modality.SetFromDataSet( ds );
+    if( modality.GetValue() != "RTSTRUCT" )
+      {
+        std::cerr << "Error. the file " << filename
+                  << " is not a Dicom Struct ? (must have 0x0008,0x0060 = RTSTRUCT [RT Structure Set Storage])"
+                  << std::endl;
+        exit(0);
+      }
+
+    gdcm::Attribute<0x20,0x10> studyid;
+    studyid.SetFromDataSet( ds );
+
+    gdcm::Tag tssroisq(0x3006,0x0020);
+    // 0x3006,0x0020 = [ Structure Set ROI Sequence ]
+    bool b = ds.FindDataElement(tssroisq);
+    if (!b) { // FIXME
+      clitkExceptionMacro("Error: tag 0x3006,0x0020 [ Structure Set ROI Sequence ] not found");
     }
-    return result;
-}
+  
+    const gdcm::DataElement &ssroisq = ds.GetDataElement( tssroisq );
+    gdcm::SmartPointer<gdcm::SequenceOfItems> roi_seq = ssroisq.GetValueAsSQ();
+    assert(roi_seq); // FIXME error message
+  
+    for(unsigned int ridx = 0; ridx < roi_seq->GetNumberOfItems(); ++ridx)
+      {
+         gdcm::Item & item = roi_seq->GetItem( ridx + 1); // Item starts at 1
 
-std::vector<std::pair<int,std::string> > vvMeshReader::GetROINames()
-{
-    assert(filename!="");
-    gdcm::File reader;
-    reader.SetFileName(filename.c_str());
-    reader.SetMaxSizeLoadEntry(16384);
-    reader.Load();
-
-    gdcm::SeqEntry * roi_info=reader.GetSeqEntry(0x3006,0x0020);
-    assert(roi_info);
-    std::vector<std::pair<int, std::string> > roi_names;
-    // DD("ici");
-    //int n=0;
-    for (gdcm::SQItem* i=roi_info->GetFirstSQItem();i!=0;i=roi_info->GetNextSQItem())
-        if (i->GetEntryValue(0x3006,0x0022)!= gdcm::GDCM_UNFOUND)
-            roi_names.push_back(make_pair(atoi(i->GetEntryValue(0x3006,0x0022).c_str()),i->GetEntryValue(0x3006,0x0026)));
-    return roi_names;
+        const gdcm::Item & sitem = roi_seq->GetItem(ridx+1); // Item start at #1   
+
+        const gdcm::DataSet& snestedds = sitem.GetNestedDataSet();
+        const gdcm::DataSet& nestedds = item.GetNestedDataSet();
+
+        if( snestedds.FindDataElement( gdcm::Tag(0x3006,0x22) ) )
+          {
+            // const gdcm::DataElement & a = nestedds.GetDataElement(gdcm::Tag(0x3006,0x26));
+            // DD(a.GetValue());
+
+            gdcm::Attribute<0x3006,0x26> roiname;
+            roiname.SetFromDataSet( snestedds );
+            std::string name = roiname.GetValue();      // 0x3006,0x0026 = [ROI Name]
+            gdcm::Attribute<0x3006,0x0022> roinumber;
+            roinumber.SetFromDataSet( snestedds );
+            int nb = roinumber.GetValue();  // 0x3006,0x0022 = [ROI Number]
+          
+            roi_names.push_back(make_pair(nb,name));
+          }
+      }
+  
+    delete reader;
+
+#else
+  gdcm::File reader;
+  reader.SetFileName(filename.c_str());
+  reader.SetMaxSizeLoadEntry(16384);
+  reader.Load();
+
+  gdcm::SeqEntry * roi_info=reader.GetSeqEntry(0x3006,0x0020);
+  assert(roi_info);
+  // DD("ici");
+  //int n=0;
+  for (gdcm::SQItem* i=roi_info->GetFirstSQItem(); i!=0; i=roi_info->GetNextSQItem())
+    if (i->GetEntryValue(0x3006,0x0022)!= gdcm::GDCM_UNFOUND)
+      roi_names.push_back(make_pair(atoi(i->GetEntryValue(0x3006,0x0022).c_str()),i->GetEntryValue(0x3006,0x0026)));
+#endif
+#endif
+
+  return roi_names;
 }
+//------------------------------------------------------------------------------
 
+
+//------------------------------------------------------------------------------
 std::vector<vvMesh::Pointer> vvMeshReader::readSelectedContours()
 {
-    gdcm::File reader;
-    reader.SetFileName(filename.c_str());
-    reader.SetMaxSizeLoadEntry(16384);
-    reader.Load();
-
-    std::vector<vvMesh::Pointer> result;
-    gdcm::SeqEntry * rois=reader.GetSeqEntry(0x3006,0x0039);
-    ///We need to iterate both on the contours themselves, and on the contour info
-    gdcm::SeqEntry * roi_info=reader.GetSeqEntry(0x3006,0x0020);
-    gdcm::SQItem* k=roi_info->GetFirstSQItem();
-    for(gdcm::SQItem* i=rois->GetFirstSQItem();i!=0;i=rois->GetNextSQItem()) //loop over ROIS
+  std::vector<vvMesh::Pointer> result;
+#if GDCM_MAJOR_VERSION >= 2
+  gdcm::Reader reader;
+  reader.SetFileName(filename.c_str());
+  reader.Read();
+
+  const gdcm::DataSet &ds = reader.GetFile().GetDataSet();
+
+  gdcm::SmartPointer<gdcm::SequenceOfItems> rois = ds.GetDataElement(gdcm::Tag(0x3006,0x39)).GetValueAsSQ();
+  gdcm::SmartPointer<gdcm::SequenceOfItems> roi_info = ds.GetDataElement(gdcm::Tag(0x3006,0x20)).GetValueAsSQ();
+  assert(rois); // TODO error message
+  assert(roi_info); // TODO error message
+  assert(rois->GetNumberOfItems() == roi_info->GetNumberOfItems());
+
+  for (unsigned ridx = 0; ridx < rois->GetNumberOfItems(); ++ridx)
     {
-        assert(k!=0);
-        vtkSmartPointer<vtkAppendPolyData> append=vtkSmartPointer<vtkAppendPolyData>::New();
-        std::istringstream ss(i->GetEntryValue(0x3006,0x0084));
-        int roi_number;ss >> roi_number;
-        if (std::find(selected_contours.begin(),selected_contours.end(),roi_number) != selected_contours.end())//Only read selected ROIs
+      vtkSmartPointer<vtkAppendPolyData> append=vtkSmartPointer<vtkAppendPolyData>::New();
+      const gdcm::DataSet& ds_rois = rois->GetItem( ridx + 1).GetNestedDataSet();
+      const gdcm::DataSet& ds_roi_info = roi_info->GetItem( ridx + 1).GetNestedDataSet();
+
+      gdcm::Attribute<0x3006,0x84> roinumber;
+      roinumber.SetFromDataSet(ds_rois);
+      if (std::find(selected_contours.begin(), selected_contours.end(), roinumber.GetValue()) != selected_contours.end()) //Only read selected ROIs
         {
-            vvMesh::Pointer current_roi=vvMesh::New();
-            std::vector<double> rgb=parse_string<double>(i->GetEntryValue(0x3006,0x002a),'\\');
-            assert(rgb.size()==3);
-            current_roi->r=rgb[0]/255; current_roi->g=rgb[1]/255; current_roi->b=rgb[2]/255;
-            current_roi->structure_name=k->GetEntryValue(0x3006,0x0026);
-            gdcm::SeqEntry * contours=i->GetSeqEntry(0x3006,0x0040);
-            double z0=-1; //Used to determine spacing between slices, assumed to be constant
-            for(gdcm::SQItem* j=contours->GetFirstSQItem();j!=0;j=contours->GetNextSQItem()) //loop over 2D contours
+          gdcm::Attribute<0x3006,0x2a> trgb;
+          trgb.SetFromDataSet(ds_rois);
+          vvMesh::Pointer current_roi=vvMesh::New();
+          current_roi->r = trgb[0] / 255.0;
+          current_roi->g = trgb[1] / 255.0;
+          current_roi->b = trgb[2] / 255.0;
+
+          gdcm::Attribute<0x3006,0x26> tstructure_name;
+          tstructure_name.SetFromDataSet(ds_roi_info);
+          current_roi->structure_name = tstructure_name.GetValue();
+
+          gdcm::SmartPointer<gdcm::SequenceOfItems> roi_seq = ds_rois.GetDataElement(gdcm::Tag(0x3006,0x40)).GetValueAsSQ();
+          double z0=-1; //Used to determine spacing between slices, assumed to be constant
+          for (unsigned j = 0; j < roi_seq->GetNumberOfItems(); ++j)
             {
-                std::string contour_type=j->GetEntryValue(0x3006,0x0042);
-                if (contour_type=="CLOSED_PLANAR ")
+              gdcm::Item & item_roi_seq = roi_seq->GetItem(j + 1); // Item starts at 1
+              const gdcm::DataSet& ds_roi_seq = item_roi_seq.GetNestedDataSet();
+              gdcm::Attribute<0x3006,0x42> tcontour_type;
+              tcontour_type.SetFromDataSet(ds_roi_seq);
+              std::string contour_type = tcontour_type.GetValue();
+              if (contour_type=="CLOSED_PLANAR ")
                 {
-                    int point_number=parse_value<int>(j->GetEntryValue(0x3006,0x0046));
-                    std::vector<float> points=parse_string<float>(j->GetEntryValue(0x3006,0x0050),'\\');
-                    assert(points.size() == static_cast<unsigned int>(point_number)*3);
-                    if (z0 == -1) //First contour
-                        z0=points[2];
-                    //2nd contour, spacing not yet set. Need to be sure we are on a different slice,
-                    //sometimes there is more than one closed contour per slice
-                    else if (current_roi->GetSpacing()==-1 && points[2] != z0 ) 
-                        current_roi->SetSpacing(points[2]-z0);
-                    vtkPolyData * contour=vtkPolyData::New();
-                    contour->Allocate(); //for cell structures
-                    contour->SetPoints(vtkPoints::New());
-                    vtkIdType ids[2];
-                    for (unsigned int idx=0;idx<points.size();idx+=3)
+                  gdcm::Attribute<0x3006,0x46> tpoint_number;
+                  tpoint_number.SetFromDataSet(ds_roi_seq);
+                  const gdcm::DataElement & points_data = ds_roi_seq.GetDataElement(gdcm::Tag(0x3006,0x50));
+                  gdcm::Attribute<0x3006,0x50> tpoints;
+                  tpoints.SetFromDataElement(points_data);
+                  assert(tpoints.GetNumberOfValues() == static_cast<unsigned int>(tpoint_number.GetValue()) * 3);
+                  const double* points = tpoints.GetValues();
+                  if (z0 == -1) //First contour
+                    z0=points[2];
+                  else
+                    if (current_roi->GetSpacing()==-1 && points[2] != z0 )
+                      current_roi->SetSpacing(points[2]-z0);
+                  vtkPolyData * contour=vtkPolyData::New();
+                  contour->Allocate(); //for cell structures
+                  contour->SetPoints(vtkPoints::New());
+                  vtkIdType ids[2];
+                  for (unsigned idx = 0; idx < tpoints.GetNumberOfValues(); idx += 3)
                     {
-                        contour->GetPoints()->InsertNextPoint(points[idx],points[idx+1],points[idx+2]);
-                        ids[0]=idx/3;ids[1]=(ids[0]+1)%point_number; //0-1,1-2,...,n-1-0
-                        contour->GetLines()->InsertNextCell(2,ids);
+                      contour->GetPoints()->InsertNextPoint(points[idx], points[idx+1], points[idx+2]);
+                      ids[0] = idx / 3;
+                      ids[1] = (ids[0] + 1) % tpoint_number.GetValue(); //0-1,1-2,...,n-1-0
+                      contour->GetLines()->InsertNextCell(2, ids);
                     }
-                    append->AddInput(contour);
+#if VTK_MAJOR_VERSION <= 5
+                  append->AddInput(contour);
+#else
+                  append->AddInputData(contour);
+#endif
                 }
-                else if (contour_type == "POINT ")
-                    ; // silently ignore POINT type since we don't need them at the moment
+              else
+                if (contour_type == "POINT ")
+                  ; // silently ignore POINT type since we don't need them at the moment
                 else
-                    std::cerr << "Warning: contour type " << contour_type << " not handled!" << std::endl;
+                  std::cerr << "Warning: contour type " << contour_type << " not handled!" << std::endl;
             }
-            append->Update();
-            current_roi->AddMesh(append->GetOutput());
-            result.push_back(current_roi);
+          append->Update();
+          current_roi->AddMesh(append->GetOutput());
+          result.push_back(current_roi);
         }
-        else
+      else
         {
-            //std::cerr << "Warning: ignoring ROI #" << roi_number << std::endl;
+          //std::cerr << "Warning: ignoring ROI #" << roi_number << std::endl;
         }
-        k=roi_info->GetNextSQItem(); //increment the second loop variable
     }
-    return result;
+#else
+  gdcm::File reader;
+  reader.SetFileName(filename.c_str());
+  reader.SetMaxSizeLoadEntry(16384);
+  reader.Load();
+
+  gdcm::SeqEntry * rois=reader.GetSeqEntry(0x3006,0x0039);
+  ///We need to iterate both on the contours themselves, and on the contour info
+  gdcm::SeqEntry * roi_info=reader.GetSeqEntry(0x3006,0x0020);
+  gdcm::SQItem* k=roi_info->GetFirstSQItem();
+  for(gdcm::SQItem* i=rois->GetFirstSQItem(); i!=0; i=rois->GetNextSQItem()) { //loop over ROIS
+    assert(k!=0);
+    vtkSmartPointer<vtkAppendPolyData> append=vtkSmartPointer<vtkAppendPolyData>::New();
+    std::istringstream ss(i->GetEntryValue(0x3006,0x0084));
+    int roi_number;
+    ss >> roi_number;
+    if (std::find(selected_contours.begin(),selected_contours.end(),roi_number) != selected_contours.end()) { //Only read selected ROIs
+      vvMesh::Pointer current_roi=vvMesh::New();
+      std::vector<double> rgb=clitk::parse_string<double>(i->GetEntryValue(0x3006,0x002a),'\\');
+      assert(rgb.size()==3);
+      current_roi->r=rgb[0]/255;
+      current_roi->g=rgb[1]/255;
+      current_roi->b=rgb[2]/255;
+      current_roi->structure_name=k->GetEntryValue(0x3006,0x0026);
+      gdcm::SeqEntry * contours=i->GetSeqEntry(0x3006,0x0040);
+      double z0=-1; //Used to determine spacing between slices, assumed to be constant
+      for(gdcm::SQItem* j=contours->GetFirstSQItem(); j!=0; j=contours->GetNextSQItem()) { //loop over 2D contours
+        std::string contour_type=j->GetEntryValue(0x3006,0x0042);
+        if (contour_type=="CLOSED_PLANAR ") {
+          int point_number=clitk::parse_value<int>(j->GetEntryValue(0x3006,0x0046));
+          std::vector<float> points=clitk::parse_string<float>(j->GetEntryValue(0x3006,0x0050),'\\');
+          assert(points.size() == static_cast<unsigned int>(point_number)*3);
+          if (z0 == -1) //First contour
+            z0=points[2];
+          //2nd contour, spacing not yet set. Need to be sure we are on a different slice,
+          //sometimes there is more than one closed contour per slice
+          else if (current_roi->GetSpacing()==-1 && points[2] != z0 )
+            current_roi->SetSpacing(points[2]-z0);
+          vtkPolyData * contour=vtkPolyData::New();
+          contour->Allocate(); //for cell structures
+          contour->SetPoints(vtkPoints::New());
+          vtkIdType ids[2];
+          for (unsigned int idx=0; idx<points.size(); idx+=3) {
+            contour->GetPoints()->InsertNextPoint(points[idx],points[idx+1],points[idx+2]);
+            ids[0]=idx/3;
+            ids[1]=(ids[0]+1)%point_number; //0-1,1-2,...,n-1-0
+            contour->GetLines()->InsertNextCell(2,ids);
+          }
+          append->AddInput(contour);
+        } else if (contour_type == "POINT ")
+          ; // silently ignore POINT type since we don't need them at the moment
+        else
+          std::cerr << "Warning: contour type " << contour_type << " not handled!" << std::endl;
+      }
+      append->Update();
+      current_roi->AddMesh(append->GetOutput());
+      result.push_back(current_roi);
+    } else {
+      //std::cerr << "Warning: ignoring ROI #" << roi_number << std::endl;
+    }
+    k=roi_info->GetNextSQItem(); //increment the second loop variable
+  }
+#endif
+  return result;
 }