common/vtkGDCMPolyDataReader.cxx

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