common/vtkGDCMPolyDataReader.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
  20   Program: GDCM (Grassroots DICOM). A DICOM library
  21
  22   Copyright (c) 2006-2011 Mathieu Malaterre
  23   All rights reserved.
  24   See Copyright.txt or http://gdcm.sourceforge.net/Copyright.html for details.
  25
  26      This software is distributed WITHOUT ANY WARRANTY; without even
  27      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  28      PURPOSE.  See the above copyright notice for more information.
  29
  30 =========================================================================*/
  31 #include "vtkGDCMPolyDataReader.h"
  32
  33 #include "vtkObjectFactory.h"
  34 #include "vtkInformation.h"
  35 #include "vtkInformationVector.h"
  36 #include "vtkPolyData.h"
  37 #include "vtkStreamingDemandDrivenPipeline.h"
  38 #include "vtkDoubleArray.h"
  39 #include "vtkCellArray.h"
  40 #include "vtkCellData.h"
  41 #include "vtkMedicalImageProperties.h"
  42 #include "vtkRTStructSetProperties.h"
  43
  44 #include "gdcmReader.h"
  45 #include "gdcmSmartPointer.h"
  46 #include "gdcmAttribute.h"
  47 #include "gdcmSequenceOfItems.h"
  48
  49 vtkCxxRevisionMacro(vtkGDCMPolyDataReader, "$Revision: 1.1 $")
  50 vtkStandardNewMacro(vtkGDCMPolyDataReader)
  51
  52 //----------------------------------------------------------------------------
  53 vtkGDCMPolyDataReader::vtkGDCMPolyDataReader()
  54 {
  55   this->FileName = NULL;
  56   this->SetNumberOfInputPorts(0);
  57   this->MedicalImageProperties = vtkMedicalImageProperties::New();
  58   this->RTStructSetProperties = vtkRTStructSetProperties::New();
  59 }
  60
  61 //----------------------------------------------------------------------------
  62 vtkGDCMPolyDataReader::~vtkGDCMPolyDataReader()
  63 {
  64   this->SetFileName(0);
  65   this->MedicalImageProperties->Delete();
  66   this->RTStructSetProperties->Delete();
  67 }
  68
  69 //----------------------------------------------------------------------------
  70 // inline keyword is needed since GetStringValueFromTag was copy/paste from vtkGDCMImageReader
  71 // FIXME: need to restructure code to avoid copy/paste
  72 inline const char *GetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds)
  73 {
  74   static std::string buffer;
  75   buffer = "";  // cleanup previous call
  76
  77   if( ds.FindDataElement( t ) )
  78     {
  79     const gdcm::DataElement& de = ds.GetDataElement( t );
  80     const gdcm::ByteValue *bv = de.GetByteValue();
  81     if( bv ) // Can be Type 2
  82       {
  83       buffer = std::string( bv->GetPointer(), bv->GetLength() );
  84       // Will be padded with at least one \0
  85       }
  86     }
  87
  88   // Since return is a const char* the very first \0 will be considered
  89   return buffer.c_str();
  90 }
  91
  92 //----------------------------------------------------------------------------
  93 void vtkGDCMPolyDataReader::FillMedicalImageInformation(const gdcm::Reader &reader)
  94 {
  95   const gdcm::File &file = reader.GetFile();
  96   const gdcm::DataSet &ds = file.GetDataSet();
  97
  98   this->RTStructSetProperties->SetStructureSetLabel( GetStringValueFromTag( gdcm::Tag(0x3006,0x0002), ds) );
  99   this->RTStructSetProperties->SetStructureSetName( GetStringValueFromTag( gdcm::Tag(0x3006,0x0004), ds) );
 100   this->RTStructSetProperties->SetStructureSetDate( GetStringValueFromTag( gdcm::Tag(0x3006,0x0008), ds) );
 101   this->RTStructSetProperties->SetStructureSetTime( GetStringValueFromTag( gdcm::Tag(0x3006,0x0009), ds) );
 102   this->RTStructSetProperties->SetSOPInstanceUID( GetStringValueFromTag( gdcm::Tag(0x0008,0x0018), ds) );
 103   this->RTStructSetProperties->SetStudyInstanceUID( GetStringValueFromTag( gdcm::Tag(0x0020,0x000d), ds) );
 104   this->RTStructSetProperties->SetSeriesInstanceUID( GetStringValueFromTag( gdcm::Tag(0x0020,0x000e), ds) );
 105
 106   // $ grep "vtkSetString\|DICOM" vtkMedicalImageProperties.h
 107   // For ex: DICOM (0010,0010) = DOE,JOHN
 108   this->MedicalImageProperties->SetPatientName( GetStringValueFromTag( gdcm::Tag(0x0010,0x0010), ds) );
 109   // For ex: DICOM (0010,0020) = 1933197
 110   this->MedicalImageProperties->SetPatientID( GetStringValueFromTag( gdcm::Tag(0x0010,0x0020), ds) );
 111   // For ex: DICOM (0010,1010) = 031Y
 112   this->MedicalImageProperties->SetPatientAge( GetStringValueFromTag( gdcm::Tag(0x0010,0x1010), ds) );
 113   // For ex: DICOM (0010,0040) = M
 114   this->MedicalImageProperties->SetPatientSex( GetStringValueFromTag( gdcm::Tag(0x0010,0x0040), ds) );
 115   // For ex: DICOM (0010,0030) = 19680427
 116   this->MedicalImageProperties->SetPatientBirthDate( GetStringValueFromTag( gdcm::Tag(0x0010,0x0030), ds) );
 117 #if ( VTK_MAJOR_VERSION == 5 && VTK_MINOR_VERSION > 0 )
 118   // For ex: DICOM (0008,0020) = 20030617
 119   this->MedicalImageProperties->SetStudyDate( GetStringValueFromTag( gdcm::Tag(0x0008,0x0020), ds) );
 120 #endif
 121   // For ex: DICOM (0008,0022) = 20030617
 122   this->MedicalImageProperties->SetAcquisitionDate( GetStringValueFromTag( gdcm::Tag(0x0008,0x0022), ds) );
 123 #if ( VTK_MAJOR_VERSION == 5 && VTK_MINOR_VERSION > 0 )
 124   // For ex: DICOM (0008,0030) = 162552.0705 or 230012, or 0012
 125   this->MedicalImageProperties->SetStudyTime( GetStringValueFromTag( gdcm::Tag(0x0008,0x0030), ds) );
 126 #endif
 127   // For ex: DICOM (0008,0032) = 162552.0705 or 230012, or 0012
 128   this->MedicalImageProperties->SetAcquisitionTime( GetStringValueFromTag( gdcm::Tag(0x0008,0x0032), ds) );
 129   // For ex: DICOM (0008,0023) = 20030617
 130   this->MedicalImageProperties->SetImageDate( GetStringValueFromTag( gdcm::Tag(0x0008,0x0023), ds) );
 131   // For ex: DICOM (0008,0033) = 162552.0705 or 230012, or 0012
 132   this->MedicalImageProperties->SetImageTime( GetStringValueFromTag( gdcm::Tag(0x0008,0x0033), ds) );
 133   // For ex: DICOM (0020,0013) = 1
 134   this->MedicalImageProperties->SetImageNumber( GetStringValueFromTag( gdcm::Tag(0x0020,0x0013), ds) );
 135   // For ex: DICOM (0020,0011) = 902
 136   this->MedicalImageProperties->SetSeriesNumber( GetStringValueFromTag( gdcm::Tag(0x0020,0x0011), ds) );
 137   // For ex: DICOM (0008,103e) = SCOUT
 138   this->MedicalImageProperties->SetSeriesDescription( GetStringValueFromTag( gdcm::Tag(0x0008,0x103e), ds) );
 139   // For ex: DICOM (0020,0010) = 37481
 140   this->MedicalImageProperties->SetStudyID( GetStringValueFromTag( gdcm::Tag(0x0020,0x0010), ds) );
 141   // For ex: DICOM (0008,1030) = BRAIN/C-SP/FACIAL
 142   this->MedicalImageProperties->SetStudyDescription( GetStringValueFromTag( gdcm::Tag(0x0008,0x1030), ds) );
 143   // For ex: DICOM (0008,0060)= CT
 144   this->MedicalImageProperties->SetModality( GetStringValueFromTag( gdcm::Tag(0x0008,0x0060), ds) );
 145   // For ex: DICOM (0008,0070) = Siemens
 146   this->MedicalImageProperties->SetManufacturer( GetStringValueFromTag( gdcm::Tag(0x0008,0x0070), ds) );
 147   // For ex: DICOM (0008,1090) = LightSpeed QX/i
 148   this->MedicalImageProperties->SetManufacturerModelName( GetStringValueFromTag( gdcm::Tag(0x0008,0x1090), ds) );
 149   // For ex: DICOM (0008,1010) = LSPD_OC8
 150   this->MedicalImageProperties->SetStationName( GetStringValueFromTag( gdcm::Tag(0x0008,0x1010), ds) );
 151   // For ex: DICOM (0008,0080) = FooCity Medical Center
 152   this->MedicalImageProperties->SetInstitutionName( GetStringValueFromTag( gdcm::Tag(0x0008,0x0080), ds) );
 153   // For ex: DICOM (0018,1210) = Bone
 154   this->MedicalImageProperties->SetConvolutionKernel( GetStringValueFromTag( gdcm::Tag(0x0018,0x1210), ds) );
 155   // For ex: DICOM (0018,0050) = 0.273438
 156   this->MedicalImageProperties->SetSliceThickness( GetStringValueFromTag( gdcm::Tag(0x0018,0x0050), ds) );
 157   // For ex: DICOM (0018,0060) = 120
 158   this->MedicalImageProperties->SetKVP( GetStringValueFromTag( gdcm::Tag(0x0018,0x0060), ds) );
 159   // For ex: DICOM (0018,1120) = 15
 160   this->MedicalImageProperties->SetGantryTilt( GetStringValueFromTag( gdcm::Tag(0x0018,0x1120), ds) );
 161   // For ex: DICOM (0018,0081) = 105
 162   this->MedicalImageProperties->SetEchoTime( GetStringValueFromTag( gdcm::Tag(0x0018,0x0081), ds) );
 163   // For ex: DICOM (0018,0091) = 35
 164   this->MedicalImageProperties->SetEchoTrainLength( GetStringValueFromTag( gdcm::Tag(0x0018,0x0091), ds) );
 165   // For ex: DICOM (0018,0080) = 2040
 166   this->MedicalImageProperties->SetRepetitionTime( GetStringValueFromTag( gdcm::Tag(0x0018,0x0080), ds) );
 167   // For ex: DICOM (0018,1150) = 5
 168   this->MedicalImageProperties->SetExposureTime( GetStringValueFromTag( gdcm::Tag(0x0018,0x1150), ds) );
 169   // For ex: DICOM (0018,1151) = 400
 170   this->MedicalImageProperties->SetXRayTubeCurrent( GetStringValueFromTag( gdcm::Tag(0x0018,0x1151), ds) );
 171   // For ex: DICOM (0018,1152) = 114
 172   this->MedicalImageProperties->SetExposure( GetStringValueFromTag( gdcm::Tag(0x0018,0x1152), ds) );
 173
 174   // virtual void AddWindowLevelPreset(double w, double l);
 175   // (0028,1050) DS [   498\  498]                           #  12, 2 WindowCenter
 176   // (0028,1051) DS [  1063\ 1063]                           #  12, 2 WindowWidth
 177   gdcm::Tag twindowcenter(0x0028,0x1050);
 178   gdcm::Tag twindowwidth(0x0028,0x1051);
 179   if( ds.FindDataElement( twindowcenter ) && ds.FindDataElement( twindowwidth) )
 180     {
 181     const gdcm::DataElement& windowcenter = ds.GetDataElement( twindowcenter );
 182     const gdcm::DataElement& windowwidth = ds.GetDataElement( twindowwidth );
 183     const gdcm::ByteValue *bvwc = windowcenter.GetByteValue();
 184     const gdcm::ByteValue *bvww = windowwidth.GetByteValue();
 185     if( bvwc && bvww ) // Can be Type 2
 186       {
 187       //gdcm::Attributes<0x0028,0x1050> at;
 188       gdcm::Element<gdcm::VR::DS,gdcm::VM::VM1_n> elwc;
 189       std::stringstream ss1;
 190       std::string swc = std::string( bvwc->GetPointer(), bvwc->GetLength() );
 191       ss1.str( swc );
 192       gdcm::VR vr = gdcm::VR::DS;
 193       unsigned int vrsize = vr.GetSizeof();
 194       unsigned int count = gdcm::VM::GetNumberOfElementsFromArray(swc.c_str(), swc.size());
 195       elwc.SetLength( count * vrsize );
 196       elwc.Read( ss1 );
 197       std::stringstream ss2;
 198       std::string sww = std::string( bvww->GetPointer(), bvww->GetLength() );
 199       ss2.str( sww );
 200       gdcm::Element<gdcm::VR::DS,gdcm::VM::VM1_n> elww;
 201       elww.SetLength( count * vrsize );
 202       elww.Read( ss2 );
 203       //assert( elww.GetLength() == elwc.GetLength() );
 204       for(unsigned int i = 0; i < elwc.GetLength(); ++i)
 205         {
 206         this->MedicalImageProperties->AddWindowLevelPreset( elww.GetValue(i), elwc.GetValue(i) );
 207         }
 208       }
 209     }
 210   gdcm::Tag twindowexplanation(0x0028,0x1055);
 211   if( ds.FindDataElement( twindowexplanation ) )
 212     {
 213     const gdcm::DataElement& windowexplanation = ds.GetDataElement( twindowexplanation );
 214     const gdcm::ByteValue *bvwe = windowexplanation.GetByteValue();
 215     if( bvwe ) // Can be Type 2
 216       {
 217       int n = this->MedicalImageProperties->GetNumberOfWindowLevelPresets();
 218       gdcm::Element<gdcm::VR::LO,gdcm::VM::VM1_n> elwe; // window explanation
 219       gdcm::VR vr = gdcm::VR::LO;
 220       std::stringstream ss;
 221       ss.str( "" );
 222       std::string swe = std::string( bvwe->GetPointer(), bvwe->GetLength() );
 223       unsigned int count = gdcm::VM::GetNumberOfElementsFromArray(swe.c_str(), swe.size()); (void)count;
 224       // I found a case with only one W/L but two comments: WINDOW1\WINDOW2
 225       // SIEMENS-IncompletePixelData.dcm
 226       //assert( count >= (unsigned int)n );
 227       elwe.SetLength( /*count*/ n * vr.GetSizeof() );
 228       ss.str( swe );
 229       elwe.Read( ss );
 230       for(int i = 0; i < n; ++i)
 231         {
 232         this->MedicalImageProperties->SetNthWindowLevelPresetComment(i, elwe.GetValue(i).c_str() );
 233         }
 234       }
 235     }
 236
 237 #if 0
 238   // gdcmData/JDDICOM_Sample4.dcm
 239   // -> (0008,0060) CS [DM  Digital microscopy]                 #  24, 1 Modality
 240   gdcm::MediaStorage ms1 = gdcm::MediaStorage::SecondaryCaptureImageStorage;
 241   ms1.GuessFromModality( this->MedicalImageProperties->GetModality(), this->FileDimensionality );
 242   gdcm::MediaStorage ms2;
 243   ms2.SetFromFile( reader.GetFile() );
 244   if( ms2 != ms1 && ms2 != gdcm::MediaStorage::SecondaryCaptureImageStorage )
 245     {
 246     vtkWarningMacro( "SHOULD NOT HAPPEN. Unrecognized Modality: " << this->MedicalImageProperties->GetModality()
 247       << " Will be set instead to the known one: " << ms2.GetModality() )
 248     this->MedicalImageProperties->SetModality( ms2.GetModality() );
 249     }
 250 #endif
 251
 252   // Add more info:
 253
 254 }
 255
 256 //(3006,0022) ?? (IS) [2 ]                                      # 2,1 ROI Number
 257 //(3006,0024) ?? (UI) [2.16.840.1.114362.1.759508.1251415878280.193]         # 44,1 Referenced Frame of Reference UID
 258 //(3006,0026) ?? (LO) [Bladder ]                                # 8,1 ROI Name
 259 //(3006,0036) ?? (CS) [MANUAL]                                  # 6,1 ROI Generation Algorithm
 260
 261 int vtkGDCMPolyDataReader::RequestData_RTStructureSetStorage(gdcm::Reader const &reader,
 262   vtkInformationVector *outputVector)
 263 {
 264 // This is done in RequestInformation
 265 //  gdcm::MediaStorage ms;
 266 //  ms.SetFromFile( reader.GetFile() );
 267 //  //std::cout << ms << std::endl;
 268 //  if( ms != gdcm::MediaStorage::RTStructureSetStorage )
 269 //    {
 270 //    return 0;
 271 //    }
 272
 273   const gdcm::DataSet& ds = reader.GetFile().GetDataSet();
 274   // (3006,0010) SQ (Sequence with undefined length #=1)     # u/l, 1 ReferencedFrameOfReferenceSequence
 275   // (3006,0020) SQ (Sequence with explicit length #=4)      # 370, 1 StructureSetROISequence
 276   // (3006,0039) SQ (Sequence with explicit length #=4)      # 24216, 1 ROIContourSequence
 277   gdcm::Tag troicsq(0x3006,0x0039);
 278   if( !ds.FindDataElement( troicsq ) )
 279     {
 280     return 0;
 281     }
 282   gdcm::Tag tssroisq(0x3006,0x0020);
 283   if( !ds.FindDataElement( tssroisq ) )
 284     {
 285     return 0;
 286     }
 287   gdcm::Tag trefframerefsq(0x3006,0x0010);
 288   if( !ds.FindDataElement( trefframerefsq ) )
 289     {
 290     return 0;
 291     }
 292   const gdcm::DataElement &refframerefsq = ds.GetDataElement( trefframerefsq );
 293   gdcm::SmartPointer<gdcm::SequenceOfItems> sqi0 = refframerefsq.GetValueAsSQ();
 294   if( !sqi0 || !sqi0->GetNumberOfItems() )
 295     {
 296     return 0;
 297     }
 298     assert( sqi0->GetNumberOfItems() == 1 );
 299   for(unsigned int pd = 0; pd < sqi0->GetNumberOfItems(); ++pd)
 300     {
 301     const gdcm::Item & item = sqi0->GetItem(pd+1); // Item start at #1
 302     const gdcm::DataSet& nestedds = item.GetNestedDataSet();
 303     // (3006,0012) SQ (Sequence with undefined length #=1)     # u/l, 1 RTReferencedStudySequence
 304     gdcm::Attribute<0x0020,0x052> frameofreferenceuid;
 305     frameofreferenceuid.SetFromDataSet( nestedds );
 306     this->RTStructSetProperties->SetReferenceFrameOfReferenceUID(
 307       frameofreferenceuid.GetValue() );
 308     gdcm::Tag trtrefstudysq(0x3006,0x0012);
 309     if( !nestedds.FindDataElement( trtrefstudysq) )
 310       {
 311       return 0;
 312       }
 313     const gdcm::DataElement &rtrefstudysq = nestedds.GetDataElement( trtrefstudysq );
 314     gdcm::SmartPointer<gdcm::SequenceOfItems> sqi00 = rtrefstudysq.GetValueAsSQ();
 315     if( !sqi00 || !sqi00->GetNumberOfItems() )
 316       {
 317       return 0;
 318       }
 319     assert( sqi00->GetNumberOfItems() == 1 );
 320     for(unsigned int pd0 = 0; pd0 < sqi00->GetNumberOfItems(); ++pd0)
 321       {
 322       const gdcm::Item & item = sqi00->GetItem(pd0+1); // Item start at #1
 323       const gdcm::DataSet& nestedds = item.GetNestedDataSet();
 324
 325         // (3006,0014) SQ (Sequence with undefined length #=1)     # u/l, 1 RTReferencedSeriesSequence
 326     gdcm::Tag trtrefseriessq(0x3006,0x0014);
 327     if( !nestedds.FindDataElement( trtrefseriessq) )
 328       {
 329       return 0;
 330       }
 331     const gdcm::DataElement &rtrefseriessq = nestedds.GetDataElement( trtrefseriessq);
 332
 333     gdcm::SmartPointer<gdcm::SequenceOfItems> sqi000 = rtrefseriessq.GetValueAsSQ();
 334     if( !sqi000 || !sqi000->GetNumberOfItems() )
 335       {
 336       return 0;
 337       }
 338     assert( sqi000->GetNumberOfItems() == 1 );
 339     for(unsigned int pd00 = 0; pd00 < sqi000->GetNumberOfItems(); ++pd00)
 340       {
 341       const gdcm::Item & item = sqi000->GetItem(pd00+1); // Item start at #1
 342       const gdcm::DataSet& nestedds = item.GetNestedDataSet();
 343
 344     gdcm::Attribute<0x0020,0x000e> seriesinstanceuid;
 345     seriesinstanceuid.SetFromDataSet( nestedds );
 346     this->RTStructSetProperties->SetReferenceSeriesInstanceUID(
 347       seriesinstanceuid.GetValue() );
 348
 349
 350             // (3006,0016) SQ (Sequence with undefined length #=162)   # u/l, 1 ContourImageSequence
 351     gdcm::Tag tcontourimageseq(0x3006,0x0016);
 352     if( !nestedds.FindDataElement( tcontourimageseq) )
 353       {
 354       return 0;
 355       }
 356     const gdcm::DataElement &contourimageseq = nestedds.GetDataElement( tcontourimageseq );
 357     gdcm::SmartPointer<gdcm::SequenceOfItems> sqi0000 = contourimageseq.GetValueAsSQ();
 358     if( !sqi0000 || !sqi0000->GetNumberOfItems() )
 359       {
 360       return 0;
 361       }
 362
 363     assert( sqi0000->GetNumberOfItems() != 1 );
 364     for(unsigned int pd000 = 0; pd000 < sqi0000->GetNumberOfItems(); ++pd000)
 365       {
 366       const gdcm::Item & item = sqi0000->GetItem(pd000+1); // Item start at #1
 367       const gdcm::DataSet& nestedds = item.GetNestedDataSet();
 368       gdcm::Attribute<0x0008,0x1150> refsopclassuid;
 369       refsopclassuid.SetFromDataSet( nestedds );
 370       gdcm::Attribute<0x0008,0x1155> refinstanceuid;
 371       refinstanceuid.SetFromDataSet( nestedds );
 372   this->RTStructSetProperties->AddReferencedFrameOfReference( refsopclassuid.GetValue().c_str(),
 373 refinstanceuid.GetValue().c_str() );
 374       }
 375
 376       }
 377
 378       }
 379     }
 380
 381   const gdcm::DataElement &roicsq = ds.GetDataElement( troicsq );
 382   //std::cout << roicsq << std::endl;
 383   //const gdcm::SequenceOfItems *sqi_debug = roicsq.GetSequenceOfItems();
 384   gdcm::SmartPointer<gdcm::SequenceOfItems> sqi = roicsq.GetValueAsSQ();
 385   if( !sqi || !sqi->GetNumberOfItems() )
 386     {
 387     return 0;
 388     }
 389   const gdcm::DataElement &ssroisq = ds.GetDataElement( tssroisq );
 390   //const gdcm::SequenceOfItems *ssqi = ssroisq.GetSequenceOfItems();
 391   gdcm::SmartPointer<gdcm::SequenceOfItems> ssqi = ssroisq.GetValueAsSQ();
 392   if( !ssqi || !ssqi->GetNumberOfItems() )
 393     {
 394     return 0;
 395     }
 396
 397   // For each Item in the DataSet create a vtkPolyData
 398   for(unsigned int pd = 0; pd < sqi->GetNumberOfItems(); ++pd)
 399     {
 400     //StructureSetROI structuresetroi;
 401     // get the info object
 402     vtkInformation *outInfo1 = outputVector->GetInformationObject(pd);
 403
 404     // get the ouptut
 405     vtkPolyData *output = vtkPolyData::SafeDownCast(
 406       outInfo1->Get(vtkDataObject::DATA_OBJECT()));
 407
 408
 409     const gdcm::Item & item = sqi->GetItem(pd+1); // Item start at #1
 410     //std::cout << item << std::endl;
 411     const gdcm::Item & sitem = ssqi->GetItem(pd+1); // Item start at #1
 412     const gdcm::DataSet& snestedds = sitem.GetNestedDataSet();
 413     // (3006,0026) ?? (LO) [date]                                    # 4,1 ROI Name
 414     gdcm::Tag stcsq(0x3006,0x0026);
 415     if( !snestedds.FindDataElement( stcsq ) )
 416       {
 417       continue;
 418       }
 419     const gdcm::DataElement &sde = snestedds.GetDataElement( stcsq );
 420     std::string s(sde.GetByteValue()->GetPointer(), sde.GetByteValue()->GetLength());
 421     //structuresetroi.ROIName = s;
 422     gdcm::Attribute<0x3006,0x0022> roinumber;
 423     roinumber.SetFromDataSet( snestedds );
 424     //structuresetroi.ROINumber = roinumber.GetValue();
 425     gdcm::Attribute<0x3006,0x0024> refframeuid;
 426     refframeuid.SetFromDataSet( snestedds );
 427     //structuresetroi.RefFrameRefUID = refframeuid.GetValue();
 428     gdcm::Attribute<0x3006,0x0026> roiname;
 429     roiname.SetFromDataSet( snestedds );
 430     assert( s == roiname.GetValue() );
 431     gdcm::Attribute<0x3006,0x0036> roigenalg;
 432     roigenalg.SetFromDataSet( snestedds );
 433     //structuresetroi.ROIGenerationAlgorithm = roigenalg.GetValue();
 434     //structuresetrois.push_back( structuresetroi );
 435
 436     this->RTStructSetProperties->AddStructureSetROI(
 437       roinumber.GetValue(),
 438       refframeuid.GetValue(),
 439       roiname.GetValue(),
 440       roigenalg.GetValue()
 441     );
 442
 443     const gdcm::DataSet& nestedds = item.GetNestedDataSet();
 444     //std::cout << nestedds << std::endl;
 445     //(3006,002a) IS [255\192\96]                              # 10,3 ROI Display Color
 446     gdcm::Tag troidc(0x3006,0x002a);
 447     gdcm::Attribute<0x3006,0x002a> color = {};
 448     if( nestedds.FindDataElement( troidc) )
 449       {
 450       const gdcm::DataElement &decolor = nestedds.GetDataElement( troidc );
 451       color.SetFromDataElement( decolor );
 452       //std::cout << "color:" << color[0] << "," << color[1] << "," << color[2] << std::endl;
 453       }
 454     //(3006,0040) SQ (Sequence with explicit length #=8)      # 4326, 1 ContourSequence
 455     gdcm::Tag tcsq(0x3006,0x0040);
 456     if( !nestedds.FindDataElement( tcsq ) )
 457       {
 458       continue;
 459       }
 460     const gdcm::DataElement& csq = nestedds.GetDataElement( tcsq );
 461     //std::cout << csq << std::endl;
 462
 463     //const gdcm::SequenceOfItems *sqi2 = csq.GetSequenceOfItems();
 464     gdcm::SmartPointer<gdcm::SequenceOfItems> sqi2 = csq.GetValueAsSQ();
 465     if( !sqi2 || !sqi2->GetNumberOfItems() )
 466       {
 467       continue;
 468       }
 469     unsigned int nitems = sqi2->GetNumberOfItems();
 470     //std::cout << nitems << std::endl;
 471     //this->SetNumberOfOutputPorts(nitems);
 472     vtkDoubleArray *scalars = vtkDoubleArray::New();
 473     scalars->SetNumberOfComponents(3);
 474
 475     vtkPoints *newPts = vtkPoints::New();
 476     //std::string s(sde.GetByteValue()->GetPointer(), sde.GetByteValue()->GetLength());
 477     //std::cout << s << std::endl;
 478     //newPts->GetData()->SetName( s.c_str() );
 479     // In VTK there is no API to specify the name of a vtkPolyData, you can only specify Name
 480     // for the scalars (pointdata or celldata), so let's do that...
 481     //scalars->SetName( structuresetroi.ROIName.c_str() );
 482     scalars->SetName( roiname.GetValue().c_str() );
 483     vtkCellArray *polys = vtkCellArray::New();
 484     for(unsigned int i = 0; i < nitems; ++i)
 485       {
 486       const gdcm::Item & item2 = sqi2->GetItem(i+1); // Item start at #1
 487
 488       const gdcm::DataSet& nestedds2 = item2.GetNestedDataSet();
 489       //std::cout << nestedds2 << std::endl;
 490       // (3006,0050) DS [43.57636\65.52504\-10.0\46.043102\62.564945\-10.0\49.126537\60.714... # 398,48 ContourData
 491       gdcm::Tag tcontourdata(0x3006,0x0050);
 492       const gdcm::DataElement & contourdata = nestedds2.GetDataElement( tcontourdata );
 493       //std::cout << contourdata << std::endl;
 494
 495       //const gdcm::ByteValue *bv = contourdata.GetByteValue();
 496       gdcm::Attribute<0x3006,0x0042> contgeotype;
 497       contgeotype.SetFromDataSet( nestedds2 );
 498       assert( contgeotype.GetValue() == "CLOSED_PLANAR " );
 499
 500       gdcm::Attribute<0x3006,0x0046> numcontpoints;
 501       numcontpoints.SetFromDataSet( nestedds2 );
 502
 503       gdcm::Attribute<0x3006,0x0050> at;
 504       at.SetFromDataElement( contourdata );
 505
 506         {
 507         assert( nestedds2.FindDataElement( gdcm::Tag(0x3006,0x0016) ) );
 508         const gdcm::DataElement &contourimagesequence = nestedds2.GetDataElement( gdcm::Tag(0x3006,0x0016) );
 509         gdcm::SmartPointer<gdcm::SequenceOfItems> contourimagesequence_sqi = contourimagesequence.GetValueAsSQ();
 510         assert( contourimagesequence_sqi && contourimagesequence_sqi->GetNumberOfItems() == 1 );
 511       const gdcm::Item & theitem = contourimagesequence_sqi->GetItem(1);
 512       const gdcm::DataSet& nestedds = theitem.GetNestedDataSet();
 513
 514       gdcm::Attribute<0x0008,0x1150> classat;
 515       classat.SetFromDataSet( nestedds );
 516       gdcm::Attribute<0x0008,0x1155> instat;
 517       instat.SetFromDataSet( nestedds );
 518
 519       this->RTStructSetProperties->AddContourReferencedFrameOfReference( pd,
 520         classat.GetValue(), instat.GetValue() );
 521         }
 522
 523       //newPts->SetNumberOfPoints( at.GetNumberOfValues() / 3 );
 524       //assert( at.GetNumberOfValues() % 3 == 0); // FIXME
 525       const double* pts = at.GetValues();
 526       vtkIdType buffer[256];
 527       vtkIdType *ptIds;
 528       unsigned int npts = at.GetNumberOfValues() / 3;
 529       assert( npts == numcontpoints.GetValue() );
 530       if(npts>256)
 531         {
 532         ptIds = new vtkIdType[npts];
 533         }
 534       else
 535         {
 536         ptIds = buffer;
 537         }
 538       for(unsigned int i = 0; i < npts * 3; i+=3)
 539         {
 540         float x[3];
 541         x[0] = pts[i+0];
 542         x[1] = pts[i+1];
 543         x[2] = pts[i+2];
 544         vtkIdType ptId = newPts->InsertNextPoint( x );
 545         ptIds[i / 3] = ptId;
 546         }
 547       // Each Contour Data is in fact a Cell:
 548       vtkIdType cellId = polys->InsertNextCell( npts , ptIds);
 549       scalars->InsertTuple3(cellId, (double)color[0]/255.0, (double)color[1]/255.0, (double)color[2]/255.0);
 550       if(npts>256)
 551         {
 552         delete[] ptIds;
 553         }
 554       }
 555     output->SetPoints(newPts);
 556     newPts->Delete();
 557     output->SetPolys(polys);
 558     polys->Delete();
 559     output->GetCellData()->SetScalars(scalars);
 560     scalars->Delete();
 561     }
 562
 563   // Add the Observations:
 564   // we can only be doing it here once all RT are loaded, since we will
 565   // attach observation to *existing* rtstruct
 566   gdcm::Tag trtroiobssq(0x3006,0x0080);
 567   if( !ds.FindDataElement( trtroiobssq ) )
 568     {
 569     return 0;
 570     }
 571   const gdcm::DataElement &rtroiobssq = ds.GetDataElement( trtroiobssq );
 572   gdcm::SmartPointer<gdcm::SequenceOfItems> rtroiobssqsqi = rtroiobssq.GetValueAsSQ();
 573   if( !rtroiobssqsqi || !rtroiobssqsqi->GetNumberOfItems() )
 574     {
 575     return 0;
 576     }
 577   for(unsigned int obs = 0; obs < rtroiobssqsqi->GetNumberOfItems(); ++obs)
 578     {
 579     const gdcm::Item & item = rtroiobssqsqi->GetItem(obs+1); // Item start at #1
 580     const gdcm::DataSet& nestedds = item.GetNestedDataSet();
 581     gdcm::Attribute<0x3006,0x0082> observationnumber;
 582     observationnumber.SetFromDataSet( nestedds );
 583     gdcm::Attribute<0x3006,0x0084> referencedroinumber;
 584     referencedroinumber.SetFromDataSet( nestedds );
 585     gdcm::Attribute<0x3006,0x00a4> rtroiinterpretedtype;
 586     rtroiinterpretedtype.SetFromDataSet( nestedds );
 587     gdcm::Attribute<0x3006,0x00a6> roiinterpreter;
 588     roiinterpreter.SetFromDataSet( nestedds );
 589     this->RTStructSetProperties->
 590       AddStructureSetROIObservation( referencedroinumber.GetValue(),
 591         observationnumber.GetValue(),
 592         rtroiinterpretedtype.GetValue(),
 593         roiinterpreter.GetValue() );
 594     }
 595
 596   return 1;
 597 }
 598
 599 int vtkGDCMPolyDataReader::RequestData_HemodynamicWaveformStorage(gdcm::Reader const &reader,
 600   vtkInformationVector *outputVector)
 601 {
 602   const gdcm::DataSet& ds = reader.GetFile().GetDataSet();
 603   // (5400,0100) SQ (Sequence with undefined length #=1)     # u/l, 1 WaveformSequence
 604   gdcm::Tag twsq(0x5400,0x0100);
 605   if( !ds.FindDataElement( twsq) )
 606     {
 607     return 0;
 608     }
 609   const gdcm::DataElement &wsq = ds.GetDataElement( twsq );
 610   //std::cout << wsq << std::endl;
 611   //const gdcm::SequenceOfItems *sqi = wsq.GetSequenceOfItems();
 612   gdcm::SmartPointer<gdcm::SequenceOfItems> sqi = wsq.GetValueAsSQ();
 613   if( !sqi || !sqi->GetNumberOfItems() )
 614     {
 615     return 0;
 616     }
 617
 618   const gdcm::Item & item = sqi->GetItem(1); // Item start at #1
 619   const gdcm::DataSet& nestedds = item.GetNestedDataSet();
 620
 621   // (5400,1004) US 16                                       #   2, 1 WaveformBitsAllocated
 622   gdcm::Tag twba(0x5400,0x1004);
 623   if( !nestedds.FindDataElement( twba ) )
 624     {
 625     return 0;
 626     }
 627   const gdcm::DataElement &wba= nestedds.GetDataElement( twba );
 628
 629   //std::cout << wba << std::endl;
 630   //  (5400,1006) CS [SS]                                     #   2, 1 WaveformSampleInterpretation
 631   // (5400,1010) OW 00ba\0030\ff76\ff8b\00a2\ffd3\ffae\ff50\0062\00c4\011e\00c2\00ba... # 57600, 1 WaveformData
 632   gdcm::Tag twd(0x5400,0x1010);
 633   if( !nestedds.FindDataElement( twd ) )
 634     {
 635     return 0;
 636     }
 637   const gdcm::DataElement &wd = nestedds.GetDataElement( twd );
 638   const gdcm::ByteValue *bv = wd.GetByteValue();
 639   size_t len = bv->GetLength();
 640   int16_t *p = (int16_t*)bv;
 641
 642   // get the info object
 643   int pd = 0;
 644   vtkInformation *outInfo1 = outputVector->GetInformationObject(pd);
 645
 646   // get the ouptut
 647   vtkPolyData *output = vtkPolyData::SafeDownCast(
 648     outInfo1->Get(vtkDataObject::DATA_OBJECT()));
 649
 650   vtkPoints *newPts = vtkPoints::New();
 651   size_t npts = len / 2;
 652   //npts = 10; // DEBUG !
 653   for(size_t i = 0; i < npts; ++i )
 654     {
 655     float x[3];
 656     x[0] = (float)p[i] / 8800;
 657     //std::cout << p[i] << std::endl;
 658     x[1] = i;
 659     x[2] = 0;
 660     vtkIdType ptId = newPts->InsertNextPoint( x );
 661     }
 662   output->SetPoints(newPts);
 663   newPts->Delete();
 664
 665   vtkCellArray* lines = vtkCellArray::New();
 666   for ( int i = 0; i < newPts->GetNumberOfPoints() - 1; ++i )
 667     {
 668     vtkIdType topol[2];
 669     topol[0] = i;
 670     topol[1] = i+1;
 671     lines->InsertNextCell( 2, topol );
 672     }
 673
 674   output->SetLines(lines);
 675   lines->Delete();
 676   output->BuildCells();
 677   //output->GetCellData()->SetScalars(scalars);
 678   //scalars->Delete();
 679
 680   return 1;
 681 }
 682
 683 //----------------------------------------------------------------------------
 684 int vtkGDCMPolyDataReader::RequestData(
 685   vtkInformation *vtkNotUsed(request),
 686   vtkInformationVector **vtkNotUsed(inputVector),
 687   vtkInformationVector *outputVector)
 688 {
 689   vtkInformation *outInfo = outputVector->GetInformationObject(0);
 690   //vtkPoints *newPts, *mergedPts;
 691   //vtkCellArray *newPolys, *mergedPolys;
 692   //vtkFloatArray *newScalars=0, *mergedScalars=0;
 693
 694   // All of the data in the first piece.
 695   if (outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER()) > 0)
 696     {
 697     return 0;
 698     }
 699
 700   if ( !this->FileName || !*this->FileName )
 701     {
 702     vtkErrorMacro(<<"A FileName must be specified.");
 703     return 0;
 704     }
 705
 706   gdcm::Reader reader;
 707   reader.SetFileName( this->FileName );
 708   if( !reader.Read() )
 709     {
 710     return 0;
 711     }
 712
 713   gdcm::MediaStorage ms;
 714   ms.SetFromFile( reader.GetFile() );
 715
 716   int ret;
 717   if( ms == gdcm::MediaStorage::RTStructureSetStorage )
 718     {
 719     ret = this->RequestData_RTStructureSetStorage(reader, outputVector);
 720     }
 721   else if( ms == gdcm::MediaStorage::HemodynamicWaveformStorage)
 722     {
 723     ret = this->RequestData_HemodynamicWaveformStorage(reader, outputVector);
 724     }
 725   else
 726     {
 727     // not handled assume error
 728     ret = 0;
 729     }
 730
 731   return ret;
 732 }
 733
 734 int vtkGDCMPolyDataReader::RequestInformation_RTStructureSetStorage(gdcm::Reader const & reader)
 735 {
 736   const gdcm::DataSet& ds = reader.GetFile().GetDataSet();
 737   // (3006,0020) SQ (Sequence with explicit length #=4)      # 370, 1 StructureSetROISequence
 738   gdcm::Tag tssroisq(0x3006,0x0020);
 739   if( !ds.FindDataElement( tssroisq ) )
 740     {
 741     return 0;
 742     }
 743
 744   const gdcm::DataElement &ssroisq = ds.GetDataElement( tssroisq );
 745   //const gdcm::SequenceOfItems *sqi = ssroisq.GetSequenceOfItems();
 746   gdcm::SmartPointer<gdcm::SequenceOfItems> sqi = ssroisq.GetValueAsSQ();
 747   if( !sqi || !sqi->GetNumberOfItems() )
 748     {
 749     return 0;
 750     }
 751   unsigned int npds = sqi->GetNumberOfItems();
 752
 753   //std::cout << "Nb pd:" << npds << std::endl;
 754   this->SetNumberOfOutputPorts( npds );
 755
 756   // Allocate
 757   for(unsigned int i = 1; i < npds; ++i) // first output is allocated for us
 758     {
 759     vtkPolyData *output2 = vtkPolyData::New();
 760     this->GetExecutive()->SetOutputData(i, output2);
 761     output2->Delete();
 762     }
 763   return 1;
 764 }
 765
 766 int vtkGDCMPolyDataReader::RequestInformation_HemodynamicWaveformStorage(gdcm::Reader const & reader)
 767 {
 768   return 1;
 769 }
 770
 771 //----------------------------------------------------------------------------
 772 int vtkGDCMPolyDataReader::RequestInformation(
 773   vtkInformation *vtkNotUsed(request),
 774   vtkInformationVector **vtkNotUsed(inputVector),
 775   vtkInformationVector *outputVector)
 776 {
 777   // get the info object
 778 //  vtkInformation *outInfo = outputVector->GetInformationObject(0);
 779 //
 780 //  outInfo->Set(vtkStreamingDemandDrivenPipeline::MAXIMUM_NUMBER_OF_PIECES(),
 781 //               -1);
 782   gdcm::Reader reader;
 783   reader.SetFileName( this->FileName );
 784   if( !reader.Read() )
 785     {
 786     return 0;
 787     }
 788
 789   gdcm::MediaStorage ms;
 790   ms.SetFromFile( reader.GetFile() );
 791
 792   int ret;
 793   if( ms == gdcm::MediaStorage::RTStructureSetStorage )
 794     {
 795     ret = this->RequestInformation_RTStructureSetStorage(reader);
 796     }
 797   else if( ms == gdcm::MediaStorage::HemodynamicWaveformStorage)
 798     {
 799     ret = this->RequestInformation_HemodynamicWaveformStorage(reader);
 800     }
 801   else
 802     {
 803     // not handled assume error
 804     ret = 0;
 805     }
 806
 807   if( ret )
 808     {
 809     // Ok let's fill in the 'extra' info:
 810     this->FillMedicalImageInformation(reader);
 811     }
 812
 813   return ret;
 814 }
 815
 816
 817 //----------------------------------------------------------------------------
 818 void vtkGDCMPolyDataReader::PrintSelf(ostream& os, vtkIndent indent)
 819 {
 820   this->Superclass::PrintSelf(os,indent);
 821
 822   os << indent << "File Name: "
 823      << (this->FileName ? this->FileName : "(none)") << "\n";
 824
 825
 826 }