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
  19 // std
  20 #include <algorithm>
  21
  22 // qt
  23 #include <QApplication>
  24
  25 // gdcm
  26 #include <gdcmFile.h>
  27 #if GDCM_MAJOR_VERSION == 2
  28 #include <gdcmReader.h>
  29 #include <gdcmTag.h>
  30 #include <gdcmAttribute.h>
  31 #else
  32 #include <gdcm.h>
  33 #include <gdcmSQItem.h>
  34 #endif
  35
  36 // vtk
  37 #include <vtkSmartPointer.h>
  38 #include <vtkAppendPolyData.h>
  39 #include <vtkCellArray.h>
  40 #include <vtkMetaImageWriter.h>
  41 #include <vtkMetaImageReader.h>
  42 #include <vtkPolyDataWriter.h>
  43 #include <vtkImageData.h>
  44 #include <vtkImageCast.h>
  45 #include <vtkImageGaussianSmooth.h>
  46
  47 // clitk
  48 #include "clitkCommon.h"
  49 #include "clitkCommon.h"
  50 #include "clitkDicomRT_StructureSet.h"
  51 #include "vvMeshReader.h"
  52 #include "vvProgressDialog.h"
  53
  54 //------------------------------------------------------------------------------
  55 vvMeshReader::vvMeshReader() :
  56 vtk_mode(false)
  57 {}
  58 //------------------------------------------------------------------------------
  59
  60
  61 //------------------------------------------------------------------------------
  62 void vvMeshReader::Update()
  63 {
  64   DD("vvMeshReader::Update");
  65
  66   //Show a progress bar only when opening a DC-struct (ie. multiple contours)
  67   vvProgressDialog progress("Opening " + filename,(!vtk_mode) && (selected_contours.size()>1));
  68   this->start();
  69   while (this->isRunning()) {
  70     progress.SetProgress(result.size(),selected_contours.size());
  71     this->wait(50);
  72     qApp->processEvents();
  73   }
  74 }
  75 //------------------------------------------------------------------------------
  76
  77
  78 //------------------------------------------------------------------------------
  79 void vvMeshReader::run()
  80 {
  81   ///Verify the important stuff has been set
  82   assert(filename != "");
  83   if (vtk_mode) { //Read vtkPolyData
  84     vvMesh::Pointer m=vvMesh::New();
  85     m->ReadFromVTK(filename.c_str());
  86     if (vf) m->propagateContour(vf);
  87     m->ComputeMasks(image->GetVTKImages()[0],false); //don't extrude the contour
  88     result.push_back(m);
  89   } else { //Read a Dicom-struct file
  90     assert(selected_contours.size() > 0);
  91     assert(image);
  92     std::vector<vvMesh::Pointer> contour_stacks=readSelectedContours();
  93     for (std::vector<vvMesh::Pointer>::iterator i=contour_stacks.begin();
  94          i!=contour_stacks.end(); i++) {
  95       (*i)->ComputeMasks(image->GetVTKImages()[0],true); //Remesh the contour
  96       (*i)->ComputeMeshes();
  97       if (vf) (*i)->propagateContour(vf);
  98       result.push_back(*i);
  99     }
 100   }
 101 }
 102 //------------------------------------------------------------------------------
 103
 104
 105 //------------------------------------------------------------------------------
 106 std::vector<std::pair<int,std::string> > vvMeshReader::GetROINames()
 107 {
 108   assert(filename!="");
 109   std::vector<std::pair<int, std::string> > roi_names;
 110   DD(GDCM_MAJOR_VERSION);
 111   DD(CLITK_USE_SYSTEM_GDCM);
 112
 113 #if GDCM_MAJOR_VERSION == 2
 114
 115   // Read RT-struct data
 116   DD("before read");
 117   vtkSmartPointer<vtkGDCMPolyDataReader> reader = vtkGDCMPolyDataReader::New();
 118   reader->SetFileName(filename.c_str());
 119   reader->Update();
 120   DD("after read");
 121
 122   // get info on roi names
 123   vtkRTStructSetProperties * p = reader->GetRTStructSetProperties();
 124   DD(p->GetNumberOfStructureSetROIs());
 125   DD(p->GetStructureSetROIName(0));
 126   DD(p->GetStructureSetROINumber(0));
 127   int n = p->GetNumberOfStructureSetROIs();
 128   DD(n);
 129
 130   for(unsigned int i=0; i<n; i++) {
 131     std::string name = p->GetStructureSetROIName(i);
 132     int nb = p->GetStructureSetROINumber(i);
 133     roi_names.push_back(make_pair(nb,name));
 134   }
 135
 136   // ====================
 137
 138   if (false) {
 139     // duplicate code from  clitk::DicomRT_StructureSet::Read
 140     gdcm::Reader * reader = new gdcm::Reader;
 141     reader->SetFileName( filename.c_str() );
 142     reader->Read();
 143     DD("after gdcm read");
 144
 145     const gdcm::DataSet &ds = reader->GetFile().GetDataSet();
 146
 147     DD("after ds");
 148     DD(ds.IsEmpty());
 149
 150     // Check file type
 151     //Verify if the file is a RT-Structure-Set dicom file
 152     gdcm::File * mFile = &(reader->GetFile());
 153     gdcm::MediaStorage ms;
 154     ms.SetFromFile(*mFile);
 155     if( ms != gdcm::MediaStorage::RTStructureSetStorage )
 156       {
 157         std::cerr << "Error. the file " << filename
 158                   << " 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])"
 159                   << std::endl;
 160         exit(0);
 161       }
 162
 163     gdcm::Attribute<0x8,0x60> modality;
 164     modality.SetFromDataSet( ds );
 165     if( modality.GetValue() != "RTSTRUCT" )
 166       {
 167         std::cerr << "Error. the file " << filename
 168                   << " is not a Dicom Struct ? (must have 0x0008,0x0060 = RTSTRUCT [RT Structure Set Storage])"
 169                   << std::endl;
 170         exit(0);
 171       }
 172
 173     gdcm::Attribute<0x20,0x10> studyid;
 174     studyid.SetFromDataSet( ds );
 175     DD(studyid.GetValue());
 176
 177     gdcm::Tag tssroisq(0x3006,0x0020);
 178     // 0x3006,0x0020 = [ Structure Set ROI Sequence ]
 179     bool b = ds.FindDataElement(tssroisq);
 180     if (!b) { // FIXME
 181       clitkExceptionMacro("Error: tag 0x3006,0x0020 [ Structure Set ROI Sequence ] not found");
 182     }
 183
 184     const gdcm::DataElement &ssroisq = ds.GetDataElement( tssroisq );
 185     DD("after ssroisq");
 186     gdcm::SmartPointer<gdcm::SequenceOfItems> roi_seq = ssroisq.GetValueAsSQ();
 187     assert(roi_seq); // FIXME error message
 188
 189     DD(roi_seq->GetNumberOfItems());
 190
 191     for(unsigned int ridx = 0; ridx < roi_seq->GetNumberOfItems(); ++ridx)
 192       {
 193         DD(ridx);
 194         gdcm::Item & item = roi_seq->GetItem( ridx + 1); // Item starts at 1
 195
 196         const gdcm::Item & sitem = roi_seq->GetItem(ridx+1); // Item start at #1
 197
 198         const gdcm::DataSet& snestedds = sitem.GetNestedDataSet();
 199         const gdcm::DataSet& nestedds = item.GetNestedDataSet();
 200
 201         DD(nestedds.IsEmpty());
 202
 203         if( snestedds.FindDataElement( gdcm::Tag(0x3006,0x22) ) )
 204           {
 205             DD("tag found");
 206             // const gdcm::DataElement & a = nestedds.GetDataElement(gdcm::Tag(0x3006,0x26));
 207             // DD(a.GetValue());
 208
 209             gdcm::Attribute<0x3006,0x26> roiname;
 210             roiname.SetFromDataSet( snestedds );
 211             std::string name = roiname.GetValue();      // 0x3006,0x0026 = [ROI Name]
 212             DD(name);
 213             gdcm::Attribute<0x3006,0x0022> roinumber;
 214             roinumber.SetFromDataSet( snestedds );
 215             int nb = roinumber.GetValue();  // 0x3006,0x0022 = [ROI Number]
 216             DD(nb);
 217
 218             roi_names.push_back(make_pair(nb,name));
 219           }
 220       }
 221
 222     delete reader;
 223
 224   }
 225
 226 #else
 227   gdcm::File reader;
 228   reader.SetFileName(filename.c_str());
 229   reader.SetMaxSizeLoadEntry(16384);
 230   reader.Load();
 231
 232   gdcm::SeqEntry * roi_info=reader.GetSeqEntry(0x3006,0x0020);
 233   assert(roi_info);
 234   // DD("ici");
 235   //int n=0;
 236   for (gdcm::SQItem* i=roi_info->GetFirstSQItem(); i!=0; i=roi_info->GetNextSQItem())
 237     if (i->GetEntryValue(0x3006,0x0022)!= gdcm::GDCM_UNFOUND)
 238       roi_names.push_back(make_pair(atoi(i->GetEntryValue(0x3006,0x0022).c_str()),i->GetEntryValue(0x3006,0x0026)));
 239 #endif
 240   return roi_names;
 241 }
 242 //------------------------------------------------------------------------------
 243
 244
 245 //------------------------------------------------------------------------------
 246 std::vector<vvMesh::Pointer> vvMeshReader::readSelectedContours()
 247 {
 248   std::vector<vvMesh::Pointer> result;
 249 #if GDCM_MAJOR_VERSION == 2
 250   gdcm::Reader reader;
 251   reader.SetFileName(filename.c_str());
 252   reader.Read();
 253
 254   const gdcm::DataSet &ds = reader.GetFile().GetDataSet();
 255
 256   gdcm::SmartPointer<gdcm::SequenceOfItems> rois = ds.GetDataElement(gdcm::Tag(0x3006,0x39)).GetValueAsSQ();
 257   gdcm::SmartPointer<gdcm::SequenceOfItems> roi_info = ds.GetDataElement(gdcm::Tag(0x3006,0x20)).GetValueAsSQ();
 258   assert(rois); // TODO error message
 259   assert(roi_info); // TODO error message
 260   assert(rois->GetNumberOfItems() == roi_info->GetNumberOfItems());
 261
 262   for (unsigned ridx = 0; ridx < rois->GetNumberOfItems(); ++ridx)
 263     {
 264       vtkSmartPointer<vtkAppendPolyData> append=vtkSmartPointer<vtkAppendPolyData>::New();
 265       const gdcm::DataSet& ds_rois = rois->GetItem( ridx + 1).GetNestedDataSet();
 266       const gdcm::DataSet& ds_roi_info = roi_info->GetItem( ridx + 1).GetNestedDataSet();
 267
 268       gdcm::Attribute<0x3006,0x84> roinumber;
 269       roinumber.SetFromDataSet(ds_rois);
 270       if (std::find(selected_contours.begin(), selected_contours.end(), roinumber.GetValue()) != selected_contours.end()) //Only read selected ROIs
 271         {
 272           gdcm::Attribute<0x3006,0x2a> trgb;
 273           trgb.SetFromDataSet(ds_rois);
 274           vvMesh::Pointer current_roi=vvMesh::New();
 275           current_roi->r = trgb[0] / 255.0;
 276           current_roi->g = trgb[1] / 255.0;
 277           current_roi->b = trgb[2] / 255.0;
 278
 279           gdcm::Attribute<0x3006,0x26> tstructure_name;
 280           tstructure_name.SetFromDataSet(ds_roi_info);
 281           current_roi->structure_name = tstructure_name.GetValue();
 282
 283           gdcm::SmartPointer<gdcm::SequenceOfItems> roi_seq = ds_rois.GetDataElement(gdcm::Tag(0x3006,0x40)).GetValueAsSQ();
 284           double z0=-1; //Used to determine spacing between slices, assumed to be constant
 285           for (unsigned j = 0; j < roi_seq->GetNumberOfItems(); ++j)
 286             {
 287               gdcm::Item & item_roi_seq = roi_seq->GetItem(j + 1); // Item starts at 1
 288               const gdcm::DataSet& ds_roi_seq = item_roi_seq.GetNestedDataSet();
 289               gdcm::Attribute<0x3006,0x42> tcontour_type;
 290               tcontour_type.SetFromDataSet(ds_roi_seq);
 291               std::string contour_type = tcontour_type.GetValue();
 292               if (contour_type=="CLOSED_PLANAR ")
 293                 {
 294                   gdcm::Attribute<0x3006,0x46> tpoint_number;
 295                   tpoint_number.SetFromDataSet(ds_roi_seq);
 296                   const gdcm::DataElement & points_data = ds_roi_seq.GetDataElement(gdcm::Tag(0x3006,0x50));
 297                   gdcm::Attribute<0x3006,0x50> tpoints;
 298                   tpoints.SetFromDataElement(points_data);
 299                   assert(tpoints.GetNumberOfValues() == static_cast<unsigned int>(tpoint_number.GetValue()) * 3);
 300                   const double* points = tpoints.GetValues();
 301                   if (z0 == -1) //First contour
 302                     z0=points[2];
 303                   else
 304                     if (current_roi->GetSpacing()==-1 && points[2] != z0 )
 305                       current_roi->SetSpacing(points[2]-z0);
 306                   vtkPolyData * contour=vtkPolyData::New();
 307                   contour->Allocate(); //for cell structures
 308                   contour->SetPoints(vtkPoints::New());
 309                   vtkIdType ids[2];
 310                   for (unsigned idx = 0; idx < tpoints.GetNumberOfValues(); idx += 3)
 311                     {
 312                       contour->GetPoints()->InsertNextPoint(points[idx], points[idx+1], points[idx+2]);
 313                       ids[0] = idx / 3;
 314                       ids[1] = (ids[0] + 1) % tpoint_number.GetValue(); //0-1,1-2,...,n-1-0
 315                       contour->GetLines()->InsertNextCell(2, ids);
 316                     }
 317                   append->AddInput(contour);
 318                 }
 319               else
 320                 if (contour_type == "POINT ")
 321                   ; // silently ignore POINT type since we don't need them at the moment
 322                 else
 323                   std::cerr << "Warning: contour type " << contour_type << " not handled!" << std::endl;
 324             }
 325           append->Update();
 326           current_roi->AddMesh(append->GetOutput());
 327           result.push_back(current_roi);
 328         }
 329       else
 330         {
 331           //std::cerr << "Warning: ignoring ROI #" << roi_number << std::endl;
 332         }
 333     }
 334 #else
 335   gdcm::File reader;
 336   reader.SetFileName(filename.c_str());
 337   reader.SetMaxSizeLoadEntry(16384);
 338   reader.Load();
 339
 340   gdcm::SeqEntry * rois=reader.GetSeqEntry(0x3006,0x0039);
 341   ///We need to iterate both on the contours themselves, and on the contour info
 342   gdcm::SeqEntry * roi_info=reader.GetSeqEntry(0x3006,0x0020);
 343   gdcm::SQItem* k=roi_info->GetFirstSQItem();
 344   for(gdcm::SQItem* i=rois->GetFirstSQItem(); i!=0; i=rois->GetNextSQItem()) { //loop over ROIS
 345     assert(k!=0);
 346     vtkSmartPointer<vtkAppendPolyData> append=vtkSmartPointer<vtkAppendPolyData>::New();
 347     std::istringstream ss(i->GetEntryValue(0x3006,0x0084));
 348     int roi_number;
 349     ss >> roi_number;
 350     if (std::find(selected_contours.begin(),selected_contours.end(),roi_number) != selected_contours.end()) { //Only read selected ROIs
 351       vvMesh::Pointer current_roi=vvMesh::New();
 352       std::vector<double> rgb=clitk::parse_string<double>(i->GetEntryValue(0x3006,0x002a),'\\');
 353       assert(rgb.size()==3);
 354       current_roi->r=rgb[0]/255;
 355       current_roi->g=rgb[1]/255;
 356       current_roi->b=rgb[2]/255;
 357       current_roi->structure_name=k->GetEntryValue(0x3006,0x0026);
 358       gdcm::SeqEntry * contours=i->GetSeqEntry(0x3006,0x0040);
 359       double z0=-1; //Used to determine spacing between slices, assumed to be constant
 360       for(gdcm::SQItem* j=contours->GetFirstSQItem(); j!=0; j=contours->GetNextSQItem()) { //loop over 2D contours
 361         std::string contour_type=j->GetEntryValue(0x3006,0x0042);
 362         if (contour_type=="CLOSED_PLANAR ") {
 363           int point_number=clitk::parse_value<int>(j->GetEntryValue(0x3006,0x0046));
 364           std::vector<float> points=clitk::parse_string<float>(j->GetEntryValue(0x3006,0x0050),'\\');
 365           assert(points.size() == static_cast<unsigned int>(point_number)*3);
 366           if (z0 == -1) //First contour
 367             z0=points[2];
 368           //2nd contour, spacing not yet set. Need to be sure we are on a different slice,
 369           //sometimes there is more than one closed contour per slice
 370           else if (current_roi->GetSpacing()==-1 && points[2] != z0 )
 371             current_roi->SetSpacing(points[2]-z0);
 372           vtkPolyData * contour=vtkPolyData::New();
 373           contour->Allocate(); //for cell structures
 374           contour->SetPoints(vtkPoints::New());
 375           vtkIdType ids[2];
 376           for (unsigned int idx=0; idx<points.size(); idx+=3) {
 377             contour->GetPoints()->InsertNextPoint(points[idx],points[idx+1],points[idx+2]);
 378             ids[0]=idx/3;
 379             ids[1]=(ids[0]+1)%point_number; //0-1,1-2,...,n-1-0
 380             contour->GetLines()->InsertNextCell(2,ids);
 381           }
 382           append->AddInput(contour);
 383         } else if (contour_type == "POINT ")
 384           ; // silently ignore POINT type since we don't need them at the moment
 385         else
 386           std::cerr << "Warning: contour type " << contour_type << " not handled!" << std::endl;
 387       }
 388       append->Update();
 389       current_roi->AddMesh(append->GetOutput());
 390       result.push_back(current_roi);
 391     } else {
 392       //std::cerr << "Warning: ignoring ROI #" << roi_number << std::endl;
 393     }
 394     k=roi_info->GetNextSQItem(); //increment the second loop variable
 395   }
 396 #endif
 397   return result;
 398 }
 399