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 /*
15  * tar is a UNIX program for archiving.
16  * Two types of operations are possible: concatenate / extract
17  * Thus the name of 'gdcmtar' to concatenate a list of 2D slices into a multi frames
18  * and the other way around: extract 2D slices from a multi frames image
19  * It also support the fake multi frame image (CSA MOSAIC)
20  */
21 
22 #include "gdcmReader.h"
23 #include "gdcmCSAHeader.h"
24 #include "gdcmVersion.h"
25 #include "gdcmImageReader.h"
26 #include "gdcmDataElement.h"
27 #include "gdcmImageWriter.h"
28 #include "gdcmSplitMosaicFilter.h"
29 #include "gdcmFilename.h"
30 #include "gdcmFilenameGenerator.h"
31 #include "gdcmDirectionCosines.h"
32 #include "gdcmImageHelper.h"
33 #include "gdcmUIDGenerator.h"
34 #include "gdcmUIDs.h"
35 #include "gdcmGlobal.h"
36 #include "gdcmDirectory.h"
37 #include "gdcmScanner.h"
38 #include "gdcmIPPSorter.h"
39 #include "gdcmAttribute.h"
40 #include "gdcmAnonymizer.h"
41 #include "gdcmTagKeywords.h"
42 #include "gdcmMrProtocol.h"
43 
44 #include <string>
45 #include <iostream>
46 
47 #include <stdio.h>     /* for printf */
48 #include <stdlib.h>    /* for exit */
49 #include <getopt.h>
50 #include <string.h>
51 
PrintVersion()52 static void PrintVersion()
53 {
54   std::cout << "gdcmtar: gdcm " << gdcm::Version::GetVersion() << " ";
55   const char date[] = "$Date$";
56   std::cout << date << std::endl;
57 }
58 
PrintHelp()59 static void PrintHelp()
60 {
61   PrintVersion();
62   std::cout << "Usage: gdcmtar [OPTION] [FILE]" << std::endl;
63   std::cout << "Concatenate/Extract DICOM files.\n";
64   std::cout << "Parameter (required):" << std::endl;
65   std::cout << "  -i --input     DICOM filename" << std::endl;
66   std::cout << "  -o --output    DICOM filename" << std::endl;
67   std::cout << "Options:" << std::endl;
68   std::cout << "     --enhance        Enhance (default)" << std::endl;
69   std::cout << "  -U --unenhance      Unenhance" << std::endl;
70   std::cout << "  -M --mosaic         Split SIEMENS Mosaic image into multiple frames." << std::endl;
71   std::cout << "     --mosaic-private When splitting SIEMENS Mosaic image into multiple frames, preserve private attributes (advanced user only)." << std::endl;
72   std::cout << "  -p --pattern        Specify trailing file pattern." << std::endl;
73   std::cout << "     --root-uid       Root UID." << std::endl;
74   //std::cout << "     --resources-path     Resources path." << std::endl;
75   std::cout << "General Options:" << std::endl;
76   std::cout << "  -V --verbose    more verbose (warning+error)." << std::endl;
77   std::cout << "  -W --warning    print warning info." << std::endl;
78   std::cout << "  -D --debug      print debug info." << std::endl;
79   std::cout << "  -E --error      print error info." << std::endl;
80   std::cout << "  -h --help       print help." << std::endl;
81   std::cout << "  -v --version    print version." << std::endl;
82   std::cout << "Env var:" << std::endl;
83   std::cout << "  GDCM_ROOT_UID Root UID" << std::endl;
84   //std::cout << "  GDCM_RESOURCES_PATH path pointing to resources files (Part3.xml, ...)" << std::endl;
85 }
86 
87 /*
88  * The following example is a basic sorted which should work in generic cases.
89  * It sort files based on:
90  * Study Instance UID
91  *   Series Instance UID
92  *     Frame of Reference UID
93  *       Image Orientation (Patient)
94  *         Image Position (Patient) (Sorting based on IPP + IOP)
95  */
96 
97 namespace gdcm {
98   const Tag T0(0x0008,0x0016); // SOP Class UID
99   const Tag T1(0x0020,0x000d); // Study Instance UID
100   const Tag T2(0x0020,0x000e); // Series Instance UID
101   const Tag T3(0x0020,0x0052); // Frame of Reference UID
102   const Tag T4(0x0020,0x0037); // Image Orientation (Patient)
103   const Tag T5(0x0008,0x0008); // Image Type
104 
105 class DiscriminateVolume
106 {
107 private:
108   static const bool debuggdcmtar = false;
109   std::vector< Directory::FilenamesType > SortedFiles;
110   std::vector< Directory::FilenamesType > UnsortedFiles;
111 
GetAllFilenamesFromTagToValue(Scanner const & s,Directory::FilenamesType const & filesubset,Tag const & t,const char * valueref)112   Directory::FilenamesType GetAllFilenamesFromTagToValue(
113     Scanner const & s, Directory::FilenamesType const &filesubset, Tag const &t, const char *valueref)
114 {
115   Directory::FilenamesType theReturn;
116   if( valueref )
117     {
118     size_t len = strlen( valueref );
119     Directory::FilenamesType::const_iterator file = filesubset.begin();
120     for(; file != filesubset.end(); ++file)
121       {
122       const char *filename = file->c_str();
123       const char * value = s.GetValue(filename, t);
124       if( value && strncmp(value, valueref, len ) == 0 )
125         {
126         theReturn.push_back( filename );
127         }
128       }
129     }
130   return theReturn;
131 }
132 
133 class frame_diff /* */
134 {
135 public:
operator <(const frame_diff & rhs) const136   bool operator<(const frame_diff &rhs) const
137     {
138     (void)rhs;
139     return true;
140     }
141 };
142 
get_frame_diff(const char * filename)143 static frame_diff get_frame_diff( const char* filename )
144 {
145   frame_diff ret;
146   gdcm::Reader reader;
147   reader.SetFileName( filename );
148 
149   gdcm::Tag dummy(0x7fe0,0x0000);
150   std::set<gdcm::Tag> skiptags;
151   skiptags.insert( dummy );
152   if ( !reader.ReadUpToTag( dummy, skiptags) )
153     {
154     assert(0);
155     }
156 
157   gdcm::CSAHeader csa;
158   const gdcm::DataSet& ds = reader.GetFile().GetDataSet();
159 
160   const gdcm::PrivateTag &t1 = csa.GetCSAImageHeaderInfoTag();
161   //const gdcm::PrivateTag &t2 = csa.GetCSASeriesHeaderInfoTag();
162   //const gdcm::PrivateTag &t3 = csa.GetCSADataInfo();
163 
164   //bool found = false;
165   if( ds.FindDataElement( t1 ) )
166     {
167     csa.LoadFromDataElement( ds.GetDataElement( t1 ) );
168     assert( csa.GetFormat() == gdcm::CSAHeader::SV10
169       || csa.GetFormat() == gdcm::CSAHeader::NOMAGIC );
170     // SIEMENS / Diffusion
171     const CSAElement &bvalue  = csa.GetCSAElementByName( "B_value" );
172     (void)bvalue;
173     const CSAElement &bmatrix = csa.GetCSAElementByName( "B_matrix" );
174     (void)bmatrix;
175     const CSAElement &dgd = csa.GetCSAElementByName( "DiffusionGradientDirection" );
176     (void)dgd;
177     //assert(0);
178     }
179   else
180     {
181     assert(0);
182     }
183 
184   return ret;
185 }
186 
ProcessAIOP(Scanner const & s,Directory::FilenamesType const & subset,const char * iopval)187 void ProcessAIOP(Scanner const & s, Directory::FilenamesType const & subset, const char *iopval)
188 {
189   if( debuggdcmtar )
190     std::cout << "IOP: " << iopval << std::endl;
191 
192   bool is_diffusion = false;
193   gdcm::Scanner::ValuesType imagetypes = s.GetValues(gdcm::T5);
194   if( imagetypes.size() == 1 )
195     {
196     const std::string & imagetype = *imagetypes.begin();
197     //gdcm::Attribute<0x0008, 0x0008> at;
198     //at.SetValue( imagetype );
199     std::string::size_type pos = imagetype.find( "DIFFUSION" );
200     is_diffusion = (pos != std::string::npos);
201     }
202 
203   if( is_diffusion )
204     {
205     std::multimap<frame_diff, std::string> mm;
206     for( Directory::FilenamesType::const_iterator file = subset.begin();
207       file != subset.end(); ++file)
208       {
209       //std::cerr << *file << std::endl;
210       mm.insert( std::make_pair(get_frame_diff( file->c_str() ), file->c_str()) );
211       }
212     //assert(0);
213     for( std::multimap<frame_diff, std::string>::const_iterator it = mm.begin();
214       it != mm.end(); ++it )
215       {
216       std::cerr << it->second << std::endl;
217       }
218       SortedFiles.push_back( subset );
219     }
220   else
221     {
222     IPPSorter ipp;
223     ipp.SetComputeZSpacing( true );
224     ipp.SetZSpacingTolerance( 1e-3 ); // ??
225     bool b = ipp.Sort( subset );
226     if( !b )
227       {
228       // If you reach here this means you need one more parameter to discriminate this
229       // series. Eg. T1 / T2 intertwined. Multiple Echo (0018,0081)
230       if( debuggdcmtar )
231         {
232         std::cerr << "Failed to sort: " << subset.begin()->c_str() << std::endl;
233         for(
234           Directory::FilenamesType::const_iterator file = subset.begin();
235           file != subset.end(); ++file)
236           {
237           std::cerr << *file << std::endl;
238           }
239         }
240       UnsortedFiles.push_back( subset );
241       return ;
242       }
243     if( debuggdcmtar )
244       ipp.Print( std::cout );
245     SortedFiles.push_back( ipp.GetFilenames() );
246     }
247 }
248 
ProcessAFrameOfRef(Scanner const & s,Directory::FilenamesType const & subset,const char * frameuid)249 void ProcessAFrameOfRef(Scanner const & s, Directory::FilenamesType const & subset, const char * frameuid)
250 {
251   // In this subset of files (belonging to same series), let's find those
252   // belonging to the same Frame ref UID:
253   Directory::FilenamesType files = GetAllFilenamesFromTagToValue(
254     s, subset, T3, frameuid);
255 
256   std::set< std::string > iopset;
257 
258   for(
259     Directory::FilenamesType::const_iterator file = files.begin();
260     file != files.end(); ++file)
261     {
262     //std::cout << *file << std::endl;
263     const char * value = s.GetValue(file->c_str(), gdcm::T4 );
264     assert( value );
265     iopset.insert( value );
266     }
267   size_t n = iopset.size();
268   if ( n == 0 )
269     {
270     assert( files.empty() );
271     return;
272     }
273 
274   if( debuggdcmtar )
275     std::cout << "Frame of Ref: " << frameuid << std::endl;
276   if ( n == 1 )
277     {
278     ProcessAIOP(s, files, iopset.begin()->c_str() );
279     }
280   else
281     {
282     const char *f = files.begin()->c_str();
283     if( debuggdcmtar )
284       std::cerr << "More than one IOP: " << f << std::endl;
285     // Make sure that there is actually 'n' different IOP
286     gdcm::DirectionCosines ref;
287     gdcm::DirectionCosines dc;
288     for(
289       std::set< std::string >::const_iterator it = iopset.begin();
290       it != iopset.end(); ++it )
291       {
292       ref.SetFromString( it->c_str() );
293       for(
294         Directory::FilenamesType::const_iterator file = files.begin();
295         file != files.end(); ++file)
296         {
297         std::string value = s.GetValue(file->c_str(), gdcm::T4 );
298         if( value != it->c_str() )
299           {
300           dc.SetFromString( value.c_str() );
301           const double crossdot = ref.CrossDot(dc);
302           const double eps = std::fabs( 1. - crossdot );
303           if( eps < 1e-6 )
304             {
305             std::cerr << "Problem with IOP discrimination: " << file->c_str()
306               << " " << it->c_str() << std::endl;
307             return;
308             }
309           }
310         }
311       }
312       // If we reach here this means there is actually 'n' different IOP
313     for(
314       std::set< std::string >::const_iterator it = iopset.begin();
315       it != iopset.end(); ++it )
316       {
317       const char *iopvalue = it->c_str();
318       Directory::FilenamesType iopfiles = GetAllFilenamesFromTagToValue(
319         s, files, T4, iopvalue );
320       ProcessAIOP(s, iopfiles, iopvalue );
321       }
322     }
323 }
324 
ProcessASeries(Scanner const & s,const char * seriesuid)325 void ProcessASeries(Scanner const & s, const char * seriesuid)
326 {
327   if( debuggdcmtar )
328     std::cout << "Series: " << seriesuid << std::endl;
329   // let's find all files belonging to this series:
330   Directory::FilenamesType seriesfiles = GetAllFilenamesFromTagToValue(
331     s, s.GetFilenames(), T2, seriesuid);
332 
333   gdcm::Scanner::ValuesType vt3 = s.GetValues(T3);
334   for(
335     gdcm::Scanner::ValuesType::const_iterator it = vt3.begin()
336     ; it != vt3.end(); ++it )
337     {
338     ProcessAFrameOfRef(s, seriesfiles, it->c_str());
339     }
340 }
341 
ProcessAStudy(Scanner const & s,const char * studyuid)342 void ProcessAStudy(Scanner const & s, const char * studyuid)
343 {
344   if( debuggdcmtar )
345     std::cout << "Study: " << studyuid << std::endl;
346   gdcm::Scanner::ValuesType vt2 = s.GetValues(T2);
347   if( vt2.empty() )
348     std::cerr << "No Series Found" << std::endl;
349   for(
350     gdcm::Scanner::ValuesType::const_iterator it = vt2.begin()
351     ; it != vt2.end(); ++it )
352     {
353     ProcessASeries(s, it->c_str());
354     }
355 }
356 public:
357 
Print(std::ostream & os)358 void Print( std::ostream & os )
359 {
360   os << "Sorted Files: " << std::endl;
361   for(
362     std::vector< Directory::FilenamesType >::const_iterator it = SortedFiles.begin();
363     it != SortedFiles.end(); ++it )
364     {
365     os << "Group: " << std::endl;
366     for(
367       Directory::FilenamesType::const_iterator file = it->begin();
368       file != it->end(); ++file)
369       {
370       os << *file << std::endl;
371       }
372     }
373   os << "Unsorted Files: " << std::endl;
374   for(
375     std::vector< Directory::FilenamesType >::const_iterator it = UnsortedFiles.begin();
376     it != UnsortedFiles.end(); ++it )
377     {
378     os << "Group: " << std::endl;
379     for(
380       Directory::FilenamesType::const_iterator file = it->begin();
381       file != it->end(); ++file)
382       {
383       os << *file << std::endl;
384       }
385     }
386 }
387 
GetSortedFiles() const388   std::vector< Directory::FilenamesType > const & GetSortedFiles() const { return SortedFiles; }
GetUnsortedFiles() const389   std::vector< Directory::FilenamesType > const & GetUnsortedFiles() const { return UnsortedFiles; }
390 
ProcessIntoVolume(Scanner const & s)391 void ProcessIntoVolume( Scanner const & s )
392 {
393   gdcm::Scanner::ValuesType vt1 = s.GetValues( gdcm::T1 );
394   for(
395     gdcm::Scanner::ValuesType::const_iterator it = vt1.begin()
396     ; it != vt1.end(); ++it )
397     {
398     ProcessAStudy( s, it->c_str() );
399     }
400 }
401 
402 };
403 
ConcatenateImages(Image & im1,Image const & im2)404 static bool ConcatenateImages(Image &im1, Image const &im2)
405 {
406   DataElement& de1 = im1.GetDataElement();
407   if( de1.GetByteValue() )
408     {
409     const ByteValue *bv1 = de1.GetByteValue();
410     std::vector<char> v1 = *bv1;
411     const DataElement& de2 = im2.GetDataElement();
412     const ByteValue *bv2 = de2.GetByteValue();
413     const std::vector<char> & v2 = *bv2;
414     v1.insert( v1.end(), v2.begin(), v2.end() );
415 
416     de1.SetByteValue(&v1[0], (uint32_t)v1.size());
417     }
418   else if( de1.GetSequenceOfFragments() )
419     {
420     SequenceOfFragments *sqf1 = de1.GetSequenceOfFragments();
421     assert( sqf1 );
422     const DataElement& de2 = im2.GetDataElement();
423     const SequenceOfFragments *sqf2 = de2.GetSequenceOfFragments();
424     assert( sqf2 );
425     assert( sqf2->GetNumberOfFragments() == 1 );
426     const Fragment& frag = sqf2->GetFragment(0);
427     sqf1->AddFragment(frag);
428     }
429   else
430     {
431     return false;
432     }
433 
434   // update meta info
435   unsigned int z = im1.GetDimension(2);
436   im1.SetDimension(2, z + 1 );
437   return true;
438 }
439 
InsertIfExist(DataSet & enhds,const DataSet & ds,const Tag & t)440 static void InsertIfExist( DataSet &enhds, const DataSet &ds, const Tag & t )
441 {
442   if( ds.FindDataElement( t ) )
443     {
444     assert( !enhds.FindDataElement( t ) );
445     enhds.Insert( ds.GetDataElement( t ) );
446     }
447 }
448 
449 
450 } // namespace gdcm
451 
452 
MakeImageEnhanced(std::string const & filename,std::string const & outfilename)453 static int MakeImageEnhanced( std::string const & filename, std::string const &outfilename )
454 {
455   if( !gdcm::System::FileIsDirectory(filename.c_str()) )
456     {
457     std::cerr << "Input needs to be directory" << std::endl;
458     return 1;
459     }
460 
461   gdcm::Directory d;
462   d.Load( filename.c_str(), true ); // recursive !
463 
464   gdcm::Scanner s;
465   s.AddTag( gdcm::T0 );
466   s.AddTag( gdcm::T1 );
467   s.AddTag( gdcm::T2 );
468   s.AddTag( gdcm::T3 );
469   s.AddTag( gdcm::T4 );
470   s.AddTag( gdcm::T5 );
471   bool b = s.Scan( d.GetFilenames() );
472   if( !b )
473     {
474     std::cerr << "Scanner failed" << std::endl;
475     return 1;
476     }
477 
478   // For now only accept MR Image Storage
479   gdcm::Scanner::ValuesType vt = s.GetValues(gdcm::T0);
480   if( vt.size() != 1 ) return 1;
481 
482   const char *sop = vt.begin()->c_str();
483   gdcm::MediaStorage msorig = gdcm::MediaStorage::GetMSType( sop );
484   if( msorig != gdcm::MediaStorage::MRImageStorage
485    && msorig != gdcm::MediaStorage::CTImageStorage )
486     {
487     std::cerr << "Sorry MediaStorage not supported: [" << sop << "]" << std::endl;
488     return 1;
489     }
490 
491   gdcm::DiscriminateVolume dv;
492   dv.ProcessIntoVolume( s );
493 //  dv.Print( std::cout );
494 
495   // gdcm::DataElement &de = im.GetImage().GetDataElement();
496   std::vector< gdcm::Directory::FilenamesType > const &sorted = dv.GetSortedFiles();
497   if( !gdcm::System::MakeDirectory( outfilename.c_str() ) )
498     {
499     std::cerr << "Could not create dir: " << outfilename << std::endl;
500     return 1;
501     }
502   for(
503     std::vector< gdcm::Directory::FilenamesType >::const_iterator it = sorted.begin();
504     it != sorted.end(); ++it )
505     {
506     gdcm::ImageWriter im0;
507 
508     gdcm::Directory::FilenamesType const & files = *it;
509     gdcm::Directory::FilenamesType::const_iterator file = files.begin();
510 
511     const char *reffile = file->c_str();
512     // construct the target dir:
513     const char* studyuid = s.GetValue(reffile, gdcm::T1);
514     const char* seriesuid = s.GetValue(reffile, gdcm::T2);
515     const char* frameuid = s.GetValue(reffile, gdcm::T3);
516     std::string targetdir = outfilename;
517     targetdir += '/';
518     targetdir += studyuid;
519     targetdir += '/';
520     targetdir += seriesuid;
521     targetdir += '/';
522     targetdir += frameuid;
523     // construct the target name:
524     std::string targetname = targetdir;
525 
526     targetdir += "/old";
527 
528     // make sure the dir exist first:
529     if( !gdcm::System::MakeDirectory( targetdir.c_str() ) )
530       {
531       std::cerr << "Could not create dir: " << targetdir << std::endl;
532       return 1;
533       }
534 
535     gdcm::FilenameGenerator fg;
536     fg.SetNumberOfFilenames( files.size() );
537     fg.SetPrefix( targetdir.c_str() );
538     fg.SetPattern( "%04d.dcm" );
539     if( !fg.Generate() )
540       {
541       assert( 0 );
542       }
543 
544     gdcm::ImageReader reader0;
545     reader0.SetFileName( reffile );
546     if( !reader0.Read() )
547       {
548       assert( 0 );
549       }
550     gdcm::Image &currentim = reader0.GetImage();
551     assert( currentim.GetNumberOfDimensions( ) == 2 );
552     currentim.SetNumberOfDimensions( 3 );
553     size_t count = 0;
554 
555     //gdcm::ImageWriter writer;
556     gdcm::Writer writer0;
557     writer0.SetFileName( fg.GetFilename( count ) );
558     writer0.SetFile( reader0.GetFile() );
559     writer0.GetFile().GetHeader().Clear();
560     if( !writer0.Write() )
561       {
562       assert( 0 );
563       }
564     ++file;
565     ++count;
566 
567     // insert Enhanced MR Image Storage
568     if( false )
569       {
570       gdcm::DataSet & enhds = im0.GetFile().GetDataSet();
571       gdcm::DataSet & ds = reader0.GetFile().GetDataSet();
572       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0005) ); // SpecificCharacterSet
573       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0008) ); // ImageType
574       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0012) ); // InstanceCreationDate
575       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0013) ); // InstanceCreationTime
576       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0020) ); // StudyDate
577       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0021) ); // StudyTime
578       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0022) ); // AcquisitionDate
579       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0023) ); // ContentDate
580       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0030) ); // StudyTime
581       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0031) ); // SeriesTime
582       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0032) ); // AcquisitionTime
583       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0033) ); // ContentTime
584       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0050) ); // AccessionNumber
585       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0070) ); // Manufacturer
586       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0080) ); // InstitutionName
587       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0081) ); // InstitutionAddress
588       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x0090) ); // ReferringPhysicianName
589       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x1010) ); // StationName
590       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x1030) ); // StudyDescription
591       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x103e) ); // SeriesDescription
592       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x1040) ); // InstitutionalDepartme
593       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x1050) ); // PerformingPhysicianNa
594       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x8,0x1090) ); // ManufacturerModelName
595 
596       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0010,0x0010)); // PatientName
597       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0010,0x0020)); //
598       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0010,0x0030)); //
599       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0010,0x0040)); //
600       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0010,0x1010)); //
601       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0010,0x1030)); //
602       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0020)); //
603       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0021)); //
604       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0022)); //
605       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0023)); //
606       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0024)); //
607       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0025)); //
608       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0050)); //
609       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0080)); //
610       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0081)); //
611       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0083)); //
612       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0084)); //
613       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0085)); //
614       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0086)); //
615       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0087)); //
616       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0088)); //
617       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0089)); //
618       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0091)); //
619       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0093)); //
620       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0094)); //
621       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x0095)); //
622       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x1000)); //
623       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x1020)); //
624       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x1030)); //
625       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x1251)); //
626       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x1310)); //
627       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x1312)); //
628       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x1314)); //
629       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x1315)); //
630       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x1316)); //
631       //gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x1318)); //
632       gdcm::InsertIfExist( enhds, ds, gdcm::Tag(0x0018,0x5100)); //
633 
634         {
635         gdcm::SmartPointer<gdcm::SequenceOfItems> sqi1 = new gdcm::SequenceOfItems;
636 
637 /*
638 (5200,9229) SQ (Sequence with undefined length)                   # u/l,1 Shared Functional Groups Sequence
639   (fffe,e000) na (Item with undefined length)
640     (0018,9006) SQ (Sequence with undefined length)               # u/l,1 MR Imaging Modifier Sequence
641       (fffe,e000) na (Item with undefined length)
642         (0018,0095) DS [109 ]                                     # 4,1 Pixel Bandwidth
643         (0018,9020) CS [NONE]                                     # 4,1 Magnetization Transfer
644         (0018,9022) CS [NO]                                       # 2,1 Blood Signal Nulling
645         (0018,9028) CS [NONE]                                     # 4,1 Tagging
646         (0018,9098) FD 63.8799                                    # 8,1-2 Transmitter Frequency
647 */
648         gdcm::Item item0;
649         gdcm::DataSet &subds0 = item0.GetNestedDataSet();
650 
651           {
652         gdcm::SmartPointer<gdcm::SequenceOfItems> sqi0 = new gdcm::SequenceOfItems;
653           const gdcm::Tag sisq0(0x0018,0x9006);
654           gdcm::DataElement de0( sisq0 );
655           de0.SetVR( gdcm::VR::SQ );
656           de0.SetValue( *sqi0 );
657           subds0.Insert( de0 );
658 
659           gdcm::Item item1;
660           gdcm::DataSet &subds1 = item1.GetNestedDataSet();
661           gdcm::InsertIfExist( subds1, ds, gdcm::Tag(0x0018,0x0095)); //
662           gdcm::InsertIfExist( subds1, ds, gdcm::Tag(0x0018,0x9020)); //
663           gdcm::InsertIfExist( subds1, ds, gdcm::Tag(0x0018,0x9022)); //
664           gdcm::InsertIfExist( subds1, ds, gdcm::Tag(0x0018,0x9028)); //
665           gdcm::InsertIfExist( subds1, ds, gdcm::Tag(0x0018,0x9098)); //
666           sqi0->AddItem( item1 );
667           }
668 
669 /*
670     (0018,9042) SQ (Sequence with undefined length)               # u/l,1 MR Receive Coil Sequence
671       (fffe,e000) na (Item with undefined length)
672         (0018,1250) SH [MULTI COIL]                               # 10,1 Receive Coil Name
673         (0018,9041) LO (no value)                                 # 0,1 Receive Coil Manufacturer Name
674         (0018,9043) CS [MULTICOIL ]                               # 10,1 Receive Coil Type
675         (0018,9044) CS [NO]                                       # 2,1 Quadrature Receive Coil
676         (0018,9045) SQ (Sequence with undefined length)           # u/l,1 Multi-Coil Definition Sequence
677           (fffe,e000) na (Item with undefined length)
678             (0018,9047) SH [MULTI ELEMENT ]                       # 14,1 Multi-Coil Element Name
679             (0018,9048) CS [YES ]                                 # 4,1 Multi-Coil Element Used
680           (fffe,e00d)
681         (fffe,e0dd)
682       (fffe,e00d)
683 */
684           {
685         gdcm::SmartPointer<gdcm::SequenceOfItems> sqi0 = new gdcm::SequenceOfItems;
686           const gdcm::Tag sisq0(0x0018,0x9042);
687           gdcm::DataElement de0( sisq0 );
688           de0.SetVR( gdcm::VR::SQ );
689           de0.SetValue( *sqi0 );
690           subds0.Insert( de0 );
691 
692           gdcm::Item item1;
693           gdcm::DataSet &subds1 = item1.GetNestedDataSet();
694           gdcm::InsertIfExist( subds1, ds, gdcm::Tag(0x0018,0x1250)); //
695           gdcm::InsertIfExist( subds1, ds, gdcm::Tag(0x0018,0x9041)); //
696           gdcm::InsertIfExist( subds1, ds, gdcm::Tag(0x0018,0x9043)); //
697           gdcm::InsertIfExist( subds1, ds, gdcm::Tag(0x0018,0x9044)); //
698           //gdcm::InsertIfExist( subds1, ds, gdcm::Tag(0x0018,0x9045)); //
699           sqi0->AddItem( item1 );
700           }
701 
702         sqi1->AddItem( item0 );
703         const gdcm::Tag sisq1(0x5200,0x9229);
704         gdcm::DataElement de1( sisq1 );
705         de1.SetVR( gdcm::VR::SQ );
706         de1.SetValue( *sqi1 );
707         enhds.Insert( de1 );
708         }
709       }
710 
711     for( ; file != files.end(); ++file, ++count )
712       {
713       gdcm::ImageReader reader;
714       reader.SetFileName( file->c_str() );
715       if( !reader.Read() )
716         {
717         assert( 0 );
718         }
719       const gdcm::Image &im = reader.GetImage();
720 
721       //gdcm::ImageWriter writer;
722       gdcm::Writer writer;
723       writer.SetFileName( fg.GetFilename( count ) );
724       writer.SetFile( reader.GetFile() );
725       writer.GetFile().GetHeader().Clear();
726       if( !writer.Write() )
727         {
728         assert( 0 );
729         }
730 
731       if( !ConcatenateImages(currentim, im) )
732         {
733         assert( 0 );
734         }
735 
736       // check consistency
737       if( false )
738         {
739         //gdcm::DataSet & enhds = im0.GetFile().GetDataSet();
740         //gdcm::DataSet & ds = reader.GetFile().GetDataSet();
741         // TODO
742         }
743       }
744 
745     im0.SetFileName( (targetname + "/new.dcm").c_str() );
746     //  im.SetFile( reader.GetFile() );
747 
748     gdcm::DataSet &ds = im0.GetFile().GetDataSet();
749 
750     gdcm::MediaStorage ms;
751     switch( msorig )
752       {
753     case gdcm::MediaStorage::CTImageStorage:
754       ms = gdcm::MediaStorage::EnhancedCTImageStorage;
755       break;
756     case gdcm::MediaStorage::MRImageStorage:
757       ms = gdcm::MediaStorage::EnhancedMRImageStorage;
758       break;
759     default:
760       return 1;
761       }
762 
763     gdcm::DataElement de( gdcm::Tag(0x0008, 0x0016) );
764     const char* msstr = gdcm::MediaStorage::GetMSString(ms);
765     de.SetByteValue( msstr, (uint32_t)strlen(msstr) );
766     de.SetVR( gdcm::Attribute<0x0008, 0x0016>::GetVR() );
767     ds.Insert( de );
768 
769     im0.SetImage( currentim );
770     if( !im0.Write() )
771       {
772       std::cerr << "Could not write: " << std::endl;
773       return 1;
774       }
775 
776     }
777 
778   std::vector< gdcm::Directory::FilenamesType > const &unsorted = dv.GetUnsortedFiles();
779   if( !unsorted.empty() )
780     {
781     std::string targetdir3 = outfilename;
782     targetdir3 += "/unhandled";
783     if( !gdcm::System::MakeDirectory( targetdir3.c_str() ) )
784       {
785       std::cerr << "Could not create dir: " << targetdir3 << std::endl;
786       return 1;
787       }
788     std::cerr << "Could not process the following files (please report): " << std::endl;
789     std::vector< gdcm::Directory::FilenamesType >::const_iterator it = unsorted.begin();
790     for( ; it != unsorted.end(); ++it )
791       {
792       gdcm::Directory::FilenamesType const & files = *it;
793       gdcm::Directory::FilenamesType::const_iterator file = files.begin();
794       for( ; file != files.end(); ++file )
795         {
796         const char *f = file->c_str();
797         std::string targetdir2 = outfilename;
798         targetdir2 += "/unhandled/";
799         gdcm::Filename fn2( f );
800         const char *outfn2 = fn2.GetName();
801         targetdir2 += outfn2;
802         //std::cerr << f << " -> " << targetdir2 << std::endl;
803         std::ifstream in( f, std::ios::binary );
804         std::ofstream out( targetdir2.c_str() , std::ios::binary );
805         out << in.rdbuf();
806         }
807       }
808     }
809 
810   return 0;
811 }
812 
813 namespace gdcm
814 {
815 
ReplaceIf(DataSet & rootds,const DataSet & ds,const Tag & t1,const Tag & t2)816 static inline void ReplaceIf(DataSet &rootds, const DataSet &ds, const Tag & t1, const Tag & t2 )
817 {
818   if( !ds.FindDataElement( t1 ) ) return ;
819   SmartPointer<SequenceOfItems> sqi1 = ds.GetDataElement( t1 ).GetValueAsSQ();
820   if( !sqi1 || sqi1->IsEmpty() ) return ;
821   const Item &item1 = sqi1->GetItem(1);
822   const DataSet & ds1 = item1.GetNestedDataSet();
823   if( !ds1.FindDataElement( t2 ) ) return ;
824   rootds.Replace(ds1.GetDataElement( t2 ));
825 }
826 
RemapSharedIntoOld(gdcm::DataSet & ds,SequenceOfItems * sfgs,SequenceOfItems * pffgs,unsigned int index)827 static bool RemapSharedIntoOld( gdcm::DataSet & ds,
828  SequenceOfItems *sfgs,
829  SequenceOfItems *pffgs,
830  unsigned int index )
831 {
832   assert( sfgs );
833   assert( pffgs );
834 
835   assert( sfgs->GetNumberOfItems() == 1 );
836   Item const &item1 = sfgs->GetItem( 1 );
837   const DataSet & sfgs_ds = item1.GetNestedDataSet();
838 #if 1
839   // Repetition Time
840   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9112), Tag(0x0018,0x0080) );
841   // Echo Train Length
842   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9112), Tag(0x0018,0x0091) );
843   // Flip Angle
844   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9112), Tag(0x0018,0x1314) );
845   // Number of Averages
846   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9119), Tag(0x0018,0x0083) );
847 
848   // Percent Sampling
849   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9125), Tag(0x0018,0x0093) );
850   // Percent Phase Field of View
851   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9125), Tag(0x0018,0x0094) );
852   // Receive Coil Name
853   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9042), Tag(0x0018,0x1250) );
854   // Transmit Coil Name
855   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9049), Tag(0x0018,0x1251) );
856   // InPlanePhaseEncodingDirection
857   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9125), Tag(0x0018,0x1312) );
858   // TransmitterFrequency
859   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9006), Tag(0x0018,0x9098) );
860   // InversionRecovery
861   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9115), Tag(0x0018,0x9009) );
862   // FlowCompensation
863   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9115), Tag(0x0018,0x9010) );
864   // ReceiveCoilType
865   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9042), Tag(0x0018,0x9043) );
866   // QuadratureReceiveCoil
867   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9042), Tag(0x0018,0x9044) );
868   // SlabThickness
869   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9107), Tag(0x0018,0x9104) );
870   // MultiCoilDefinitionSequence
871   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9042), Tag(0x0018,0x9045) );
872   // SlabOrientation
873   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9107), Tag(0x0018,0x9105) );
874   // MidSlabPosition
875   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9107), Tag(0x0018,0x9106) );
876   // OperatingModeSequence
877   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9112), Tag(0x0018,0x9176) );
878   // MRAcquisitionPhaseEncodingStepsOutOf
879   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9125), Tag(0x0018,0x9232) );
880   // SpecificAbsorptionRateSequence
881   ReplaceIf(ds, sfgs_ds, Tag(0x0018,0x9112), Tag(0x0018,0x9239) );
882   // AnatomicRegionSequence
883   ReplaceIf(ds, sfgs_ds, Tag(0x0020,0x9071), Tag(0x0008,0x2218) );
884   // Purpose of Reference Code Sequence
885   // FIXME what if there is multiple purpose of rcs ?
886   ReplaceIf(ds, sfgs_ds, Tag(0x0008,0x1140), Tag(0x0040,0xa170) );
887 #else
888   for(
889     DataSet::ConstIterator it = sfgs_ds.Begin();
890     it != sfgs_ds.End(); ++it )
891     {
892     ds.Replace( *it );
893     }
894 #endif
895 
896   Item const &item2 = pffgs->GetItem( index + 1 );
897   const DataSet & pffgs_ds = item2.GetNestedDataSet();
898 
899 #if 1
900   // Effective Echo Time
901   ReplaceIf(ds, pffgs_ds, Tag(0x0018,0x9114), Tag(0x0018,0x9082) );
902   // -> should also be Echo Time
903   // Nominal Cardiac Trigger Delay Time
904   ReplaceIf(ds, pffgs_ds, Tag(0x0018,0x9118), Tag(0x0020,0x9153) );
905   // Metabolite Map Description
906   ReplaceIf(ds, pffgs_ds, Tag(0x0018,0x9152), Tag(0x0018,0x9080) );
907   // IPP
908   ReplaceIf(ds, pffgs_ds, Tag(0x0020,0x9113), Tag(0x0020,0x0032) );
909   // IOP
910   ReplaceIf(ds, pffgs_ds, Tag(0x0020,0x9116), Tag(0x0020,0x0037) );
911   // Slice Thickness
912   ReplaceIf(ds, pffgs_ds, Tag(0x0028,0x9110), Tag(0x0018,0x0050) );
913   // Pixel Spacing
914   ReplaceIf(ds, pffgs_ds, Tag(0x0028,0x9110), Tag(0x0028,0x0030) );
915 
916   // window level
917   ReplaceIf(ds, pffgs_ds, Tag(0x0028,0x9132), Tag(0x0028,0x1050) );
918   ReplaceIf(ds, pffgs_ds, Tag(0x0028,0x9132), Tag(0x0028,0x1051) );
919 
920   // rescale slope/intercept
921   ReplaceIf(ds, pffgs_ds, Tag(0x0028,0x9145), Tag(0x0028,0x1052) );
922   ReplaceIf(ds, pffgs_ds, Tag(0x0028,0x9145), Tag(0x0028,0x1053) );
923   ReplaceIf(ds, pffgs_ds, Tag(0x0028,0x9145), Tag(0x0028,0x1054) );
924 
925   // FrameReferenceDateTime
926   ReplaceIf(ds, pffgs_ds, Tag(0x0020,0x9111), Tag(0x0018,0x9151) );
927   // FrameAcquisitionDuration
928   ReplaceIf(ds, pffgs_ds, Tag(0x0020,0x9111), Tag(0x0018,0x9220) );
929   // TemporalPositionIndex
930   ReplaceIf(ds, pffgs_ds, Tag(0x0020,0x9111), Tag(0x0020,0x9128) );
931   // InStackPositionNumber
932   ReplaceIf(ds, pffgs_ds, Tag(0x0020,0x9111), Tag(0x0020,0x9057) );
933   // FrameType
934   ReplaceIf(ds, pffgs_ds, Tag(0x0018,0x9226), Tag(0x0008,0x9007) );
935   // DimensionIndexValues
936   ReplaceIf(ds, pffgs_ds, Tag(0x0020,0x9111), Tag(0x0020,0x9157) );
937   // FrameAcquisitionDateTime
938   ReplaceIf(ds, pffgs_ds, Tag(0x0020,0x9111), Tag(0x0018,0x9074) );
939   // Nominal Cardiac Trigger Delay Time -> Trigger Time
940   //const DataElement &NominalCardiacTriggerDelayTime =
941   //  GetNestedDataElement(pffgs_ds, Tag(0x0018,0x9226), Tag(0x0008,0x9007) );
942 #endif
943 
944   // (0020,9228) UL 158                                      #   4, 1 ConcatenationFrameOffsetNumber
945   gdcm::Attribute<0x0020,0x9228> at = { index };
946   ds.Replace( at.GetAsDataElement() );
947 
948   return true;
949 }
950 
951 } // namespace gdcm
952 
main(int argc,char * argv[])953 int main (int argc, char *argv[])
954 {
955   int c;
956   //int digit_optind = 0;
957 
958   int rootuid = 0;
959   std::string filename;
960   std::string outfilename;
961   std::string root;
962   int resourcespath = 0;
963   int mosaic = 0;
964   int mosaic_private = 0;
965   int enhance = 1;
966   int unenhance = 0;
967   std::string xmlpath;
968 
969   int verbose = 0;
970   int warning = 0;
971   int debug = 0;
972   int error = 0;
973   int help = 0;
974   int version = 0;
975 
976   std::string pattern;
977   while (true) {
978     //int this_option_optind = optind ? optind : 1;
979     int option_index = 0;
980     static struct option long_options[] = {
981         {"input", 1, nullptr, 0},
982         {"output", 1, nullptr, 0},
983         {"mosaic", 0, &mosaic, 1}, // split siemens mosaic into multiple frames
984         {"pattern", 1, nullptr, 0},               // p
985         {"enhance", 0, &enhance, 1},               // unenhance
986         {"unenhance", 0, &unenhance, 1},               // unenhance
987         {"root-uid", 1, &rootuid, 1}, // specific Root (not GDCM)
988         //{"resources-path", 0, &resourcespath, 1},
989         {"mosaic-private", 0, &mosaic_private, 1}, // keep private attributes
990 
991 // General options !
992         {"verbose", 0, &verbose, 1},
993         {"warning", 0, &warning, 1},
994         {"debug", 0, &debug, 1},
995         {"error", 0, &error, 1},
996         {"help", 0, &help, 1},
997         {"version", 0, &version, 1},
998 
999         {nullptr, 0, nullptr, 0}
1000     };
1001 
1002     c = getopt_long (argc, argv, "i:o:MUp:VWDEhv",
1003       long_options, &option_index);
1004     if (c == -1)
1005       {
1006       break;
1007       }
1008 
1009     switch (c)
1010       {
1011     case 0:
1012         {
1013         const char *s = long_options[option_index].name; (void)s;
1014         //printf ("option %s", s);
1015         if (optarg)
1016           {
1017           if( option_index == 0 ) /* input */
1018             {
1019             assert( strcmp(s, "input") == 0 );
1020             assert( filename.empty() );
1021             filename = optarg;
1022             }
1023           else if( option_index == 3 ) /* pattern */
1024             {
1025             assert( strcmp(s, "pattern") == 0 );
1026             assert( pattern.empty() );
1027             pattern = optarg;
1028             }
1029            else if( option_index == 6 ) /* root uid */
1030             {
1031             assert( strcmp(s, "root-uid") == 0 );
1032             root = optarg;
1033             }
1034             else if( option_index == 7 ) /* resourcespath */
1035             {
1036             assert( strcmp(s, "resources-path") == 0 );
1037             assert( xmlpath.empty() );
1038             xmlpath = optarg;
1039             }
1040           else
1041             {
1042             printf (" with arg %s, index = %d", optarg, option_index);
1043             assert(0);
1044             }
1045           //printf (" with arg %s, index = %d", optarg, option_index);
1046           }
1047         //printf ("\n");
1048         }
1049       break;
1050 
1051     case 'i':
1052       assert( filename.empty() );
1053       filename = optarg;
1054       break;
1055 
1056     case 'o':
1057       assert( outfilename.empty() );
1058       outfilename = optarg;
1059       break;
1060 
1061     case 'U':
1062       //assert( outfilename.empty() );
1063       //outfilename = optarg;
1064       //printf ("option unenhance \n");
1065       unenhance = 1;
1066       break;
1067 
1068     case 'M':
1069       //assert( outfilename.empty() );
1070       //outfilename = optarg;
1071       mosaic = 1;
1072       break;
1073 
1074     case 'p':
1075       assert( pattern.empty() );
1076       pattern = optarg;
1077       break;
1078 
1079     case 'V':
1080       verbose = 1;
1081       break;
1082 
1083     case 'W':
1084       warning = 1;
1085       break;
1086 
1087     case 'D':
1088       debug = 1;
1089       break;
1090 
1091     case 'E':
1092       error = 1;
1093       break;
1094 
1095     case 'h':
1096       help = 1;
1097       break;
1098 
1099     case 'v':
1100       version = 1;
1101       break;
1102 
1103     case '?':
1104       break;
1105 
1106     default:
1107       printf ("?? getopt returned character code 0%o ??\n", c);
1108       }
1109   }
1110 
1111   // For now only support one input / one output
1112   if (optind < argc)
1113     {
1114     std::vector<std::string> files;
1115     while (optind < argc)
1116       {
1117       //printf ("%s\n", argv[optind++]);
1118       files.emplace_back(argv[optind++] );
1119       }
1120     //printf ("\n");
1121     if( files.size() == 2
1122       && filename.empty()
1123       && outfilename.empty()
1124     )
1125       {
1126       filename = files[0];
1127       outfilename = files[1];
1128       }
1129     else
1130       {
1131       PrintHelp();
1132       return 1;
1133       }
1134     }
1135 
1136   if( version )
1137     {
1138     //std::cout << "version" << std::endl;
1139     PrintVersion();
1140     return 0;
1141     }
1142 
1143   if( help )
1144     {
1145     //std::cout << "help" << std::endl;
1146     PrintHelp();
1147     return 0;
1148     }
1149 
1150   if( filename.empty() )
1151     {
1152     //std::cerr << "Need input file (-i)\n";
1153     PrintHelp();
1154     return 1;
1155     }
1156   if( outfilename.empty() )
1157     {
1158     //std::cerr << "Need output file (-o)\n";
1159     PrintHelp();
1160     return 1;
1161     }
1162 
1163   // Debug is a little too verbose
1164   gdcm::Trace::SetDebug( (debug  > 0 ? true : false));
1165   gdcm::Trace::SetWarning(  (warning  > 0 ? true : false));
1166   gdcm::Trace::SetError(  (error  > 0 ? true : false));
1167   // when verbose is true, make sure warning+error are turned on:
1168   if( verbose )
1169     {
1170     gdcm::Trace::SetWarning( (verbose  > 0 ? true : false) );
1171     gdcm::Trace::SetError( (verbose  > 0 ? true : false) );
1172     }
1173 
1174   gdcm::FileMetaInformation::SetSourceApplicationEntityTitle( "gdcmtar" );
1175   if( !rootuid )
1176     {
1177     // only read the env var is no explicit cmd line option
1178     // maybe there is an env var defined... let's check
1179     const char *rootuid_env = getenv("GDCM_ROOT_UID");
1180     if( rootuid_env )
1181       {
1182       rootuid = 1;
1183       root = rootuid_env;
1184       }
1185     }
1186   if( rootuid )
1187     {
1188     if( !gdcm::UIDGenerator::IsValid( root.c_str() ) )
1189       {
1190       std::cerr << "specified Root UID is not valid: " << root << std::endl;
1191       return 1;
1192       }
1193     gdcm::UIDGenerator::SetRoot( root.c_str() );
1194     }
1195 
1196   if( unenhance && false )
1197     {
1198     gdcm::Global& g = gdcm::Global::GetInstance();
1199     // First thing we need to locate the XML dict
1200     // did the user requested to look XML file in a particular directory ?
1201     if( !resourcespath )
1202       {
1203       const char *xmlpathenv = getenv("GDCM_RESOURCES_PATH");
1204       if( xmlpathenv )
1205         {
1206         // Make sure to look for XML dict in user explicitly specified dir first:
1207         xmlpath = xmlpathenv;
1208         resourcespath = 1;
1209         }
1210       }
1211     if( resourcespath )
1212       {
1213       // xmlpath is set either by the cmd line option or the env var
1214       if( !g.Prepend( xmlpath.c_str() ) )
1215         {
1216         std::cerr << "specified Resources Path is not valid: " << xmlpath << std::endl;
1217         return 1;
1218         }
1219       }
1220 
1221     // All set, then load the XML files:
1222     if( !g.LoadResourcesFiles() )
1223       {
1224       return 1;
1225       }
1226 
1227     //const gdcm::Defs &defs = g.GetDefs();
1228     }
1229 
1230 
1231   if( mosaic )
1232     {
1233     gdcm::ImageReader reader;
1234     reader.SetFileName( filename.c_str() );
1235     if( !reader.Read() )
1236       {
1237       std::cerr << "could not read: " << filename << std::endl;
1238       return 1;
1239       }
1240 
1241     gdcm::SplitMosaicFilter filter;
1242     filter.SetImage( reader.GetImage() );
1243     filter.SetFile( reader.GetFile() );
1244     bool b = filter.Split();
1245     if( !b )
1246       {
1247       std::cerr << "Could not split : " << filename << std::endl;
1248       return 1;
1249       }
1250 
1251     const gdcm::Image &image = filter.GetImage();
1252     const unsigned int *dims = image.GetDimensions();
1253     const gdcm::DataElement &pixeldata = image.GetDataElement();
1254     const gdcm::ByteValue *bv = pixeldata.GetByteValue();
1255     size_t slice_len = image.GetBufferLength() / dims[2];
1256 
1257     gdcm::FilenameGenerator fg;
1258     fg.SetNumberOfFilenames( dims[2] );
1259     fg.SetPrefix( outfilename.c_str() );
1260     fg.SetPattern( pattern.c_str() );
1261     if( !fg.Generate() )
1262       {
1263       std::cerr << "could not generate filenames" << std::endl;
1264       return 1;
1265       }
1266     const double *cosines = image.GetDirectionCosines();
1267     gdcm::DirectionCosines dc( cosines );
1268     double normal[3];
1269     dc.Cross( normal );
1270     const double *origin = image.GetOrigin();
1271     double zspacing = image.GetSpacing(2);
1272 
1273     gdcm::CSAHeader csa;
1274     gdcm::DataSet & ds = reader.GetFile().GetDataSet();
1275 
1276     gdcm::MrProtocol mrprot;
1277     if( !csa.GetMrProtocol(ds, mrprot) ) return 1;
1278 
1279     gdcm::MrProtocol::SliceArray sa;
1280     b = mrprot.GetSliceArray(sa);
1281     if( !b ) return 1;
1282 
1283     size_t size = sa.Slices.size();
1284     if( !size ) return 1;
1285 
1286     if( !mosaic_private )
1287     {
1288       gdcm::Anonymizer ano;
1289       ano.SetFile( reader.GetFile() );
1290       // Remove CSA header
1291       ano.RemovePrivateTags();
1292     }
1293 
1294     double slicePos[3];
1295     double sliceNor[3];
1296     namespace kwd = gdcm::Keywords;
1297     gdcm::UIDGenerator ug;
1298 
1299     kwd::InstanceNumber instart;
1300     instart.Set(ds);
1301     int istart = instart.GetValue();
1302 
1303     for(unsigned int i = 0; i < dims[2]; ++i)
1304       {
1305       gdcm::MrProtocol::Slice & protSlice = sa.Slices[i];
1306       gdcm::MrProtocol::Vector3 & protV = protSlice.Position;
1307       gdcm::MrProtocol::Vector3 & protN = protSlice.Normal;
1308       slicePos[0] = protV.dSag;
1309       slicePos[1] = protV.dCor;
1310       slicePos[2] = protV.dTra;
1311       sliceNor[0] = protN.dSag;
1312       sliceNor[1] = protN.dCor;
1313       sliceNor[2] = protN.dTra;
1314 
1315       double new_origin[3];
1316       for (int j = 0; j < 3; j++)
1317         {
1318         // the n'th slice is n * z-spacing aloung the IOP-derived
1319         // z-axis
1320         new_origin[j] = origin[j] + normal[j] * i * zspacing;
1321         if( std::fabs(slicePos[j] - new_origin[j]) > 1e-3 )
1322           {
1323           gdcmErrorMacro("Invalid position found");
1324           return 1;
1325           }
1326         const double snv_dot = gdcm::DirectionCosines::Dot( normal, sliceNor );
1327         if( std::fabs(1. - snv_dot) > 1e-6 )
1328           {
1329           gdcmErrorMacro("Invalid direction found");
1330           return 1;
1331           }
1332         }
1333 
1334       kwd::SOPInstanceUID sid;
1335       sid.SetValue( ug.Generate() );
1336       ds.Replace( sid.GetAsDataElement() );
1337 
1338       const char *outfilenamei = fg.GetFilename(i);
1339       kwd::SliceLocation sl;
1340       sl.SetValue( new_origin[2] );
1341       ds.Replace( sl.GetAsDataElement() );
1342       kwd::InstanceNumber in;
1343       in.SetValue( istart + i ); // Start at mosaic instance number
1344       ds.Replace( in.GetAsDataElement() );
1345       gdcm::ImageWriter writer;
1346       writer.SetFileName( outfilenamei );
1347       //writer.SetFile( filter.GetFile() );
1348       writer.SetFile( reader.GetFile() );
1349 
1350       //
1351       //writer.SetImage( filter.GetImage() );
1352       gdcm::Image &slice = writer.GetImage();
1353       slice = filter.GetImage();
1354       slice.SetOrigin( new_origin );
1355       slice.SetNumberOfDimensions( 2 );
1356       assert( slice.GetPixelFormat() == filter.GetImage().GetPixelFormat() );
1357       slice.SetSpacing(2, filter.GetImage().GetSpacing(2) );
1358       //slice.Print( std::cout );
1359       gdcm::DataElement &pd = slice.GetDataElement();
1360       const char *sliceptr = bv->GetPointer() + i * slice_len;
1361       pd.SetByteValue( sliceptr, (uint32_t)slice_len);
1362 
1363       if( !writer.Write() )
1364         {
1365         std::cerr << "Failed to write: " << outfilename << std::endl;
1366         return 1;
1367         }
1368       }
1369 
1370     return 0;
1371     }
1372   else if ( unenhance )
1373     {
1374     gdcm::ImageReader reader;
1375     reader.SetFileName( filename.c_str() );
1376     if( !reader.Read() )
1377       {
1378       std::cerr << "could not read: " << filename << std::endl;
1379       return 1;
1380       }
1381 
1382     gdcm::File &file = reader.GetFile();
1383     gdcm::DataSet &ds = file.GetDataSet();
1384     gdcm::MediaStorage ms;
1385     ms.SetFromFile(file);
1386     if( ms.IsUndefined() )
1387       {
1388       std::cerr << "Unknown MediaStorage" << std::endl;
1389       return 1;
1390       }
1391 
1392     gdcm::UIDs uid;
1393     uid.SetFromUID( ms.GetString() );
1394 
1395     if( uid != gdcm::UIDs::uid_1_2_840_10008_5_1_4_1_1_4_1 ) //NOTE: uid_1_2_840_10008_5_1_4_1_1_4_1 = 121, // Enhanced MR Image Storage
1396       {
1397       std::cerr << "MediaStorage is not handled " << ms << " [" << uid.GetName() << "]" << std::endl;
1398       return 1;
1399       }
1400 
1401   // Preserve info:
1402   gdcm::DataElement oldsopclassuid = ds.GetDataElement( gdcm::Tag(0x8,0x16) );
1403   gdcm::DataElement oldinstanceuid = ds.GetDataElement( gdcm::Tag(0x8,0x18) );
1404 
1405   // Ok then change it old Old MR Image Storage
1406     gdcm::DataElement de( gdcm::Tag(0x0008, 0x0016) );
1407     ms = gdcm::MediaStorage::MRImageStorage;
1408     const char* msstr = gdcm::MediaStorage::GetMSString(ms);
1409     de.SetByteValue( msstr, (uint32_t)strlen(msstr) );
1410     de.SetVR( gdcm::Attribute<0x0008, 0x0016>::GetVR() );
1411     ds.Replace( de );
1412 
1413     const gdcm::Image &image = reader.GetImage();
1414     const unsigned int *dims = image.GetDimensions();
1415     //std::cout << image << std::endl;
1416     const gdcm::DataElement &pixeldata = image.GetDataElement();
1417     //const gdcm::ByteValue *bv = pixeldata.GetByteValue();
1418     gdcm::SmartPointer<gdcm::ByteValue> bv = const_cast<gdcm::ByteValue*>(pixeldata.GetByteValue());
1419     if( !bv )
1420     {
1421       std::cerr << "decompress first" << std::endl;
1422       return 1;
1423     }
1424     unsigned long slice_len = image.GetBufferLength() / dims[2];
1425     assert( slice_len * dims[2] == image.GetBufferLength() );
1426     //assert( image.GetBufferLength() == bv->GetLength() );
1427 
1428     gdcm::FilenameGenerator fg;
1429     fg.SetNumberOfFilenames( dims[2] );
1430     fg.SetPrefix( outfilename.c_str() );
1431     fg.SetPattern( pattern.c_str() );
1432     if( !fg.Generate() )
1433       {
1434       std::cerr << "could not generate" << std::endl;
1435       return 1;
1436       }
1437     const double *cosines = image.GetDirectionCosines();
1438     gdcm::DirectionCosines dc( cosines );
1439     double normal[3];
1440     dc.Cross( normal );
1441     //const double *origin = image.GetOrigin();
1442     //double zspacing = image.GetSpacing(2);
1443 
1444     // Remove SharedFunctionalGroupsSequence
1445     gdcm::SmartPointer<gdcm::SequenceOfItems> sfgs =
1446       ds.GetDataElement( gdcm::Tag( 0x5200,0x9229 ) ).GetValueAsSQ();
1447     ds.Remove( gdcm::Tag( 0x5200,0x9229 ) );
1448     assert( ds.FindDataElement( gdcm::Tag( 0x5200,0x9229 ) ) == false );
1449     // Remove PerFrameFunctionalGroupsSequence
1450     gdcm::SmartPointer<gdcm::SequenceOfItems> pffgs =
1451       ds.GetDataElement( gdcm::Tag( 0x5200,0x9230 ) ).GetValueAsSQ();
1452     ds.Remove( gdcm::Tag( 0x5200,0x9230 ) );
1453     assert( ds.FindDataElement( gdcm::Tag( 0x5200,0x9230 ) ) == false );
1454     ds.Remove( gdcm::Tag( 0x28,0x8) );
1455     ds.Remove( gdcm::Tag( 0x7fe0,0x0010) );
1456     assert( ds.FindDataElement( gdcm::Tag( 0x7fe0,0x0010) ) == false );
1457     //ds.Remove( gdcm::Tag( 0x0008,0x0012) );
1458     //ds.Remove( gdcm::Tag( 0x0008,0x0013) );
1459 
1460   // reference the old instance:
1461   // PS 3.3-2009 C.7.6.16.1.3
1462 #if 0
1463   assert( ds.FindDataElement( gdcm::Tag(0x0008,0x1150) ) == false );
1464   assert( ds.FindDataElement( gdcm::Tag(0x0008,0x1155) ) == false );
1465   assert( ds.FindDataElement( gdcm::Tag(0x0008,0x1160) ) == false );
1466   oldsopclassuid.SetTag( gdcm::Tag(0x8,0x1150) );
1467   oldinstanceuid.SetTag( gdcm::Tag(0x8,0x1155) );
1468   ds.Insert( oldsopclassuid );
1469   ds.Insert( oldinstanceuid );
1470 #endif
1471 
1472   char date[22];
1473   const size_t datelen = 8;
1474   gdcm::System::GetCurrentDateTime(date);
1475   gdcm::Attribute<0x8,0x12> instcreationdate;
1476   instcreationdate.SetValue( gdcm::DTComp( date, datelen ) );
1477   ds.Replace( instcreationdate.GetAsDataElement() );
1478   gdcm::Attribute<0x8,0x13> instcreationtime;
1479   instcreationtime.SetValue( gdcm::DTComp( date + datelen, 13 ) );
1480   ds.Replace( instcreationtime.GetAsDataElement() );
1481   const char *offset = gdcm::System::GetTimezoneOffsetFromUTC();
1482   gdcm::Attribute<0x8,0x201> timezoneoffsetfromutc;
1483   timezoneoffsetfromutc.SetValue( offset );
1484   ds.Replace( timezoneoffsetfromutc.GetAsDataElement() );
1485 
1486     for(unsigned int i = 0; i < dims[2]; ++i)
1487       {
1488 #if 0
1489       double new_origin[3];
1490       for (int j = 0; j < 3; j++)
1491         {
1492         // the n'th slice is n * z-spacing aloung the IOP-derived
1493         // z-axis
1494         new_origin[j] = origin[j] + normal[j] * i * zspacing;
1495         }
1496 #endif
1497 
1498       const char *outfilenamei = fg.GetFilename(i);
1499       //gdcm::ImageWriter writer;
1500       gdcm::Writer writer;
1501       writer.SetFileName( outfilenamei );
1502       //writer.SetFile( filter.GetFile() );
1503       writer.SetFile( reader.GetFile() );
1504 
1505       if ( !gdcm::RemapSharedIntoOld( ds, sfgs, pffgs, i ) )
1506         {
1507         return 1;
1508         }
1509 
1510       //
1511       //writer.SetImage( filter.GetImage() );
1512       //gdcm::Image & //slice = writer.GetImage();
1513       //slice = reader.GetImage();
1514 //      slice.SetOrigin( new_origin );
1515 //      slice.SetNumberOfDimensions( 2 );
1516 //      assert( slice.GetPixelFormat() == reader.GetImage().GetPixelFormat() );
1517 //      slice.SetSpacing(2, reader.GetImage().GetSpacing(2) );
1518       //slice.Print( std::cout );
1519 //      gdcm::DataElement &pd = slice.GetDataElement();
1520       const char *sliceptr = bv->GetPointer() + i * slice_len;
1521       gdcm::DataElement newpixeldata( gdcm::Tag(0x7fe0,0x0010) );
1522       newpixeldata.SetVR( pixeldata.GetVR() );
1523       newpixeldata.SetByteValue( sliceptr, (uint32_t)slice_len); // slow !
1524       ds.Replace( newpixeldata );
1525 
1526       if( !writer.Write() )
1527         {
1528         std::cerr << "Failed to write: " << outfilenamei << std::endl;
1529         return 1;
1530         }
1531       //else
1532       //  {
1533       //  std::cout << "Success to write: " << outfilenamei << std::endl;
1534       //  }
1535       }
1536 
1537     return 0;
1538     }
1539   else
1540     {
1541     assert( enhance );
1542     return MakeImageEnhanced( filename, outfilename );
1543 #if 0
1544     std::cerr << "Not implemented" << std::endl;
1545     return 1;
1546     gdcm::ImageReader reader;
1547     reader.SetFileName( filename.c_str() );
1548     if( !reader.Read() )
1549       {
1550       std::cerr << "could not read: " << filename << std::endl;
1551       return 1;
1552       }
1553 
1554     const gdcm::Image &image = reader.GetImage();
1555     const unsigned int *dims = image.GetDimensions();
1556     std::cout << image << std::endl;
1557     const gdcm::DataElement &pixeldata = image.GetDataElement();
1558     const gdcm::ByteValue *bv = pixeldata.GetByteValue();
1559     unsigned long slice_len = image.GetBufferLength() / dims[2];
1560     //assert( image.GetBufferLength() == bv->GetLength() );
1561 
1562     gdcm::FilenameGenerator fg;
1563     fg.SetNumberOfFilenames( dims[2] );
1564     fg.SetPrefix( outfilename.c_str() );
1565     fg.SetPattern( pattern.c_str() );
1566     if( !fg.Generate() )
1567       {
1568       std::cerr << "could not generate" << std::endl;
1569       return 1;
1570       }
1571     const double *cosines = image.GetDirectionCosines();
1572     gdcm::DirectionCosines dc( cosines );
1573     double normal[3];
1574     dc.Cross( normal );
1575     const double *origin = image.GetOrigin();
1576     double zspacing = image.GetSpacing(2);
1577 
1578     for(unsigned int i = 0; i < dims[2]; ++i)
1579       {
1580       double new_origin[3];
1581       for (int j = 0; j < 3; j++)
1582         {
1583         // the n'th slice is n * z-spacing aloung the IOP-derived
1584         // z-axis
1585         new_origin[j] = origin[j] + normal[j] * i * zspacing;
1586         }
1587 
1588       const char *outfilenamei = fg.GetFilename(i);
1589       gdcm::ImageWriter writer;
1590       writer.SetFileName( outfilenamei );
1591       //writer.SetFile( filter.GetFile() );
1592       writer.SetFile( reader.GetFile() );
1593 
1594       //
1595       //writer.SetImage( filter.GetImage() );
1596       gdcm::Image &slice = writer.GetImage();
1597       slice = reader.GetImage();
1598       slice.SetOrigin( new_origin );
1599       slice.SetNumberOfDimensions( 2 );
1600       assert( slice.GetPixelFormat() == reader.GetImage().GetPixelFormat() );
1601       slice.SetSpacing(2, reader.GetImage().GetSpacing(2) );
1602       //slice.Print( std::cout );
1603       gdcm::DataElement &pd = slice.GetDataElement();
1604       const char *sliceptr = bv->GetPointer() + i * slice_len;
1605       pd.SetByteValue( sliceptr, slice_len); // slow !
1606 
1607       if( !writer.Write() )
1608         {
1609         std::cerr << "Failed to write: " << outfilenamei << std::endl;
1610         return 1;
1611         }
1612       else
1613         {
1614         std::cout << "Success to write: " << outfilenamei << std::endl;
1615         }
1616       }
1617 
1618     return 0;
1619 #endif
1620     }
1621 }
1622