]> Creatis software - clitk.git/blob - common/vtkGDCMPolyDataReader.cxx
obsolete since the user manually update his repo
[clitk.git] / 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 }