1 /*============================================================================
2   MetaIO
3   Copyright 2000-2010 Insight Software Consortium
4 
5   Distributed under the OSI-approved BSD License (the "License");
6   see accompanying file Copyright.txt for details.
7 
8   This software is distributed WITHOUT ANY WARRANTY; without even the
9   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
10   See the License for more information.
11 ============================================================================*/
12 #include "metaImage.h"
13 
14 #include <algorithm> // for std::min & std::max
15 #include <cctype>
16 #include <cmath>
17 #include <cstdio>
18 #include <cstdlib> // for atoi
19 #include <cstring> // for memcpy
20 #include <string>
21 
22 #if defined (__BORLANDC__) && (__BORLANDC__ >= 0x0580)
23 #include <mem.h>
24 #endif
25 
26 #if defined(_WIN32)
27 # include <io.h>
28 #endif
29 
30 // support for access
31 #ifndef _WIN32
32 #include <climits>
33 #include <csignal>    /* sigprocmask */
34 #include <pwd.h>
35 #include <sys/ioctl.h>
36 #include <sys/param.h>
37 #include <sys/wait.h>
38 #include <termios.h>
39 #include <unistd.h>
40 #endif
41 
42 namespace {
43 
openReadStream(std::ifstream & inputStream,const std::string & fname)44 void openReadStream(std::ifstream & inputStream, const std::string& fname)
45 {
46 #ifdef __sgi
47   inputStream.open( fname, std::ios::in );
48 #else
49   inputStream.open( fname, std::ios::in |
50                            std::ios::binary );
51 #endif
52 }
53 
openWriteStream(std::ofstream & outputStream,const std::string & fname,bool append)54 void openWriteStream(std::ofstream & outputStream, const std::string& fname, bool append)
55 {
56 // Some older sgi compilers have a error in the ofstream constructor
57 // that requires a file to exist for output
58 #ifdef __sgi
59 {
60   std::ofstream tFile(fname, std::ios::out);
61   tFile.close();
62 }
63 #endif
64 
65   if(!append)
66     {
67 #ifdef __sgi
68     outputStream.open(fname, std::ios::out);
69 #else
70     outputStream.open(fname, std::ios::binary |
71                              std::ios::out);
72 #endif
73     }
74   else
75     {
76 #ifdef __sgi
77     outputStream.open(fname, std::ios::app |
78                              std::ios::out);
79 #else
80     outputStream.open(fname, std::ios::binary |
81                              std::ios::app |
82                              std::ios::out);
83 #endif
84     }
85 }
86 
87 } // end anonymous namespace
88 
89 #if (METAIO_USE_NAMESPACE)
90 namespace METAIO_NAMESPACE {
91 #endif
92 
93 // 1 Gigabyte is the maximum chunk to read/write in on function call
94 static const std::streamoff MaxIOChunk = 1024*1024*1024;
95 
96 //
97 // MetaImage Constructors
98 //
99 MetaImage::
MetaImage()100 MetaImage()
101 :MetaObject()
102 {
103   if(META_DEBUG)
104     {
105     std::cout << "MetaImage()" << std::endl;
106     }
107 
108   m_CompressionTable = new MET_CompressionTableType;
109   m_CompressionTable->compressedStream = nullptr;
110   m_CompressionTable->buffer = nullptr;
111   Clear();
112 }
113 
114 //
115 MetaImage::
MetaImage(const char * _headerName)116 MetaImage(const char *_headerName)
117 :MetaObject()
118 {
119   if(META_DEBUG)
120     {
121     std::cout << "MetaImage()" << std::endl;
122     }
123 
124   m_CompressionTable = new MET_CompressionTableType;
125   m_CompressionTable->compressedStream = nullptr;
126   m_CompressionTable->buffer = nullptr;
127   Clear();
128 
129   Read(_headerName);
130 }
131 
132 //
133 MetaImage::
MetaImage(MetaImage * _im)134 MetaImage(MetaImage *_im)
135 :MetaObject()
136 {
137   if(META_DEBUG)
138     {
139     std::cout << "MetaImage()" << std::endl;
140     }
141 
142   m_CompressionTable = new MET_CompressionTableType;
143   m_CompressionTable->compressedStream = nullptr;
144   m_CompressionTable->buffer = nullptr;
145   Clear();
146 
147   InitializeEssential(_im->NDims(),
148                       _im->DimSize(),
149                       _im->ElementSpacing(),
150                       _im->ElementType(),
151                       _im->ElementNumberOfChannels(),
152                       _im->ElementData(),
153                       false);
154   CopyInfo(_im);
155 }
156 
157 //
158 void MetaImage::
InitHelper(int _nDims,const int * _dimSize,const double * _elementSpacing,MET_ValueEnumType _elementType,int _elementNumberOfChannels,void * _elementData)159 InitHelper(int _nDims,
160           const int * _dimSize,
161           const double * _elementSpacing,
162           MET_ValueEnumType _elementType,
163           int _elementNumberOfChannels,
164           void *_elementData)
165 {
166   if(META_DEBUG)
167     {
168     std::cout << "MetaImage()" << std::endl;
169     }
170 
171   m_CompressionTable = new MET_CompressionTableType;
172   m_CompressionTable->buffer = nullptr;
173   m_CompressionTable->compressedStream = nullptr;
174   Clear();
175 
176   if(_elementData == nullptr)
177     {
178     InitializeEssential(_nDims,
179                         _dimSize,
180                         _elementSpacing,
181                         _elementType,
182                         _elementNumberOfChannels,
183                         nullptr, true);
184     }
185   else
186     {
187     InitializeEssential(_nDims,
188                         _dimSize,
189                         _elementSpacing,
190                         _elementType,
191                         _elementNumberOfChannels,
192                         _elementData, false);
193     }
194 
195 }
196 
197 //
198 MetaImage::
MetaImage(int _nDims,const int * _dimSize,const float * _elementSpacing,MET_ValueEnumType _elementType,int _elementNumberOfChannels,void * _elementData)199 MetaImage(int _nDims,
200           const int * _dimSize,
201           const float * _elementSpacing,
202           MET_ValueEnumType _elementType,
203           int _elementNumberOfChannels,
204           void *_elementData)
205 :MetaObject()
206 {
207   // Only consider at most 10 element of spacing:
208   // See MetaObject::InitializeEssential(_nDims)
209   double tmpElementSpacing[10];
210   int ndims = std::max( std::min( _nDims, 10 ), 0);
211   for( int i = 0; i < ndims; ++i )
212     {
213     tmpElementSpacing[i] = static_cast<double>(_elementSpacing[i]);
214     }
215    InitHelper(_nDims, _dimSize, tmpElementSpacing, _elementType,
216      _elementNumberOfChannels, _elementData);
217 }
218 
219 //
220 MetaImage::
MetaImage(int _nDims,const int * _dimSize,const double * _elementSpacing,MET_ValueEnumType _elementType,int _elementNumberOfChannels,void * _elementData)221 MetaImage(int _nDims,
222           const int * _dimSize,
223           const double * _elementSpacing,
224           MET_ValueEnumType _elementType,
225           int _elementNumberOfChannels,
226           void *_elementData)
227 :MetaObject()
228 {
229   InitHelper(_nDims, _dimSize, _elementSpacing, _elementType,
230     _elementNumberOfChannels, _elementData);
231 }
232 
233 //
234 MetaImage::
MetaImage(int _x,int _y,double _elementSpacingX,double _elementSpacingY,MET_ValueEnumType _elementType,int _elementNumberOfChannels,void * _elementData)235 MetaImage(int _x, int _y,
236           double _elementSpacingX, double _elementSpacingY,
237           MET_ValueEnumType _elementType,
238           int _elementNumberOfChannels, void *_elementData)
239 :MetaObject()
240 {
241   if(META_DEBUG)
242     {
243     std::cout << "MetaImage()" << std::endl;
244     }
245 
246   m_CompressionTable = new MET_CompressionTableType;
247   m_CompressionTable->compressedStream = nullptr;
248   m_CompressionTable->buffer = nullptr;
249   Clear();
250 
251   int ds[2];
252   ds[0] = _x;
253   ds[1] = _y;
254 
255   double es[2];
256   es[0] = _elementSpacingX;
257   es[1] = _elementSpacingY;
258 
259   if(_elementData == nullptr)
260     {
261     InitializeEssential(2,
262                         ds,
263                         es,
264                         _elementType,
265                         _elementNumberOfChannels,
266                         nullptr,
267                         true);
268     }
269   else
270     {
271     InitializeEssential(2,
272                         ds,
273                         es,
274                         _elementType,
275                         _elementNumberOfChannels,
276                         _elementData,
277                         false);
278     }
279 }
280 
281 //
282 MetaImage::
MetaImage(int _x,int _y,int _z,double _elementSpacingX,double _elementSpacingY,double _elementSpacingZ,MET_ValueEnumType _elementType,int _elementNumberOfChannels,void * _elementData)283 MetaImage(int _x, int _y, int _z,
284           double _elementSpacingX,
285           double _elementSpacingY,
286           double _elementSpacingZ,
287           MET_ValueEnumType _elementType,
288           int _elementNumberOfChannels,
289           void *_elementData)
290 :MetaObject()
291 {
292   if(META_DEBUG)
293     {
294     std::cout << "MetaImage()" << std::endl;
295     }
296 
297   m_CompressionTable = new MET_CompressionTableType;
298   m_CompressionTable->compressedStream = nullptr;
299   m_CompressionTable->buffer = nullptr;
300   Clear();
301 
302   int ds[3];
303   ds[0] = _x;
304   ds[1] = _y;
305   ds[2] = _z;
306 
307   double es[3];
308   es[0] = _elementSpacingX;
309   es[1] = _elementSpacingY;
310   es[2] = _elementSpacingZ;
311 
312   if(_elementData == nullptr)
313     {
314     InitializeEssential(3,
315                         ds,
316                         es,
317                         _elementType,
318                         _elementNumberOfChannels,
319                         nullptr,
320                         true);
321     }
322   else
323     {
324     InitializeEssential(3,
325                         ds,
326                         es,
327                         _elementType,
328                         _elementNumberOfChannels,
329                         _elementData,
330                         false);
331     }
332 }
333 
334 //
335 MetaImage::
~MetaImage()336 ~MetaImage()
337 {
338   M_Destroy();
339 }
340 
341 //
342 void MetaImage::
PrintInfo() const343 PrintInfo() const
344 {
345   int i;
346 
347   MetaObject::PrintInfo();
348 
349   std::string s;
350   MET_ImageModalityToString(m_Modality, s);
351   std::cout << "Modality = " << s << std::endl;
352 
353   std::cout << "DimSize = ";
354   for(i=0; i<m_NDims; i++)
355     {
356     std::cout << m_DimSize[i] << " ";
357     }
358   std::cout << std::endl;
359   std::cout << "SubQuantity = ";
360   for(i=0; i<m_NDims; i++)
361     {
362     std::cout << m_SubQuantity[i] << " ";
363     }
364   std::cout << std::endl;
365 
366   std::cout << "Quantity = " << m_Quantity << std::endl;
367 
368 
369   std::cout << "HeaderSize = " << m_HeaderSize << std::endl;
370 
371   std::cout << "SequenceID = ";
372   for(i=0; i<m_NDims; i++)
373     {
374     std::cout << m_SequenceID[i] << " ";
375     }
376   std::cout << std::endl;
377 
378   std::cout << "ElementSizeValid = " << (int)m_ElementSizeValid
379             << std::endl;
380   std::cout << "ElementSize = ";
381   for(i=0; i<m_NDims; i++)
382     {
383     std::cout << m_ElementSize[i] << " ";
384     }
385   std::cout << std::endl;
386 
387   char str[22];
388   MET_TypeToString(m_ElementType, str);
389   std::cout << "ElementType = " << str << std::endl;
390 
391   std::cout << "ElementNumberOfChannels = "
392                       << m_ElementNumberOfChannels << std::endl;
393 
394   if(m_ElementMinMaxValid)
395     {
396     std::cout << "Min and Max are valid" << std::endl;
397     std::cout << "   Min = " << m_ElementMin << std::endl;
398     std::cout << "   Max = " << m_ElementMax << std::endl;
399     }
400   else
401     {
402     std::cout << "Min and Max are not valid" << std::endl;
403     }
404 
405   std::cout << "ElementToIntensityFunctionSlope = "
406                       << m_ElementToIntensityFunctionSlope
407                       << std::endl;
408   std::cout << "ElementToIntensityFunctionOffset = "
409                       << m_ElementToIntensityFunctionOffset
410                       << std::endl;
411 
412 
413   std::cout << "AutoFreeElementData = "
414                       << ((m_AutoFreeElementData)?"True":"False")
415                       << std::endl;
416 
417   std::cout << "ElementData = "
418                       << ((m_ElementData==nullptr)?"NULL":"Valid")
419                       << std::endl;
420 
421   std::cout << "ElementDataFileName = "
422                       << m_ElementDataFileName << std::endl;
423 
424 }
425 
426 void MetaImage::
CopyInfo(const MetaObject * _object)427 CopyInfo(const MetaObject * _object)
428 {
429   MetaObject::CopyInfo(_object);
430 
431   if(_object)
432     {
433     const MetaImage * im;
434     try
435       {
436       im = static_cast<const MetaImage *>(_object);
437       }
438     catch( ... )
439       {
440       return;
441       }
442 
443     if( im )
444       {
445       Modality(im->Modality());
446 
447       HeaderSize(im->HeaderSize());
448 
449       SequenceID(im->SequenceID());
450 
451       ElementSizeValid(im->ElementSizeValid());
452       if(im->ElementSizeValid())
453         {
454         ElementSize(im->ElementSize());
455         }
456 
457       ElementMinMaxValid(im->ElementMinMaxValid());
458       if(im->ElementMinMaxValid())
459         {
460         ElementMin(im->ElementMin());
461         ElementMax(im->ElementMax());
462         }
463 
464       ElementToIntensityFunctionSlope(im->ElementToIntensityFunctionSlope());
465       ElementToIntensityFunctionOffset(im->ElementToIntensityFunctionOffset());
466       }
467     }
468 }
469 
470 /** Clear function */
Clear()471 void MetaImage::Clear()
472 {
473   if(META_DEBUG)
474     {
475     std::cout << "MetaImage: Clear" << std::endl;
476     }
477 
478   m_Modality = MET_MOD_UNKNOWN;
479 
480   m_DimSize[0] = 0;
481   m_SubQuantity[0] = 0;
482   m_Quantity = 0;
483 
484   m_HeaderSize = 0;
485 
486   memset(m_SequenceID, 0, sizeof(m_SequenceID));
487 
488   m_ElementSizeValid = false;
489   memset(m_ElementSize, 0, sizeof(m_ElementSize));
490 
491   m_ElementType = MET_NONE;
492 
493   m_ElementNumberOfChannels = 1;
494 
495   m_ElementMinMaxValid = false;
496   m_ElementMin = 0;
497   m_ElementMax = 0;
498 
499   m_ElementToIntensityFunctionSlope = 1;
500   m_ElementToIntensityFunctionOffset = 0;
501 
502   m_AutoFreeElementData = true;
503 
504   m_ElementData = nullptr;
505 
506   m_ElementDataFileName = "";
507 
508   MetaObject::Clear();
509 
510   strcpy(m_ObjectTypeName,"Image");
511 
512   // Change the default for this object
513   m_BinaryData = true;
514 
515   if(m_CompressionTable)
516     {
517     if(m_CompressionTable->compressedStream)
518       {
519       inflateEnd(m_CompressionTable->compressedStream);
520       delete m_CompressionTable->compressedStream;
521       delete [] m_CompressionTable->buffer;
522       m_CompressionTable->buffer = nullptr;
523       }
524     m_CompressionTable->compressedStream = nullptr;
525     m_CompressionTable->offsetList.clear();
526     }
527   else
528     {
529     m_CompressionTable = new MET_CompressionTableType;
530     m_CompressionTable->compressedStream = nullptr;
531     }
532 
533 }
534 
535 bool MetaImage::
InitializeEssential(int _nDims,const int * _dimSize,const float * _elementSpacing,MET_ValueEnumType _elementType,int _elementNumberOfChannels,void * _elementData,bool _allocElementMemory)536 InitializeEssential(int _nDims,
537                     const int * _dimSize,
538                     const float * _elementSpacing,
539                     MET_ValueEnumType _elementType,
540                     int _elementNumberOfChannels,
541                     void * _elementData,
542                     bool _allocElementMemory)
543 {
544   // Only consider at most 10 element of spacing:
545   // See MetaObject::InitializeEssential(_nDims)
546   double tmpElementSpacing[10];
547   int ndims = std::max( std::min( _nDims, 10 ), 0);
548   for( int i = 0; i < ndims; ++i )
549     {
550     tmpElementSpacing[i] = static_cast<double>(_elementSpacing[i]);
551     }
552   return InitializeEssential(_nDims, _dimSize, tmpElementSpacing, _elementType,
553     _elementNumberOfChannels, _elementData, _allocElementMemory);
554 }
555 
556 bool MetaImage::
InitializeEssential(int _nDims,const int * _dimSize,const double * _elementSpacing,MET_ValueEnumType _elementType,int _elementNumberOfChannels,void * _elementData,bool _allocElementMemory)557 InitializeEssential(int _nDims,
558                     const int * _dimSize,
559                     const double * _elementSpacing,
560                     MET_ValueEnumType _elementType,
561                     int _elementNumberOfChannels,
562                     void * _elementData,
563                     bool _allocElementMemory)
564 {
565   if(META_DEBUG)
566     {
567     std::cout << "MetaImage: Initialize" << std::endl;
568     }
569 
570   MetaObject::InitializeEssential(_nDims);
571 
572   int i;
573   if(!m_CompressionTable)
574     {
575     m_CompressionTable = new MET_CompressionTableType;
576     m_CompressionTable->buffer = nullptr;
577     m_CompressionTable->compressedStream = nullptr;
578     }
579   m_SubQuantity[0] = 1;
580   m_Quantity = 1;
581   m_ElementSizeValid = false;
582   for(i=0; i<m_NDims; i++)
583     {
584     m_DimSize[i] = _dimSize[i];
585     m_Quantity *= _dimSize[i];
586     if(i>0)
587       {
588       m_SubQuantity[i] = m_SubQuantity[i-1]*m_DimSize[i-1];
589       }
590     m_ElementSpacing[i] = _elementSpacing[i];
591     if(m_ElementSize[i] == 0)
592       {
593       m_ElementSize[i] = m_ElementSpacing[i];
594       }
595     else
596       {
597       m_ElementSizeValid = true;
598       }
599     }
600 
601   m_ElementType = _elementType;
602 
603   m_ElementNumberOfChannels = _elementNumberOfChannels;
604 
605   if(_elementData != nullptr)
606     {
607     m_AutoFreeElementData = false;
608     m_ElementData = (void *)_elementData;
609     }
610   else if(_allocElementMemory)
611     {
612     m_AutoFreeElementData = true;
613     MET_SizeOfType(m_ElementType, &i);
614     m_ElementData = new char[static_cast<size_t>(m_Quantity*m_ElementNumberOfChannels*i)];
615     if(m_ElementData == nullptr)
616       {
617       std::cerr << "MetaImage:: M_Allocate:: Insufficient memory"
618                           << std::endl;
619       return false;
620       }
621     }
622   else
623     {
624     m_AutoFreeElementData = true;
625     m_ElementData = nullptr;
626     }
627 
628   return true;
629 }
630 
631 
632 int MetaImage::
HeaderSize() const633 HeaderSize() const
634 {
635   return m_HeaderSize;
636 }
637 
638 void MetaImage::
HeaderSize(int _headerSize)639 HeaderSize(int _headerSize)
640 {
641   m_HeaderSize = _headerSize;
642 }
643 
644 MET_ImageModalityEnumType MetaImage::
Modality() const645 Modality() const
646 {
647   return m_Modality;
648 }
649 
650 void MetaImage::
Modality(MET_ImageModalityEnumType _modality)651 Modality(MET_ImageModalityEnumType _modality)
652 {
653   m_Modality = _modality;
654 }
655 
656 const int * MetaImage::
DimSize() const657 DimSize() const
658 {
659   return m_DimSize;
660 }
661 
662 int MetaImage::
DimSize(int _i) const663 DimSize(int _i) const
664 {
665   return m_DimSize[_i];
666 }
667 
668 std::streamoff MetaImage::
Quantity() const669 Quantity() const
670 {
671   return m_Quantity;
672 }
673 
674 const std::streamoff * MetaImage::
SubQuantity() const675 SubQuantity() const
676 {
677   return m_SubQuantity;
678 }
679 
680 std::streamoff MetaImage::
SubQuantity(int _i) const681 SubQuantity(int _i) const
682 {
683   return m_SubQuantity[_i];
684 }
685 
686 const float * MetaImage::
SequenceID() const687 SequenceID() const
688 {
689   return m_SequenceID;
690 }
691 
692 float MetaImage::
SequenceID(int _i) const693 SequenceID(int _i) const
694 {
695   return m_SequenceID[_i];
696 }
697 
698 void MetaImage::
SequenceID(const float * _sequenceID)699 SequenceID(const float *_sequenceID)
700 {
701   memcpy(m_SequenceID, _sequenceID, m_NDims*sizeof(float));
702 }
703 
704 void MetaImage::
SequenceID(int _i,float _value)705 SequenceID(int _i, float _value)
706 {
707   m_SequenceID[_i] = _value;
708 }
709 
710 bool MetaImage::
ElementSizeValid() const711 ElementSizeValid() const
712 {
713   return m_ElementSizeValid;
714 }
715 
716 void MetaImage::
ElementSizeValid(bool _elementSizeValid)717 ElementSizeValid(bool _elementSizeValid)
718 {
719   m_ElementSizeValid = _elementSizeValid;
720 }
721 
722 const double * MetaImage::
ElementSize() const723 ElementSize() const
724 {
725   return m_ElementSize;
726 }
727 
728 double MetaImage::
ElementSize(int _i) const729 ElementSize(int _i) const
730 {
731   return m_ElementSize[_i];
732 }
733 
734 void MetaImage::
ElementSize(const double * _elementSize)735 ElementSize(const double *_elementSize)
736 {
737   memcpy(m_ElementSize, _elementSize, m_NDims*sizeof(*m_ElementSize));
738   m_ElementSizeValid = true;
739 }
740 
741 void MetaImage::
ElementSize(const float * _elementSize)742 ElementSize(const float *_elementSize)
743 {
744   for(int i = 0; i < m_NDims; ++i)
745     {
746     m_ElementSize[i] = static_cast<double>(_elementSize[i]);
747     }
748   m_ElementSizeValid = true;
749 }
750 
751 
752 void MetaImage::
ElementSize(int _i,double _value)753 ElementSize(int _i, double _value)
754 {
755   m_ElementSize[_i] = _value;
756   m_ElementSizeValid = true;
757 }
758 
759 MET_ValueEnumType MetaImage::
ElementType() const760 ElementType() const
761 {
762   return m_ElementType;
763 }
764 
765 void MetaImage::
ElementType(MET_ValueEnumType _elementType)766 ElementType(MET_ValueEnumType _elementType)
767 {
768   m_ElementType = _elementType;
769 }
770 
771 int MetaImage::
ElementNumberOfChannels() const772 ElementNumberOfChannels() const
773 {
774   return m_ElementNumberOfChannels;
775 }
776 
777 void MetaImage::
ElementNumberOfChannels(int _elementNumberOfChannels)778 ElementNumberOfChannels(int _elementNumberOfChannels)
779 {
780   m_ElementNumberOfChannels = _elementNumberOfChannels;
781 }
782 
783 void MetaImage::
ElementByteOrderSwap(std::streamoff _quantity)784 ElementByteOrderSwap(std::streamoff _quantity)
785 {
786 
787   // use the user provided value if provided or the internal ivar
788   std::streamoff quantity = _quantity ? _quantity : m_Quantity;
789 
790   if(META_DEBUG)
791     {
792     std::cout << "MetaImage: ElementByteOrderSwap"
793                         << std::endl;
794     }
795 
796   int eSize;
797   MET_SizeOfType(m_ElementType, &eSize);
798   switch(eSize)
799     {
800     default:
801     case 0:
802     case 1:
803       {
804       break;
805       }
806     case 2:
807       {
808       for(size_t i=0;
809         i<static_cast<size_t>(quantity*m_ElementNumberOfChannels); i++)
810         {
811         ((MET_USHORT_TYPE *)m_ElementData)[i] =
812               MET_ByteOrderSwapShort(((MET_USHORT_TYPE *)m_ElementData)[i]);
813         }
814       break;
815       }
816     case 4:
817       {
818       for(size_t i=0;
819         i<static_cast<size_t>(quantity*m_ElementNumberOfChannels); i++)
820         {
821         ((MET_UINT_TYPE *)m_ElementData)[i] =
822               MET_ByteOrderSwapLong(((MET_UINT_TYPE *)m_ElementData)[i]);
823         }
824       break;
825       }
826     case 8:
827       {
828       char* data = (char*)m_ElementData;
829       for(size_t i=0;
830         i<static_cast<size_t>(quantity*m_ElementNumberOfChannels); i++)
831         {
832         MET_ByteOrderSwap8(data);
833         data += 8;
834         }
835       break;
836       }
837     }
838   m_BinaryDataByteOrderMSB = !m_BinaryDataByteOrderMSB;
839 }
840 
841 bool MetaImage::
ElementByteOrderFix(std::streamoff _quantity)842 ElementByteOrderFix(std::streamoff _quantity)
843 {
844   if(m_BinaryDataByteOrderMSB != MET_SystemByteOrderMSB())
845     {
846     ElementByteOrderSwap(_quantity);
847     return true;
848     }
849   return true;
850 }
851 
852 bool MetaImage::
ElementMinMaxValid() const853 ElementMinMaxValid() const
854 {
855   return m_ElementMinMaxValid;
856 }
857 
858 void MetaImage::
ElementMinMaxValid(bool _elementMinMaxValid)859 ElementMinMaxValid(bool _elementMinMaxValid)
860 {
861   m_ElementMinMaxValid = _elementMinMaxValid;
862 }
863 
864 void MetaImage::
ElementMinMaxRecalc()865 ElementMinMaxRecalc()
866 {
867   double tf;
868 
869   if(m_ElementData == nullptr)
870     return;
871 
872   ElementByteOrderFix();
873 
874   MET_ValueToDouble(m_ElementType, m_ElementData, 0, &tf);
875   m_ElementMin = tf;
876   m_ElementMax = tf;
877   for(size_t i=1;
878     i<static_cast<size_t>(m_Quantity*m_ElementNumberOfChannels); i++)
879     {
880     MET_ValueToDouble(m_ElementType, m_ElementData, i, &tf);
881     if(tf<m_ElementMin)
882       {
883       m_ElementMin = tf;
884       }
885     else if(tf>m_ElementMax)
886       {
887       m_ElementMax = tf;
888       }
889     }
890 
891   m_ElementMinMaxValid = true;
892 }
893 
894 double MetaImage::
ElementMin() const895 ElementMin() const
896 {
897   return m_ElementMin;
898 }
899 
900 void MetaImage::
ElementMin(double _elementMin)901 ElementMin(double _elementMin)
902 {
903   m_ElementMin = _elementMin;
904 }
905 
906 double MetaImage::
ElementMax() const907 ElementMax() const
908 {
909   return m_ElementMax;
910 }
911 
912 void MetaImage::
ElementMax(double _elementMax)913 ElementMax(double _elementMax)
914 {
915   m_ElementMax = _elementMax;
916 }
917 
918 double MetaImage::
ElementToIntensityFunctionSlope() const919 ElementToIntensityFunctionSlope() const
920 {
921   return m_ElementToIntensityFunctionSlope;
922 }
923 
924 void MetaImage::
ElementToIntensityFunctionSlope(double _elementToIntensityFunctionSlope)925 ElementToIntensityFunctionSlope(double _elementToIntensityFunctionSlope)
926 {
927   m_ElementToIntensityFunctionSlope = _elementToIntensityFunctionSlope;
928 }
929 
930 double MetaImage::
ElementToIntensityFunctionOffset() const931 ElementToIntensityFunctionOffset() const
932 {
933   return m_ElementToIntensityFunctionOffset;
934 }
935 
936 void MetaImage::
ElementToIntensityFunctionOffset(double _elementOffset)937 ElementToIntensityFunctionOffset(double _elementOffset)
938 {
939   m_ElementToIntensityFunctionOffset = _elementOffset;
940 }
941 
942 bool MetaImage::
AutoFreeElementData() const943 AutoFreeElementData() const
944 {
945   return m_AutoFreeElementData;
946 }
947 
948 void MetaImage::
AutoFreeElementData(bool _autoFreeElementData)949 AutoFreeElementData(bool _autoFreeElementData)
950 {
951   m_AutoFreeElementData = _autoFreeElementData;
952 }
953 
954 const char * MetaImage::
ElementDataFileName() const955 ElementDataFileName() const
956 {
957   return m_ElementDataFileName.c_str();
958 }
959 
960 void MetaImage::
ElementDataFileName(const char * _elementDataFileName)961 ElementDataFileName(const char * _elementDataFileName)
962 {
963   m_ElementDataFileName = _elementDataFileName;
964 }
965 
966 void * MetaImage::
ElementData()967 ElementData()
968 {
969   return m_ElementData;
970 }
971 
972 double MetaImage::
ElementData(std::streamoff _i) const973 ElementData(std::streamoff _i) const
974 {
975   double tf = 0;
976   MET_ValueToDouble(m_ElementType, m_ElementData, _i, &tf);
977 
978   return tf;
979 }
980 
981 bool MetaImage::
ElementData(std::streamoff _i,double _v)982 ElementData(std::streamoff _i, double _v)
983 {
984   if(_i<m_Quantity)
985     {
986     MET_DoubleToValue(_v, m_ElementType, m_ElementData, _i);
987     return true;
988     }
989   return false;
990 }
991 
992 void MetaImage::
ElementData(void * _elementData,bool _autoFreeElementData)993 ElementData(void * _elementData, bool _autoFreeElementData)
994 {
995   if(m_AutoFreeElementData)
996     {
997     delete [] (char *)m_ElementData;
998     }
999   m_ElementData = _elementData;
1000   m_AutoFreeElementData = _autoFreeElementData;
1001 }
1002 
1003 bool MetaImage::
ConvertElementDataTo(MET_ValueEnumType _elementType,double _toMin,double _toMax)1004 ConvertElementDataTo(MET_ValueEnumType _elementType,
1005                      double _toMin, double _toMax)
1006 {
1007   int eSize;
1008   MET_SizeOfType(_elementType, &eSize);
1009   void * newElementData = new char[ static_cast<size_t>(
1010     m_Quantity*m_ElementNumberOfChannels*eSize) ];
1011 
1012   ElementByteOrderFix();
1013   if(!ElementMinMaxValid())
1014     {
1015     ElementMinMaxRecalc();
1016     }
1017 
1018   for(size_t i=0;
1019     i<static_cast<size_t>(m_Quantity*m_ElementNumberOfChannels); i++)
1020     {
1021     MET_ValueToValue(m_ElementType, m_ElementData, i, _elementType,
1022                      newElementData, m_ElementMin, m_ElementMax,
1023                      _toMin, _toMax);
1024     }
1025 
1026   if(m_AutoFreeElementData)
1027     {
1028     delete [] (char *)m_ElementData;
1029     }
1030   m_ElementData = newElementData;
1031   m_ElementType = _elementType;
1032   m_ElementMinMaxValid = true;
1033   m_ElementMin = _toMin;
1034   m_ElementMax = _toMax;
1035   m_AutoFreeElementData = true;
1036 
1037   return true;
1038 }
1039 
1040 bool MetaImage::
ConvertElementDataToIntensityData(MET_ValueEnumType _elementType)1041 ConvertElementDataToIntensityData(MET_ValueEnumType _elementType)
1042 {
1043   ElementByteOrderFix();
1044   if(!ElementMinMaxValid())
1045     {
1046     ElementMinMaxRecalc();
1047     }
1048 
1049   double toMin = m_ElementMin + m_ElementToIntensityFunctionOffset;
1050   double toMax = (m_ElementMax-m_ElementMin)
1051                    * m_ElementToIntensityFunctionSlope
1052                    + m_ElementMin;
1053 
1054   return ConvertElementDataTo(_elementType, toMin, toMax);
1055 }
1056 
1057 bool MetaImage::
ConvertIntensityDataToElementData(MET_ValueEnumType _elementType)1058 ConvertIntensityDataToElementData(MET_ValueEnumType _elementType)
1059 {
1060   ElementByteOrderFix();
1061   if(!ElementMinMaxValid())
1062     {
1063     ElementMinMaxRecalc();
1064     }
1065 
1066   double toMin = m_ElementMin - m_ElementToIntensityFunctionOffset;
1067   double toMax = (m_ElementMax - m_ElementMin)
1068                    / m_ElementToIntensityFunctionSlope
1069                    + toMin;
1070 
1071   return ConvertElementDataTo(_elementType, toMin, toMax);
1072 }
1073 
1074 // return true if the file exists
M_FileExists(const char * filename) const1075 bool MetaImage::M_FileExists(const char* filename) const
1076 {
1077 #ifdef _MSC_VER
1078 # define access _access
1079 #endif
1080 #ifndef R_OK
1081 # define R_OK 04
1082 #endif
1083   if ( access(filename, R_OK) != 0 )
1084     {
1085     return false;
1086     }
1087   else
1088     {
1089     return true;
1090     }
1091 }
1092 
FileIsFullPath(const char * in_name) const1093 bool MetaImage::FileIsFullPath(const char* in_name) const
1094 {
1095 #if defined(_WIN32) || defined(__CYGWIN__)
1096   // On Windows, the name must be at least two characters long.
1097   if(strlen(in_name) < 2)
1098     {
1099     return false;
1100     }
1101   if(in_name[1] == ':')
1102     {
1103     return true;
1104     }
1105   if(in_name[0] == '\\')
1106     {
1107     return true;
1108     }
1109 #else
1110   // On UNIX, the name must be at least one character long.
1111   if(strlen(in_name) < 1)
1112     {
1113     return false;
1114     }
1115 #endif
1116 #if !defined(_WIN32)
1117   if(in_name[0] == '~')
1118     {
1119     return true;
1120     }
1121 #endif
1122   // On UNIX, the name must begin in a '/'.
1123   // On Windows, if the name begins in a '/', then it is a full
1124   // network path.
1125   if(in_name[0] == '/')
1126     {
1127     return true;
1128     }
1129   return false;
1130 }
1131 
1132 // Return the value of a tag
1133 std::string
1134 MetaImage::
M_GetTagValue(const std::string & buffer,const char * tag) const1135 M_GetTagValue(const std::string & buffer, const char* tag) const
1136 {
1137   size_t stringPos = buffer.find(tag);
1138   if( stringPos == std::string::npos )
1139     {
1140     return "";
1141     }
1142 
1143   size_t pos2 = buffer.find('=',stringPos);
1144   if(pos2 == std::string::npos )
1145     {
1146     pos2 = buffer.find(':',stringPos);
1147     }
1148 
1149   if(pos2 == std::string::npos )
1150     {
1151     return "";
1152     }
1153 
1154   size_t posend = buffer.find('\r',pos2);
1155   if(posend == std::string::npos )
1156     {
1157     posend = buffer.find('\n',pos2);
1158     }
1159 
1160   // Get the element data filename
1161   std::string value = "";
1162   bool firstspace = true;
1163   size_t index = pos2+1;
1164   while(index<buffer.size()
1165         && buffer[index] != '\r'
1166         && buffer[index] != '\n'
1167         )
1168     {
1169     if(buffer[index] != ' ')
1170       {
1171       firstspace = false;
1172       }
1173     if(!firstspace)
1174       {
1175       value += buffer[index];
1176       }
1177     index++;
1178     }
1179 
1180   return value;
1181 }
1182 
1183 bool MetaImage::
CanRead(const char * _headerName) const1184 CanRead(const char *_headerName) const
1185 {
1186   // First check the extension
1187   std::string fname = _headerName;
1188   if(  fname == "" )
1189     {
1190     return false;
1191     }
1192 
1193   bool extensionFound = false;
1194 
1195   std::string::size_type stringPos = fname.rfind(".mhd");
1196   if ((stringPos != std::string::npos)
1197       && (stringPos == fname.length() - 4))
1198     {
1199     extensionFound = true;
1200     }
1201 
1202   stringPos = fname.rfind(".mha");
1203   if ((stringPos != std::string::npos)
1204       && (stringPos == fname.length() - 4))
1205     {
1206     extensionFound = true;
1207     }
1208 
1209   if( !extensionFound )
1210     {
1211     return false;
1212     }
1213 
1214   // Now check the file content
1215   std::ifstream inputStream;
1216 
1217   openReadStream(inputStream, fname.c_str());
1218 
1219   if( inputStream.fail() )
1220     {
1221     return false;
1222     }
1223 
1224   char* buf = new char[8001];
1225   inputStream.read(buf,8000);
1226   unsigned long fileSize = static_cast<unsigned long>(inputStream.gcount());
1227   buf[fileSize] = 0;
1228   std::string header(buf);
1229   header.resize(fileSize);
1230   delete [] buf;
1231   inputStream.close();
1232 
1233   stringPos = header.find("NDims");
1234   if( stringPos == std::string::npos )
1235     {
1236     return false;
1237     }
1238 
1239   std::string elementDataFileName = M_GetTagValue(header,"ElementDataFile");
1240 
1241   return true;
1242 }
1243 
1244 bool MetaImage::
Read(const char * _headerName,bool _readElements,void * _buffer)1245 Read(const char *_headerName, bool _readElements, void * _buffer)
1246 {
1247   M_Destroy();
1248 
1249   Clear();
1250 
1251   M_SetupReadFields();
1252 
1253   if(_headerName != nullptr)
1254     {
1255     m_FileName = _headerName;
1256     }
1257 
1258   M_PrepareNewReadStream();
1259 
1260   std::ifstream * tmpReadStream = new std::ifstream;
1261 
1262   openReadStream(*tmpReadStream, m_FileName);
1263 
1264   if(!tmpReadStream->is_open())
1265     {
1266     delete tmpReadStream;
1267     return false;
1268     }
1269 
1270   if( !this->ReadStream(0, tmpReadStream, _readElements, _buffer) )
1271     {
1272     tmpReadStream->close();
1273     delete tmpReadStream;
1274     return false;
1275     }
1276 
1277   tmpReadStream->close();
1278 
1279   delete tmpReadStream;
1280 
1281   return true;
1282 }
1283 
1284 bool MetaImage::
CanReadStream(std::ifstream * _stream) const1285 CanReadStream(std::ifstream * _stream) const
1286 {
1287   if(!strncmp(MET_ReadType(*_stream).c_str(), "Image", 5))
1288     {
1289     return true;
1290     }
1291   return false;
1292 }
1293 
1294 
1295 bool MetaImage::
ReadStream(int _nDims,std::ifstream * _stream,bool _readElements,void * _buffer)1296 ReadStream(int _nDims,
1297            std::ifstream * _stream,
1298            bool _readElements,
1299            void * _buffer)
1300 {
1301   if(!MetaObject::ReadStream(_nDims, _stream))
1302     {
1303     std::cerr << "MetaImage: Read: Cannot parse file"
1304                         << std::endl;
1305     return false;
1306     }
1307 
1308   if(_readElements)
1309     {
1310     if(_buffer == nullptr)
1311       {
1312       InitializeEssential(m_NDims,
1313                           m_DimSize,
1314                           m_ElementSpacing,
1315                           m_ElementType,
1316                           m_ElementNumberOfChannels,
1317                           nullptr, true);
1318       }
1319     else
1320       {
1321       InitializeEssential(m_NDims,
1322                           m_DimSize,
1323                           m_ElementSpacing,
1324                           m_ElementType,
1325                           m_ElementNumberOfChannels,
1326                           _buffer, false);
1327       }
1328 
1329     int i;
1330     size_t j;
1331     bool usePath;
1332     std::string pathName;
1333     std::string fName;
1334     usePath = MET_GetFilePath(m_FileName, pathName);
1335 
1336     if("Local" == m_ElementDataFileName ||
1337        "LOCAL" == m_ElementDataFileName ||
1338        "local" == m_ElementDataFileName)
1339       {
1340       M_ReadElements(_stream, m_ElementData, m_Quantity);
1341       }
1342     else if("LIST" == m_ElementDataFileName.substr(0, 4))
1343       {
1344       int fileImageDim = m_NDims - 1;
1345       int nWrds;
1346       char **wrds;
1347       MET_StringToWordArray(m_ElementDataFileName.c_str(), &nWrds, &wrds);
1348       if(nWrds > 1)
1349         {
1350         fileImageDim = (int)atof(wrds[1]);
1351         }
1352       for(i=0; i<nWrds; i++)
1353         {
1354         delete [] wrds[i];
1355         }
1356       delete [] wrds;
1357       if ( (fileImageDim == 0) || (fileImageDim > m_NDims) )
1358         {
1359         // if optional file dimension size is not given or is larger than
1360         // overall dimension then default to a size of m_NDims - 1.
1361         fileImageDim = m_NDims-1;
1362         }
1363       std::string s;
1364       std::ifstream* readStreamTemp = new std::ifstream;
1365       int elementSize;
1366       MET_SizeOfType(m_ElementType, &elementSize);
1367       elementSize *= m_ElementNumberOfChannels;
1368       int totalFiles = 1;
1369       for(i = m_NDims; i > fileImageDim; i--)
1370         {
1371         totalFiles *= m_DimSize[i-1];
1372         }
1373       for(i=0; i< totalFiles && !_stream->eof(); i++)
1374         {
1375         std::getline(*_stream, s);
1376         if(!_stream->eof())
1377           {
1378           j = s.length() - 1;
1379           while(j>0 && (isspace(s[j]) || !isprint(s[j])))
1380             {
1381             s[j--] = '\0';
1382             }
1383           if(usePath && !FileIsFullPath(s.c_str()))
1384             {
1385             fName = pathName + s;
1386             }
1387           else
1388             {
1389             fName = s;
1390             }
1391 
1392           openReadStream(*readStreamTemp, fName);
1393           if(!readStreamTemp->is_open())
1394             {
1395             std::cerr << "MetaImage: Read: cannot open slice"
1396                                 << std::endl;
1397             continue;
1398             }
1399           M_ReadElements(readStreamTemp,
1400                        &(((char *)m_ElementData)[i*m_SubQuantity[fileImageDim]*
1401                                                  elementSize]),
1402                        m_SubQuantity[fileImageDim]);
1403           readStreamTemp->close();
1404           }
1405         }
1406       delete readStreamTemp;
1407       }
1408     else if(m_ElementDataFileName.find('%') != std::string::npos)
1409       {
1410       int elementSize;
1411       MET_SizeOfType(m_ElementType, &elementSize);
1412       elementSize *= m_ElementNumberOfChannels;
1413 
1414       int nWrds;
1415       char **wrds;
1416       int minV = 1;
1417       int maxV = m_DimSize[m_NDims-1];
1418       int stepV = 1;
1419       std::string s;
1420       std::ifstream* readStreamTemp = new std::ifstream;
1421       MET_StringToWordArray(m_ElementDataFileName.c_str(), &nWrds, &wrds);
1422       if(nWrds >= 2)
1423         {
1424         minV = (int)atof(wrds[1]);
1425         maxV = minV + m_DimSize[m_NDims-1] - 1;
1426         }
1427       if(nWrds >= 3)
1428         {
1429         maxV = (int)atof(wrds[2]);
1430         stepV = (maxV-minV)/(m_DimSize[m_NDims-1]);
1431         }
1432       if(nWrds >= 4)
1433         {
1434         stepV = (int)atof(wrds[3]);
1435         }
1436       if(nWrds >= 5 )
1437       {
1438         // In this case, the filename must have had spaces in the
1439         // name.  The filename was parsed into multiple pieces by the
1440         // MET_StringToWordArray, which parses based on spaces.
1441         // Thus, we need to reconstruct the filename in this case.
1442         // The last three wrds must be numbers.  If they are not, we give an error.
1443         for( i = nWrds-3; i < nWrds; i++ )
1444         {
1445           for( j = 0; j < strlen(wrds[i]); j++ )
1446           {
1447             if( !isdigit(wrds[i][j]) )
1448             {
1449               std::cerr << "MetaImage: Read: Last three arguments must be numbers!"
1450                   << std::endl;
1451               continue;
1452             }
1453           }
1454         }
1455         stepV = (int)atof(wrds[nWrds-1]);
1456         maxV =  (int)atof(wrds[nWrds-2]);
1457         minV =  (int)atof(wrds[nWrds-3]);
1458         for( i = 1; i < nWrds-3; i++ )
1459         {
1460           strcat(wrds[0]," ");
1461           strcat(wrds[0],wrds[i]);
1462         }
1463       }
1464       // If the specified size of the third dimension is less than the size
1465       // specified by the regular expression, we should only read a volume with the specified
1466       // size.  Otherwise, the code will crash when trying to fill m_ElementData more than it can hold.
1467       // Therefore, we modify maxV to ensure that the images spanned by minV:stepV:maxV are less than or equal
1468       // to the size in the last dimension.
1469       int numberOfImages = 1 + (maxV - minV)/stepV;
1470       if( numberOfImages > m_DimSize[m_NDims-1] )
1471       {
1472         maxV = (m_DimSize[m_NDims-1]-1)*stepV + minV;
1473       }
1474       int cnt = 0;
1475       for(i=minV; i<=maxV; i += stepV)
1476         {
1477         s = string_format(wrds[0], i);
1478         if(usePath && !FileIsFullPath(s.c_str()))
1479           {
1480           fName = pathName + s;
1481           }
1482         else
1483           {
1484           fName = s;
1485           }
1486         openReadStream(*readStreamTemp,fName);
1487         if(!readStreamTemp->is_open())
1488           {
1489           std::cerr << "MetaImage: Read: cannot construct file"
1490                               << std::endl;
1491           continue;
1492           }
1493 
1494         M_ReadElements(readStreamTemp,
1495                        &(((char *)m_ElementData)[cnt*m_SubQuantity[m_NDims-1]*
1496                                                  elementSize]),
1497                        m_SubQuantity[m_NDims-1]);
1498         cnt++;
1499 
1500         readStreamTemp->close();
1501         }
1502       delete readStreamTemp;
1503       for(i=0; i<nWrds; i++)
1504         {
1505         delete [] wrds[i];
1506         }
1507       delete [] wrds;
1508       }
1509     else
1510       {
1511       if(usePath && !FileIsFullPath(m_ElementDataFileName.c_str()))
1512         {
1513         fName = pathName + m_ElementDataFileName;
1514         }
1515       else
1516         {
1517         fName = m_ElementDataFileName;
1518         }
1519 
1520       std::ifstream* readStreamTemp = new std::ifstream;
1521 
1522       const char *extensions[] = { "", ".gz", ".Z", nullptr };
1523       for(unsigned ii = 0; extensions[ii] != nullptr; ii++)
1524         {
1525         std::string tempFName(fName);
1526         tempFName += extensions[ii];
1527         openReadStream(*readStreamTemp,tempFName.c_str());
1528         if(readStreamTemp->is_open())
1529           {
1530           if(ii > 0)
1531             {
1532             this->CompressedData(true);
1533             this->BinaryData(true);
1534             }
1535           break;
1536           }
1537         }
1538 
1539       if(!readStreamTemp->is_open())
1540         {
1541         std::cerr << "MetaImage: Read: Cannot open data file"
1542                             << std::endl;
1543         if(m_ReadStream)
1544           {
1545           m_ReadStream->close();
1546           }
1547         return false;
1548         }
1549       M_ReadElements(readStreamTemp, m_ElementData, m_Quantity);
1550 
1551       readStreamTemp->close();
1552       delete readStreamTemp;
1553       }
1554     }
1555 
1556   return true;
1557 }
1558 
1559 
1560 bool MetaImage::
Write(const char * _headName,const char * _dataName,bool _writeElements,const void * _constElementData,bool _append)1561 Write(const char *_headName,
1562       const char *_dataName,
1563       bool _writeElements,
1564       const void * _constElementData,
1565       bool _append)
1566 {
1567   if(_headName != nullptr)
1568     {
1569     FileName(_headName);
1570     }
1571 
1572   bool userDataFileName = true;
1573   if(_dataName == nullptr && m_ElementDataFileName.empty())
1574     {
1575     userDataFileName = false;
1576     int sPtr = 0;
1577     MET_GetFileSuffixPtr(m_FileName, &sPtr);
1578     if(!strcmp(&m_FileName[sPtr], "mha"))
1579       {
1580       ElementDataFileName("LOCAL");
1581       }
1582     else
1583       {
1584       if(!_append)
1585         {
1586         MET_SetFileSuffix(m_FileName, "mhd");
1587         }
1588       m_ElementDataFileName = m_FileName;
1589       if(m_CompressedData)
1590         {
1591         MET_SetFileSuffix(m_ElementDataFileName, "zraw");
1592         }
1593       else
1594         {
1595         MET_SetFileSuffix(m_ElementDataFileName, "raw");
1596         }
1597       }
1598     }
1599   else if(_dataName != nullptr)
1600     {
1601     userDataFileName = false;
1602     ElementDataFileName(_dataName);
1603     }
1604 
1605   // make sure suffix is valid
1606   if(!_append)
1607     {
1608     if (m_ElementDataFileName == "LOCAL")
1609       {
1610       MET_SetFileSuffix(m_FileName, "mha");
1611       }
1612     else
1613       {
1614       MET_SetFileSuffix(m_FileName, "mhd");
1615       }
1616     }
1617 
1618   std::string pathName;
1619   bool usePath = MET_GetFilePath(m_FileName, pathName);
1620   if(usePath)
1621     {
1622     std::string elementPathName;
1623     MET_GetFilePath(m_ElementDataFileName, elementPathName);
1624     if (pathName == elementPathName)
1625       {
1626       elementPathName = m_ElementDataFileName.substr(pathName.length());
1627       m_ElementDataFileName = elementPathName;
1628       }
1629     }
1630 
1631   std::ofstream * tmpWriteStream = new std::ofstream;
1632 
1633   openWriteStream(*tmpWriteStream, m_FileName, _append);
1634 
1635   if(!tmpWriteStream->is_open())
1636     {
1637     if(!userDataFileName)
1638       {
1639       ElementDataFileName("");
1640       }
1641 
1642     delete tmpWriteStream;
1643 
1644     return false;
1645     }
1646 
1647   bool result = MetaImage::WriteStream(tmpWriteStream,
1648                                        _writeElements,
1649                                        _constElementData);
1650 
1651   if(!userDataFileName)
1652     {
1653     ElementDataFileName("");
1654     }
1655 
1656   tmpWriteStream->close();
1657   delete tmpWriteStream;
1658 
1659   return result;
1660 }
1661 
1662 bool MetaImage::
WriteStream(std::ofstream * _stream,bool _writeElements,const void * _constElementData)1663 WriteStream(std::ofstream * _stream,
1664             bool _writeElements,
1665             const void * _constElementData)
1666 {
1667   if(m_WriteStream != nullptr)
1668     {
1669     std::cerr << "MetaArray: WriteStream: two files open?"
1670                         << std::endl;
1671     delete m_WriteStream;
1672     }
1673 
1674   m_WriteStream = _stream;
1675 
1676   unsigned char * compressedElementData = nullptr;
1677   if(m_BinaryData && m_CompressedData && m_ElementDataFileName.find('%') == std::string::npos)
1678     // compressed & !slice/file
1679     {
1680     int elementSize;
1681     MET_SizeOfType(m_ElementType, &elementSize);
1682     int elementNumberOfBytes = elementSize*m_ElementNumberOfChannels;
1683 
1684     if(_constElementData == nullptr)
1685       {
1686       compressedElementData = MET_PerformCompression(
1687                                   (const unsigned char *)m_ElementData,
1688                                   m_Quantity * elementNumberOfBytes,
1689                                   & m_CompressedDataSize,
1690                                   m_CompressionLevel );
1691       }
1692     else
1693       {
1694       compressedElementData = MET_PerformCompression(
1695                                   (const unsigned char *)_constElementData,
1696                                   m_Quantity * elementNumberOfBytes,
1697                                   & m_CompressedDataSize,
1698                                   m_CompressionLevel );
1699       }
1700     }
1701 
1702   M_SetupWriteFields();
1703 
1704   M_Write();
1705 
1706   if(_writeElements)
1707     {
1708     if(m_BinaryData && m_CompressedData && m_ElementDataFileName.find('%') == std::string::npos)
1709       // compressed & !slice/file
1710       {
1711       M_WriteElements(m_WriteStream,
1712                       compressedElementData,
1713                       m_CompressedDataSize);
1714 
1715       delete [] compressedElementData;
1716       m_CompressedDataSize = 0;
1717       }
1718     else
1719       {
1720       if(_constElementData == nullptr)
1721         {
1722         M_WriteElements(m_WriteStream,
1723                         m_ElementData,
1724                         m_Quantity);
1725         }
1726       else
1727         {
1728         M_WriteElements(m_WriteStream,
1729                         _constElementData,
1730                         m_Quantity);
1731         }
1732       }
1733     }
1734 
1735   m_WriteStream = nullptr;
1736 
1737   return true;
1738 }
1739 
1740 
1741 /** Write a portion of an image */
WriteROI(int * _indexMin,int * _indexMax,const char * _headName,const char * _dataName,bool _writeElements,const void * _constElementData,bool _append)1742 bool MetaImage::WriteROI( int * _indexMin, int * _indexMax,
1743                           const char *_headName,
1744                           const char *_dataName,
1745                           bool _writeElements,
1746                           const void * _constElementData,
1747                           bool _append
1748                           )
1749 {
1750   if( _headName != nullptr )
1751     {
1752     FileName( _headName );
1753     }
1754 
1755   if( !_writeElements )
1756     {
1757     return false;
1758     }
1759 
1760   // Check if the file exists
1761   if( M_FileExists(_headName) )
1762     {
1763     char* elementData = const_cast<char*>(
1764                             static_cast<const char*>(_constElementData) );
1765     if( elementData == nullptr )
1766       {
1767       elementData = (char*)m_ElementData;
1768       }
1769     if( elementData == nullptr )
1770       {
1771       std::cerr << "Element data is NULL" << std::endl;
1772       return false;
1773       }
1774 
1775     // Find the start of the data
1776     std::ifstream * readStream = new std::ifstream;
1777     readStream->open( m_FileName, std::ios::binary |
1778                                   std::ios::in);
1779 
1780     // File must be readable
1781     if( !MetaObject::ReadStream( m_NDims, readStream ) )
1782       {
1783       std::cerr << "MetaImage: Read: Cannot parse file"
1784                           << std::endl;
1785       delete readStream;
1786       return false;
1787       }
1788 
1789     // File must not be compressed
1790     if(m_CompressedData)
1791       {
1792       std::cerr
1793                << "MetaImage cannot insert ROI into a compressed file."
1794                << std::endl;
1795       readStream->close();
1796       delete readStream;
1797       return false;
1798       }
1799 
1800     InitializeEssential( m_NDims,
1801                          m_DimSize,
1802                          m_ElementSpacing,
1803                          m_ElementType,
1804                          m_ElementNumberOfChannels,
1805                          nullptr, false ); // no memory allocation
1806 
1807     std::string  filename = ElementDataFileName();
1808     std::streampos dataPos = 0;
1809 
1810     // local file
1811     if( filename == "LOCAL" )
1812       {
1813       filename = m_FileName;
1814       dataPos = readStream->tellg();
1815       }
1816     else if( filename == "LIST"
1817              || strstr(filename.c_str(), "%") )
1818       {
1819       std::cerr
1820                << "MetaImage cannot insert ROI into a list of files."
1821                << std::endl;
1822       readStream->close();
1823       delete readStream;
1824       return false;
1825       }
1826 
1827     readStream->close();
1828     delete readStream;
1829 
1830     // Write the region
1831     if( !M_FileExists(filename.c_str()) )
1832       {
1833       std::string pathName;
1834       MET_GetFilePath(_headName, pathName);
1835       filename = pathName+filename;
1836       }
1837 
1838     std::ofstream * tmpWriteStream = new std::ofstream;
1839     tmpWriteStream->open( filename.c_str(),
1840                           std::ios::binary |
1841                           std::ios::in |
1842                           std::ios::out );
1843 
1844     if( !tmpWriteStream->is_open() )
1845       {
1846       std::cerr << "Cannot open ROI file: "
1847                           << filename.c_str()
1848                           << std::endl;
1849       delete tmpWriteStream;
1850       return false;
1851       }
1852 
1853     int elementSize;
1854     MET_SizeOfType( m_ElementType, &elementSize );
1855     int elementNumberOfBytes = elementSize*m_ElementNumberOfChannels;
1856 
1857     // seek to the end and write one byte to allocate the entire file size
1858     std::streamoff seekoff = m_Quantity*elementNumberOfBytes;
1859     tmpWriteStream->seekp(0, std::ios::end);
1860     if (tmpWriteStream->tellp() != (dataPos+seekoff))
1861       {
1862       seekoff = seekoff - 1;
1863       tmpWriteStream->seekp(dataPos+seekoff, std::ios::beg);
1864       const char zerobyte = 0;
1865       tmpWriteStream->write(&zerobyte, 1);
1866       }
1867 
1868     if( !elementData )
1869       {
1870       std::cerr << "Element data is NULL" << std::endl;
1871       delete tmpWriteStream;
1872       return false;
1873       }
1874 
1875     M_WriteElementsROI(tmpWriteStream, elementData, dataPos,
1876                        _indexMin, _indexMax);
1877 
1878     tmpWriteStream->close();
1879     delete tmpWriteStream;
1880     }
1881   else // the file doesn't exist
1882     {
1883     if(m_CompressedData)
1884       {
1885       std::cerr
1886                << "MetaImage cannot write an ROI using compression."
1887                << std::endl;
1888       return false;
1889       }
1890 
1891     // Get the data filename right...
1892     bool userDataFileName = true;
1893     if( _dataName == nullptr && m_ElementDataFileName.empty() )
1894       {
1895       userDataFileName = false;
1896       int sPtr = 0;
1897       MET_GetFileSuffixPtr(m_FileName, &sPtr);
1898       if( !strcmp(&m_FileName[sPtr], "mha") )
1899         {
1900         ElementDataFileName( "LOCAL" );
1901         }
1902       else
1903         {
1904         if(!_append)
1905           {
1906           MET_SetFileSuffix(m_FileName, "mhd");
1907           }
1908         m_ElementDataFileName = m_FileName;
1909         if(m_CompressedData)
1910           {
1911           MET_SetFileSuffix(m_ElementDataFileName, "zraw");
1912           }
1913         else
1914           {
1915           MET_SetFileSuffix(m_ElementDataFileName, "raw");
1916           }
1917         }
1918       }
1919     else if(_dataName != nullptr)
1920       {
1921       userDataFileName = false;
1922       ElementDataFileName(_dataName);
1923       }
1924 
1925     if( m_ElementDataFileName == "LIST"
1926         || m_ElementDataFileName.find('%') != std::string::npos)
1927       {
1928       std::cerr
1929                << "MetaImage cannot insert ROI into a list of files."
1930                << std::endl;
1931       return false;
1932       }
1933 
1934     // make sure the header suffix is valid, unless forcing to match an
1935     // existing file via the append bool argument.
1936     if(!_append)
1937       {
1938       if (m_ElementDataFileName == "LOCAL")
1939         {
1940         MET_SetFileSuffix(m_FileName, "mha");
1941         }
1942       else
1943         {
1944         MET_SetFileSuffix(m_FileName, "mhd");
1945         }
1946       }
1947 
1948     std::string pathName;
1949     bool usePath = MET_GetFilePath(m_FileName, pathName);
1950     if(usePath)
1951       {
1952       std::string elementPathName;
1953       MET_GetFilePath(m_ElementDataFileName, elementPathName);
1954       if (pathName == elementPathName)
1955         {
1956         elementPathName = m_ElementDataFileName.substr(pathName.length());
1957         m_ElementDataFileName = elementPathName;
1958         }
1959       }
1960 
1961     std::ofstream * tmpWriteStream = new std::ofstream;
1962 
1963     openWriteStream(*tmpWriteStream, m_FileName, _append);
1964 
1965     if(!tmpWriteStream->is_open())
1966       {
1967       if(!userDataFileName)
1968         {
1969         ElementDataFileName("");
1970         }
1971       delete tmpWriteStream;
1972       return false;
1973       }
1974 
1975     // Write the ROI header file
1976     char* elementData = const_cast<char*>(
1977                           static_cast<const char*>(_constElementData) );
1978     if(elementData == nullptr)
1979       {
1980       elementData = (char*)m_ElementData;
1981       }
1982 
1983     m_WriteStream = tmpWriteStream;
1984     M_SetupWriteFields();
1985     M_Write();
1986 
1987     std::streampos dataPos = m_WriteStream->tellp();
1988 
1989     // If data is in a separate file, set dataPos and point to that file.
1990     //   ( we've already verified the name isn't LIST and doesn't
1991     //     contain % )
1992     if( m_ElementDataFileName != "LOCAL" )
1993       {
1994       m_WriteStream = nullptr;
1995       tmpWriteStream->close();
1996 
1997       dataPos = 0;
1998 
1999       std::string dataFileName;
2000       if(usePath&& !FileIsFullPath(m_ElementDataFileName.c_str()))
2001         {
2002         dataFileName = pathName + m_ElementDataFileName;
2003         }
2004       else
2005         {
2006         dataFileName = m_ElementDataFileName;
2007         }
2008 
2009       openWriteStream(*tmpWriteStream, dataFileName, _append);
2010       m_WriteStream = tmpWriteStream;
2011       }
2012 
2013     int elementSize;
2014     MET_SizeOfType( m_ElementType, &elementSize );
2015     int elementNumberOfBytes = elementSize*m_ElementNumberOfChannels;
2016 
2017     // write the last byte in the file to allocate it
2018     std::streamoff seekoff = m_Quantity * elementNumberOfBytes;
2019     seekoff -= 1;
2020     m_WriteStream->seekp(seekoff, std::ios::cur);
2021     const char zerobyte = 0;
2022     m_WriteStream->write(&zerobyte, 1);
2023 
2024     M_WriteElementsROI(m_WriteStream, elementData, dataPos,
2025                        _indexMin, _indexMax);
2026 
2027     m_WriteStream = nullptr;
2028 
2029     if(!userDataFileName)
2030       {
2031       ElementDataFileName("");
2032       }
2033 
2034     tmpWriteStream->close();
2035     delete tmpWriteStream;
2036     }
2037 
2038   return true;
2039 }
2040 
2041 bool  MetaImage::
M_WriteElementsROI(std::ofstream * _fstream,const void * _data,std::streampos _dataPos,int * _indexMin,int * _indexMax)2042 M_WriteElementsROI(std::ofstream * _fstream,
2043                    const void * _data,
2044                    std::streampos _dataPos,
2045                    int * _indexMin,
2046                    int* _indexMax )
2047 {
2048   const char* data = static_cast<const char*>(_data);
2049 
2050   int elementSize;
2051   MET_SizeOfType(m_ElementType, &elementSize);
2052   const int elementNumberOfBytes = elementSize*m_ElementNumberOfChannels;
2053 
2054   // Write the IO region line by line
2055   int * currentIndex = new int[m_NDims];
2056   for(int i=0; i<m_NDims; i++)
2057     {
2058     currentIndex[i] = _indexMin[i];
2059     }
2060 
2061   // Optimize the size of the buffer to written depending on the
2062   // region shape
2063   // This calculate the number of continuous bytes in the file
2064   // which can be written
2065   std::streamoff elementsToWrite = 1;
2066   int movingDirection = 0;
2067   do
2068     {
2069     elementsToWrite *= _indexMax[movingDirection] - _indexMin[movingDirection] + 1;
2070     ++movingDirection;
2071     }
2072   while(movingDirection < m_NDims
2073         && _indexMin[movingDirection-1] == 0
2074         && _indexMax[movingDirection-1] == m_DimSize[movingDirection-1]-1);
2075 
2076   // write line by line
2077   bool done = false;
2078   while(!done)
2079     {
2080     // Seek to the right position
2081     std::streamoff seekoff = _dataPos;
2082     for(int i=0; i<m_NDims; i++)
2083       {
2084       seekoff += m_SubQuantity[i] * currentIndex[i] * elementNumberOfBytes;
2085       }
2086     _fstream->seekp( seekoff, std::ios::beg );
2087 
2088     M_WriteElementData( _fstream, data, elementsToWrite );
2089     data += elementsToWrite * elementNumberOfBytes;
2090 
2091     // check if there is only one write needed
2092     if( movingDirection >= m_NDims )
2093       {
2094       break;
2095       }
2096 
2097     ++currentIndex[movingDirection];
2098 
2099     // Check if we are still in the region
2100     for( int j=movingDirection; j<m_NDims; j++ )
2101       {
2102       if( currentIndex[j] > _indexMax[j] )
2103         {
2104         if( j == m_NDims-1 )
2105           {
2106           done = true;
2107           break;
2108           }
2109         else
2110           {
2111           currentIndex[j] = _indexMin[j];
2112           currentIndex[j+1]++;
2113           }
2114         }
2115       }
2116     } // end writing  loop
2117 
2118   delete [] currentIndex;
2119 
2120   return true;
2121 }
2122 
2123 bool MetaImage::
Append(const char * _headName)2124 Append(const char *_headName)
2125 {
2126   if(META_DEBUG)
2127    {
2128    std::cout << "MetaImage: Append" << std::endl;
2129    }
2130 
2131   return this->Write(_headName, nullptr, true, nullptr, true);
2132 }
2133 
2134 void MetaImage::
M_Destroy()2135 M_Destroy()
2136 {
2137   if(m_AutoFreeElementData && m_ElementData != nullptr)
2138     {
2139     delete [] (char *)m_ElementData;
2140     }
2141 
2142   m_ElementData = nullptr;
2143 
2144   if(m_CompressionTable && m_CompressionTable->compressedStream)
2145     {
2146     inflateEnd(m_CompressionTable->compressedStream);
2147     delete m_CompressionTable->compressedStream;
2148     delete [] m_CompressionTable->buffer;
2149     m_CompressionTable->buffer = nullptr;
2150     }
2151   delete m_CompressionTable;
2152   m_CompressionTable = nullptr;
2153 
2154   MetaObject::M_Destroy();
2155 }
2156 
2157 void MetaImage::
M_SetupReadFields()2158 M_SetupReadFields()
2159 {
2160   if(META_DEBUG)
2161     {
2162     std::cout << "MetaImage: M_SetupReadFields"
2163                         << std::endl;
2164     }
2165 
2166   MetaObject::M_SetupReadFields();
2167 
2168   MET_FieldRecordType * mF;
2169 
2170   int nDimsRecNum = MET_GetFieldRecordNumber("NDims", &m_Fields);
2171 
2172   mF = new MET_FieldRecordType;
2173   MET_InitReadField(mF, "DimSize", MET_INT_ARRAY, true, nDimsRecNum);
2174   mF->required = true;
2175   m_Fields.push_back(mF);
2176 
2177   mF = new MET_FieldRecordType;
2178   MET_InitReadField(mF, "HeaderSize", MET_INT, false);
2179   m_Fields.push_back(mF);
2180 
2181   mF = new MET_FieldRecordType;
2182   MET_InitReadField(mF, "Modality", MET_STRING, false);
2183   m_Fields.push_back(mF);
2184 
2185   mF = new MET_FieldRecordType;
2186   MET_InitReadField(mF, "ImagePosition", MET_FLOAT_ARRAY, false, nDimsRecNum);
2187   m_Fields.push_back(mF);
2188 
2189   mF = new MET_FieldRecordType;
2190   MET_InitReadField(mF, "SequenceID", MET_INT_ARRAY, false, nDimsRecNum);
2191   m_Fields.push_back(mF);
2192 
2193   mF = new MET_FieldRecordType;
2194   MET_InitReadField(mF, "ElementMin", MET_FLOAT, false);
2195   m_Fields.push_back(mF);
2196 
2197   mF = new MET_FieldRecordType;
2198   MET_InitReadField(mF, "ElementMax", MET_FLOAT, false);
2199   m_Fields.push_back(mF);
2200 
2201   mF = new MET_FieldRecordType;
2202   MET_InitReadField(mF, "ElementNumberOfChannels", MET_INT, false);
2203   m_Fields.push_back(mF);
2204 
2205   mF = new MET_FieldRecordType;
2206   MET_InitReadField(mF, "ElementSize", MET_FLOAT_ARRAY, false, nDimsRecNum);
2207   m_Fields.push_back(mF);
2208 
2209   mF = new MET_FieldRecordType;  // Set but not used...
2210   MET_InitReadField(mF, "ElementNBits", MET_INT, false);
2211   m_Fields.push_back(mF);
2212 
2213   mF = new MET_FieldRecordType;  // Used by ConvertElementToIntensity funcs
2214   MET_InitReadField(mF, "ElementToIntensityFunctionSlope", MET_FLOAT, false);
2215   m_Fields.push_back(mF);
2216 
2217   mF = new MET_FieldRecordType;  // Used by ConvertElementToIntensity funcs
2218   MET_InitReadField(mF, "ElementToIntensityFunctionOffset", MET_FLOAT, false);
2219   m_Fields.push_back(mF);
2220 
2221   mF = new MET_FieldRecordType;
2222   MET_InitReadField(mF, "ElementType", MET_STRING, true);
2223   mF->required = true;
2224   m_Fields.push_back(mF);
2225 
2226   mF = new MET_FieldRecordType;
2227   MET_InitReadField(mF, "ElementDataFile", MET_STRING, true);
2228   mF->required = true;
2229   mF->terminateRead = true;
2230   m_Fields.push_back(mF);
2231 }
2232 
2233 void MetaImage::
M_SetupWriteFields()2234 M_SetupWriteFields()
2235 {
2236   MetaObject::M_SetupWriteFields();
2237 
2238   MET_FieldRecordType * mF;
2239 
2240   mF = new MET_FieldRecordType;
2241   MET_InitWriteField(mF, "DimSize", MET_INT_ARRAY, m_NDims, m_DimSize);
2242   m_Fields.push_back(mF);
2243 
2244   char s[22];
2245   if(m_HeaderSize > 0 || m_HeaderSize == -1)
2246     {
2247     mF = new MET_FieldRecordType;
2248     MET_InitWriteField(mF, "HeaderSize", MET_INT);
2249     m_Fields.push_back(mF);
2250     }
2251 
2252   int i;
2253   if(m_Modality != MET_MOD_UNKNOWN)
2254     {
2255     mF = new MET_FieldRecordType;
2256     strcpy(s, MET_ValueTypeName[m_Modality]);
2257     MET_InitWriteField(mF, "Modality", MET_STRING, strlen(s), s);
2258     m_Fields.push_back(mF);
2259     }
2260 
2261   i = MET_GetFieldRecordNumber("AnatomicalOrientation", &m_Fields);
2262   if(i < 0)
2263     {
2264     const char * str = AnatomicalOrientationAcronym();
2265     mF = new MET_FieldRecordType;
2266     MET_InitWriteField(mF, "AnatomicalOrientation",
2267                        MET_STRING, strlen(str), str);
2268     m_Fields.push_back(mF);
2269     }
2270 
2271   bool valid = false;
2272   for(i=0; i<4; i++)
2273     {
2274     if(m_SequenceID[i] != 0)
2275       {
2276       valid = true;
2277       break;
2278       }
2279     }
2280   if(valid)
2281     {
2282     mF = new MET_FieldRecordType;
2283     MET_InitWriteField(mF, "SequenceID", MET_FLOAT_ARRAY, m_NDims,
2284                        m_SequenceID);
2285     m_Fields.push_back(mF);
2286     }
2287 
2288   if(m_ElementMinMaxValid)
2289     {
2290     mF = new MET_FieldRecordType;
2291     MET_InitWriteField(mF, "ElementMin", MET_FLOAT, m_ElementMin);
2292     m_Fields.push_back(mF);
2293 
2294     mF = new MET_FieldRecordType;
2295     MET_InitWriteField(mF, "ElementMax", MET_FLOAT, m_ElementMax);
2296     m_Fields.push_back(mF);
2297     }
2298 
2299   if(m_ElementNumberOfChannels>1)
2300     {
2301     mF = new MET_FieldRecordType;
2302     MET_InitWriteField(mF, "ElementNumberOfChannels", MET_INT,
2303                        m_ElementNumberOfChannels);
2304     m_Fields.push_back(mF);
2305     }
2306 
2307   if(m_ElementSizeValid)
2308     {
2309     mF = new MET_FieldRecordType;
2310     MET_InitWriteField(mF, "ElementSize", MET_FLOAT_ARRAY, m_NDims,
2311                        m_ElementSize);
2312     m_Fields.push_back(mF);
2313     }
2314 
2315   if(m_ElementToIntensityFunctionSlope != 1 ||
2316      m_ElementToIntensityFunctionOffset != 0)
2317     {
2318     mF = new MET_FieldRecordType;
2319     MET_InitWriteField(mF, "ElementToIntensityFunctionSlope",
2320                        MET_FLOAT, m_ElementToIntensityFunctionSlope);
2321     m_Fields.push_back(mF);
2322     mF = new MET_FieldRecordType;
2323     MET_InitWriteField(mF, "ElementToIntensityFunctionOffset",
2324                        MET_FLOAT, m_ElementToIntensityFunctionOffset);
2325     m_Fields.push_back(mF);
2326     }
2327 
2328   mF = new MET_FieldRecordType;
2329   MET_TypeToString(m_ElementType, s);
2330   MET_InitWriteField(mF, "ElementType", MET_STRING, strlen(s), s);
2331   m_Fields.push_back(mF);
2332 
2333   mF = new MET_FieldRecordType;
2334   MET_InitWriteField(mF, "ElementDataFile", MET_STRING,
2335                      m_ElementDataFileName.length(),
2336                      m_ElementDataFileName.c_str());
2337   mF->terminateRead = true;
2338   m_Fields.push_back(mF);
2339 }
2340 
2341 bool MetaImage::
M_Read()2342 M_Read()
2343 {
2344   if(META_DEBUG)
2345     {
2346     std::cout << "MetaImage: M_Read: Loading Header"
2347                         << std::endl;
2348     }
2349   if(!MetaObject::M_Read())
2350     {
2351     std::cerr << "MetaImage: M_Read: Error parsing file"
2352                         << std::endl;
2353     return false;
2354     }
2355 
2356   if(META_DEBUG)
2357     {
2358     std::cout << "MetaImage: M_Read: Parsing Header"
2359                         << std::endl;
2360     }
2361   MET_FieldRecordType * mF;
2362 
2363   if(META_DEBUG)
2364     {
2365     std::cout << "metaImage: M_Read: elementSpacing[" << 0 << "] = "
2366                         << m_ElementSpacing[0] << std::endl;
2367     }
2368   mF = MET_GetFieldRecord("DimSize", &m_Fields);
2369   if(mF && mF->defined)
2370     {
2371     int i;
2372     for(i=0; i<m_NDims; i++)
2373       {
2374       m_DimSize[i] = (int)mF->value[i];
2375       }
2376     }
2377 
2378   mF = MET_GetFieldRecord("HeaderSize", &m_Fields);
2379   if(mF && mF->defined)
2380     {
2381     m_HeaderSize = (int)mF->value[0];
2382     }
2383 
2384   mF = MET_GetFieldRecord("Modality", &m_Fields);
2385   if(mF && mF->defined)
2386     {
2387     MET_StringToImageModality((char *)mF->value, &m_Modality);
2388     }
2389 
2390   mF = MET_GetFieldRecord("SequenceID", &m_Fields);
2391   if(mF && mF->defined)
2392     {
2393     int i;
2394     for(i=0; i<m_NDims; i++)
2395       {
2396       m_SequenceID[i] = (float)(mF->value[i]);
2397       }
2398     }
2399 
2400   mF = MET_GetFieldRecord("ImagePosition", &m_Fields);
2401   if(mF && mF->defined)
2402     {
2403     int i;
2404     for(i=0; i<m_NDims; i++)
2405       {
2406       m_Offset[i] = static_cast<double>(mF->value[i]);
2407       }
2408     }
2409 
2410   mF = MET_GetFieldRecord("ElementMin", &m_Fields);
2411   if(mF && mF->defined)
2412     {
2413     m_ElementMin = mF->value[0];
2414     }
2415 
2416   mF = MET_GetFieldRecord("ElementMax", &m_Fields);
2417   if(mF && mF->defined)
2418     {
2419     m_ElementMax = mF->value[0];
2420     }
2421 
2422   mF = MET_GetFieldRecord("ElementNumberOfChannels", &m_Fields);
2423   if(mF && mF->defined)
2424     {
2425     m_ElementNumberOfChannels = (int)mF->value[0];
2426     }
2427 
2428 
2429   mF = MET_GetFieldRecord("ElementSize", &m_Fields);
2430   if(mF && mF->defined)
2431     {
2432     m_ElementSizeValid = true;
2433     int i;
2434     for(i=0; i<m_NDims; i++)
2435       {
2436       m_ElementSize[i] = mF->value[i];
2437       }
2438     mF = MET_GetFieldRecord("ElementSpacing", &m_Fields);
2439     if(mF && !mF->defined)
2440       {
2441       for(i=0; i<m_NDims; i++)
2442         {
2443         m_ElementSpacing[i] = m_ElementSize[i];
2444         }
2445       }
2446     }
2447   else
2448     {
2449     int i;
2450     m_ElementSizeValid = false;
2451     for(i=0; i<m_NDims; i++)
2452       {
2453       m_ElementSize[i] = m_ElementSpacing[i];
2454       }
2455     }
2456 
2457   m_ElementToIntensityFunctionSlope = 1;
2458   m_ElementToIntensityFunctionOffset = 0;
2459   mF = MET_GetFieldRecord("ElementToIntensityFunctionSlope", &m_Fields);
2460   if(mF && mF->defined)
2461     {
2462     m_ElementToIntensityFunctionSlope = mF->value[0];
2463     }
2464   mF = MET_GetFieldRecord("ElementToIntensityFunctionOffset", &m_Fields);
2465   if(mF && mF->defined)
2466     {
2467     m_ElementToIntensityFunctionOffset = mF->value[0];
2468     }
2469 
2470   mF = MET_GetFieldRecord("ElementType", &m_Fields);
2471   if(mF && mF->defined)
2472     {
2473     MET_StringToType((char *)(mF->value), &m_ElementType);
2474     }
2475 
2476   mF = MET_GetFieldRecord("ElementDataFile", &m_Fields);
2477   if(mF && mF->defined)
2478     {
2479     m_ElementDataFileName = (char *)(mF->value);
2480     }
2481 
2482   return true;
2483 }
2484 
2485 bool MetaImage::
M_ReadElements(std::ifstream * _fstream,void * _data,std::streamoff _dataQuantity)2486 M_ReadElements(std::ifstream * _fstream, void * _data,
2487                std::streamoff _dataQuantity)
2488 {
2489   if(META_DEBUG)
2490     {
2491     std::cout << "MetaImage: M_ReadElements" << std::endl;
2492     }
2493 
2494   if(m_HeaderSize>(int)0)
2495     {
2496     _fstream->seekg(m_HeaderSize, std::ios::beg);
2497     if(!_fstream->good())
2498       {
2499       std::cerr << "MetaImage: Read: header not read correctly"
2500                           << std::endl;
2501       return false;
2502       }
2503     }
2504 
2505   int elementSize;
2506   MET_SizeOfType(m_ElementType, &elementSize);
2507   std::streamoff readSize = _dataQuantity*m_ElementNumberOfChannels*elementSize;
2508   if(META_DEBUG)
2509     {
2510     std::cout << "MetaImage: M_ReadElements: ReadSize = "
2511                         << readSize << std::endl;
2512     }
2513 
2514   if(m_HeaderSize == -1)
2515     {
2516     if(META_DEBUG)
2517       {
2518       std::cout << "MetaImage: M_ReadElements: Skipping header"
2519                           << std::endl;
2520       }
2521     _fstream->seekg(-readSize, std::ios::end);
2522     }
2523 
2524   // If compressed we inflate
2525   if(m_BinaryData && m_CompressedData)
2526     {
2527     // if m_CompressedDataSize is not defined we assume the size of the
2528     // file is the size of the compressed data
2529     bool compressedDataDeterminedFromFile = false;
2530     if(m_CompressedDataSize==0)
2531       {
2532       compressedDataDeterminedFromFile = true;
2533       _fstream->seekg(0, std::ios::end);
2534       m_CompressedDataSize = _fstream->tellg();
2535       _fstream->seekg(0, std::ios::beg);
2536       }
2537 
2538     unsigned char* compr = new unsigned char[static_cast<size_t>(m_CompressedDataSize)];
2539 
2540     M_ReadElementData( _fstream, compr, m_CompressedDataSize );
2541 
2542     MET_PerformUncompression(compr, m_CompressedDataSize,
2543                              (unsigned char *)_data, readSize);
2544 
2545     if (compressedDataDeterminedFromFile)
2546       {
2547       m_CompressedDataSize = 0;
2548       }
2549 
2550     delete [] compr;
2551     }
2552   else // if not compressed
2553     {
2554     if(!m_BinaryData)
2555       {
2556 
2557       M_ReadElementData( _fstream, _data, _dataQuantity );
2558 
2559       }
2560     else
2561       {
2562 
2563       if ( !M_ReadElementData( _fstream, _data, _dataQuantity ) )
2564         return false;
2565 
2566       }
2567     }
2568 
2569   return true;
2570 }
2571 
2572 bool MetaImage::
M_WriteElements(std::ofstream * _fstream,const void * _data,std::streamoff _dataQuantity)2573 M_WriteElements(std::ofstream * _fstream,
2574                 const void * _data,
2575                 std::streamoff _dataQuantity)
2576 {
2577 
2578   if (m_ElementDataFileName == "LOCAL")
2579     {
2580     MetaImage::M_WriteElementData(_fstream, _data, _dataQuantity);
2581     }
2582   else // write the data in a separate file
2583     {
2584     std::string dataFileName;
2585     std::string pathName;
2586     bool usePath = MET_GetFilePath(m_FileName, pathName);
2587     if(usePath&& !FileIsFullPath(m_ElementDataFileName.c_str()))
2588       {
2589       dataFileName = pathName + m_ElementDataFileName;
2590       }
2591     else
2592       {
2593       dataFileName = m_ElementDataFileName;
2594       }
2595 
2596     if(dataFileName.find('%') != std::string::npos) // write slice by slice
2597       {
2598       int i;
2599       std::string fName;
2600       int elementSize;
2601       MET_SizeOfType(m_ElementType, &elementSize);
2602       std::streamoff elementNumberOfBytes = elementSize*m_ElementNumberOfChannels;
2603       std::streamoff sliceNumberOfBytes = m_SubQuantity[m_NDims-1]*elementNumberOfBytes;
2604 
2605       std::ofstream* writeStreamTemp = new std::ofstream;
2606       for(i=1; i<=m_DimSize[m_NDims-1]; i++)
2607         {
2608         fName = string_format(dataFileName, i);
2609 
2610         openWriteStream(*writeStreamTemp, fName, false);
2611 
2612         if(!m_CompressedData)
2613           {
2614           // BUG? This looks wrong to me as the third parameter should
2615           // contain the number of elements/quantity, not number of bytes -BCL
2616           MetaImage::M_WriteElementData(writeStreamTemp,
2617                              &(((const char *)_data)[(i-1)*sliceNumberOfBytes]),
2618                              sliceNumberOfBytes);
2619           }
2620         else
2621           {
2622           unsigned char * compressedData = nullptr;
2623           std::streamoff compressedDataSize = 0;
2624 
2625           // Compress the data slice by slice
2626           compressedData = MET_PerformCompression(
2627                   &(((const unsigned char *)_data)[(i-1)*sliceNumberOfBytes]),
2628                   sliceNumberOfBytes,
2629                   & compressedDataSize,
2630                   m_CompressionLevel );
2631 
2632           // Write the compressed data
2633           MetaImage::M_WriteElementData( writeStreamTemp,
2634                               compressedData,
2635                               compressedDataSize );
2636 
2637           delete [] compressedData;
2638           }
2639 
2640         writeStreamTemp->close();
2641         }
2642 
2643       delete writeStreamTemp;
2644       }
2645     else // write the image in one unique other file
2646       {
2647       std::ofstream* writeStreamTemp = new std::ofstream;
2648       openWriteStream(*writeStreamTemp, dataFileName, false);
2649 
2650       MetaImage::M_WriteElementData(writeStreamTemp, _data, _dataQuantity);
2651 
2652       writeStreamTemp->close();
2653       delete writeStreamTemp;
2654       }
2655     }
2656 
2657   return true;
2658 }
2659 
2660 
2661 bool MetaImage::
M_WriteElementData(std::ofstream * _fstream,const void * _data,std::streamoff _dataQuantity)2662 M_WriteElementData(std::ofstream * _fstream,
2663                    const void * _data,
2664                    std::streamoff _dataQuantity)
2665 {
2666   if(!m_BinaryData)
2667     {
2668 
2669     double tf;
2670     for(std::streamoff i=0; i<_dataQuantity; i++)
2671       {
2672       MET_ValueToDouble(m_ElementType, _data, i, &tf);
2673       if((i+1)/10 == (double)(i+1.0)/10.0)
2674         {
2675         (*_fstream) << tf << std::endl;
2676         }
2677       else
2678         {
2679         (*_fstream) << tf << " ";
2680         }
2681       }
2682     }
2683   else
2684     {
2685     if(m_CompressedData)
2686       {
2687       // the data is written in writes no bigger then MaxIOChunk
2688       std::streamoff bytesRemaining = _dataQuantity;
2689       while ( bytesRemaining )
2690         {
2691         std::streamoff chunkToWrite = bytesRemaining > MaxIOChunk ? MaxIOChunk : bytesRemaining;
2692         _fstream->write( (const char *)_data, (size_t)chunkToWrite );
2693         _data = (const char *)(_data) + chunkToWrite; // <- Note: data is changed
2694         bytesRemaining -= chunkToWrite;
2695         }
2696       }
2697     else
2698       {
2699       int elementSize;
2700       MET_SizeOfType(m_ElementType, &elementSize);
2701       std::streamoff elementNumberOfBytes = elementSize*m_ElementNumberOfChannels;
2702 
2703       // the data is written in writes no bigger then MaxIOChunk
2704       std::streamoff bytesRemaining = _dataQuantity * elementNumberOfBytes;
2705       while ( bytesRemaining )
2706         {
2707         std::streamoff chunkToWrite = bytesRemaining > MaxIOChunk ? MaxIOChunk : bytesRemaining;
2708         _fstream->write( (const char *)_data, (size_t)chunkToWrite );
2709         _data = (const char *)(_data) + chunkToWrite; // <- Note: _data is changed
2710         bytesRemaining -= chunkToWrite;
2711         }
2712       }
2713     }
2714 
2715   // check if the io stream did not fail in the process of writing
2716   if ( _fstream->fail() )
2717     {
2718     std::cerr
2719       << "MetaImage: M_WriteElementsData: file stream is fail after write"
2720       << std::endl;
2721     return false;
2722     }
2723 
2724   return true;
2725 }
2726 
2727 /** Streaming related functions */
2728 bool MetaImage::
ReadROI(int * _indexMin,int * _indexMax,const char * _headerName,bool _readElements,void * _buffer,unsigned int subSamplingFactor)2729 ReadROI(int * _indexMin, int * _indexMax,
2730         const char *_headerName,
2731         bool _readElements,
2732         void * _buffer,
2733         unsigned int subSamplingFactor)
2734 {
2735   M_Destroy();
2736 
2737   Clear();
2738 
2739   M_SetupReadFields();
2740 
2741   if(_headerName != nullptr)
2742     {
2743     m_FileName = _headerName;
2744     }
2745 
2746   M_PrepareNewReadStream();
2747 
2748   std::ifstream * tmpReadStream = new std::ifstream;
2749 
2750   openReadStream(*tmpReadStream, m_FileName);
2751 
2752   if(!tmpReadStream->is_open())
2753     {
2754     delete tmpReadStream;
2755     return false;
2756     }
2757 
2758   if( !this->ReadROIStream(_indexMin, _indexMax,
2759                            0, tmpReadStream, _readElements, _buffer,subSamplingFactor) )
2760     {
2761     tmpReadStream->close();
2762     delete tmpReadStream;
2763     return false;
2764     }
2765 
2766   tmpReadStream->close();
2767 
2768   delete tmpReadStream;
2769 
2770   return true;
2771 }
2772 
2773 /** Read the ROI Stream */
ReadROIStream(int * _indexMin,int * _indexMax,int _nDims,std::ifstream * _stream,bool _readElements,void * _buffer,unsigned int subSamplingFactor)2774 bool MetaImage::ReadROIStream(int * _indexMin, int * _indexMax,
2775                               int _nDims,
2776                               std::ifstream * _stream,
2777                               bool _readElements,
2778                               void * _buffer,
2779                               unsigned int subSamplingFactor)
2780 {
2781   if(!MetaObject::ReadStream(_nDims, _stream))
2782     {
2783     std::cerr << "MetaImage: Read: Cannot parse file"
2784                         << std::endl;
2785     return false;
2786     }
2787 
2788   if(_readElements)
2789     {
2790     if(_buffer == nullptr)
2791       {
2792       InitializeEssential(m_NDims,
2793                           m_DimSize,
2794                           m_ElementSpacing,
2795                           m_ElementType,
2796                           m_ElementNumberOfChannels,
2797                           nullptr, true);
2798       }
2799     else
2800       {
2801       InitializeEssential(m_NDims,
2802                           m_DimSize,
2803                           m_ElementSpacing,
2804                           m_ElementType,
2805                           m_ElementNumberOfChannels,
2806                           _buffer, false);
2807       }
2808 
2809     // Streaming related. We need to update some of the fields
2810     std::streamoff quantity = 1;
2811     int i;
2812     size_t j;
2813     for(i=0; i<m_NDims; i++)
2814       {
2815       quantity *= (_indexMax[i] - _indexMin[i] + 1);
2816       }
2817 
2818     bool usePath;
2819     std::string pathName;
2820     std::string fName;
2821     usePath = MET_GetFilePath(m_FileName, pathName);
2822 
2823     if("Local" == m_ElementDataFileName ||
2824        "LOCAL" == m_ElementDataFileName ||
2825        "local" == m_ElementDataFileName)
2826       {
2827       M_ReadElementsROI(_stream, m_ElementData, quantity,
2828                         _indexMin, _indexMax, subSamplingFactor,
2829                         m_Quantity);
2830       }
2831     else if("LIST" == m_ElementDataFileName.substr(0, 4))
2832       {
2833       int fileImageDim = m_NDims - 1;
2834       int nWrds;
2835       char **wrds;
2836       MET_StringToWordArray(m_ElementDataFileName.c_str(), &nWrds, &wrds);
2837       if(nWrds > 1)
2838         {
2839         fileImageDim = (int)atof(wrds[1]);
2840         }
2841       for(i=0; i<nWrds; i++)
2842         {
2843         delete [] wrds[i];
2844         }
2845       delete [] wrds;
2846       if ( (fileImageDim == 0) || (fileImageDim > m_NDims) )
2847         {
2848         // if optional file dimension size is not given or is larger than
2849         // overall dimension then default to a size of m_NDims - 1.
2850         fileImageDim = m_NDims-1;
2851         }
2852       char s[1024];
2853       std::ifstream* readStreamTemp = new std::ifstream;
2854       int elementSize;
2855       MET_SizeOfType(m_ElementType, &elementSize);
2856       elementSize *= m_ElementNumberOfChannels;
2857 
2858       int minV = _indexMin[m_NDims-1];
2859       int maxV = minV + (_indexMax[m_NDims-1]-_indexMin[m_NDims-1]);
2860 
2861       int cnt=0;
2862 
2863       // Read the previous lines
2864       for(i=0;i<minV;i++)
2865         {
2866         _stream->getline(s, 1024);
2867         }
2868 
2869       for(i=minV; i<=maxV; i+=1)
2870         {
2871         _stream->getline(s, 1024);
2872         if(!_stream->eof())
2873           {
2874           j = strlen(s)-1;
2875           while(j>0 && (isspace(s[j]) || !isprint(s[j])))
2876             {
2877             s[j--] = '\0';
2878             }
2879           if(usePath && !FileIsFullPath(s))
2880             {
2881             fName = pathName + s;
2882             }
2883           else
2884             {
2885             fName = s;
2886             }
2887 
2888           openReadStream(*readStreamTemp, fName);
2889           if(!readStreamTemp->is_open())
2890             {
2891             std::cerr << "MetaImage: Read: cannot open slice"
2892                                 << std::endl;
2893             continue;
2894             }
2895 
2896           // read only one slice
2897           int * indexMin = new int[m_NDims];
2898           int * indexMax = new int[m_NDims];
2899           quantity = 1;
2900           for(int k = 0;k<m_NDims-1;k++)
2901             {
2902             quantity *= _indexMax[k]-_indexMin[k]+1;
2903             indexMin[k]= _indexMin[k];
2904             indexMax[k]= _indexMax[k];
2905             }
2906           indexMin[m_NDims-1]=0;
2907           indexMax[m_NDims-1]=0;
2908 
2909           M_ReadElementsROI(readStreamTemp,
2910                              &(((char *)m_ElementData)[cnt*quantity*
2911                                                      elementSize]),
2912                              quantity, indexMin, indexMax,
2913                              subSamplingFactor,
2914                              m_SubQuantity[m_NDims-1]);
2915 
2916           cnt++;
2917           readStreamTemp->close();
2918           }
2919         }
2920       delete readStreamTemp;
2921       }
2922     else if(m_ElementDataFileName.find('%') != std::string::npos)
2923       {
2924       int elementSize;
2925       MET_SizeOfType(m_ElementType, &elementSize);
2926       elementSize *= m_ElementNumberOfChannels;
2927 
2928       int nWrds;
2929       char **wrds;
2930       int minV = 1;
2931       int maxV = m_DimSize[m_NDims-1];
2932       int stepV = 1;
2933       std::string s;
2934       std::ifstream* readStreamTemp = new std::ifstream;
2935       MET_StringToWordArray(m_ElementDataFileName.c_str(), &nWrds, &wrds);
2936       if(nWrds >= 2)
2937         {
2938         minV = (int)atof(wrds[1]);
2939         maxV = minV + m_DimSize[m_NDims-1] - 1;
2940         }
2941       if(nWrds >= 3)
2942         {
2943         maxV = (int)atof(wrds[2]);
2944         stepV = (maxV-minV)/(m_DimSize[m_NDims-1]);
2945         }
2946       if(nWrds >= 4)
2947         {
2948         stepV = (int)atof(wrds[3]);
2949         }
2950       if(nWrds >= 5 )
2951       {
2952         // In this case, the filename must have had spaces in the
2953         // name.  The filename was parsed into multiple pieces by the
2954         // MET_StringToWordArray, which parses based on spaces.
2955         // Thus, we need to reconstruct the filename in this case.
2956         // The last three wrds must be numbers.  If they are not, we give an error.
2957         for( i = nWrds-3; i < nWrds; i++ )
2958         {
2959           for( j = 0; j < strlen(wrds[i]); j++ )
2960           {
2961             if( !isdigit(wrds[i][j]) )
2962             {
2963               std::cerr << "MetaImage: Read: Last three arguments must be numbers!"
2964                   << std::endl;
2965               continue;
2966             }
2967           }
2968         }
2969         stepV = (int)atof(wrds[nWrds-1]);
2970         maxV =  (int)atof(wrds[nWrds-2]);
2971         minV =  (int)atof(wrds[nWrds-3]);
2972         for( i = 1; i < nWrds-3; i++ )
2973         {
2974           strcat(wrds[0]," ");
2975           strcat(wrds[0],wrds[i]);
2976         }
2977       }
2978       // If the specified size of the third dimension is less than the size
2979       // specified by the regular expression, we should only read a volume with the specified
2980       // size.  Otherwise, the code will crash when trying to fill m_ElementData more than it can hold.
2981       // Therefore, we modify maxV to ensure that the images spanned by minV:stepV:maxV are less than or equal
2982       // to the size in the last dimension.
2983       int numberOfImages = 1 + (maxV - minV)/stepV;
2984       if( numberOfImages > m_DimSize[m_NDims-1] )
2985       {
2986         maxV = (m_DimSize[m_NDims-1]-1)*stepV + minV;
2987       }
2988 
2989       int cnt = 0;
2990 
2991       // Uses the _indexMin and _indexMax
2992       minV += _indexMin[m_NDims-1];
2993       maxV = minV + (_indexMax[m_NDims-1]-_indexMin[m_NDims-1])*stepV;
2994 
2995       for(i=minV; i<=maxV; i += stepV)
2996         {
2997         s = string_format(wrds[0], i);
2998         if(usePath && !FileIsFullPath(s.c_str()))
2999           {
3000           fName = pathName + s;
3001           }
3002         else
3003           {
3004           fName = s;
3005           }
3006 
3007 
3008         openReadStream(*readStreamTemp, fName);
3009         if(!readStreamTemp->is_open())
3010           {
3011           std::cerr << "MetaImage: Read: cannot construct file"
3012                               << std::endl;
3013           continue;
3014           }
3015 
3016         // read only one slice
3017         int * indexMin = new int[m_NDims];
3018         int * indexMax = new int[m_NDims];
3019         quantity = 1;
3020         for(int k = 0;k<m_NDims-1;k++)
3021           {
3022           quantity *= _indexMax[k]-_indexMin[k]+1;
3023           indexMin[k]= _indexMin[k];
3024           indexMax[k]= _indexMax[k];
3025           }
3026         indexMin[m_NDims-1]=0;
3027         indexMax[m_NDims-1]=0;
3028 
3029         M_ReadElementsROI(readStreamTemp,
3030                        &(((char *)m_ElementData)[cnt*quantity*
3031                                                  elementSize]),
3032                        quantity, indexMin, indexMax,
3033                        subSamplingFactor,
3034                        m_SubQuantity[m_NDims-1]);
3035 
3036         cnt++;
3037 
3038         delete [] indexMin;
3039         delete [] indexMax;
3040 
3041         readStreamTemp->close();
3042         }
3043 
3044       for(i=0; i<nWrds; i++)
3045         {
3046         delete [] wrds[i];
3047         }
3048       delete [] wrds;
3049 
3050       delete readStreamTemp;
3051       }
3052     else
3053       {
3054       if(usePath && !FileIsFullPath(m_ElementDataFileName.c_str()))
3055         {
3056         fName = pathName + m_ElementDataFileName;
3057         }
3058       else
3059         {
3060         fName = m_ElementDataFileName;
3061         }
3062 
3063       std::ifstream* readStreamTemp = new std::ifstream;
3064 
3065       const char *extensions[] = { "", ".gz", ".Z", nullptr };
3066       for(unsigned ii = 0; extensions[ii] != nullptr; ii++)
3067         {
3068         std::string tempFName(fName);
3069         tempFName += extensions[ii];
3070         openReadStream(*readStreamTemp,tempFName.c_str());
3071         if(readStreamTemp->is_open())
3072           {
3073           if(ii > 0)
3074             {
3075             this->CompressedData(true);
3076             this->BinaryData(true);
3077             }
3078           break;
3079           }
3080         }
3081 
3082       if(!readStreamTemp->is_open())
3083         {
3084         std::cerr << "MetaImage: ReadROI: Cannot open data file"
3085                             << std::endl;
3086         if(m_ReadStream)
3087           {
3088           m_ReadStream->close();
3089           }
3090         delete readStreamTemp;
3091         return false;
3092         }
3093 
3094       M_ReadElementsROI(readStreamTemp, m_ElementData, quantity,
3095                         _indexMin, _indexMax, subSamplingFactor,
3096                         m_Quantity);
3097 
3098       readStreamTemp->close();
3099       delete readStreamTemp;
3100       }
3101     }
3102   return true;
3103 }
3104 
3105 /** Read an ROI */
3106 bool MetaImage::
M_ReadElementsROI(std::ifstream * _fstream,void * _data,std::streamoff _dataQuantity,int * _indexMin,int * _indexMax,unsigned int subSamplingFactor,std::streamoff _totalDataQuantity)3107 M_ReadElementsROI(std::ifstream * _fstream, void * _data,
3108                   std::streamoff _dataQuantity,
3109                   int* _indexMin, int* _indexMax,unsigned int subSamplingFactor,
3110                   std::streamoff _totalDataQuantity)
3111 {
3112   if(_totalDataQuantity ==0)
3113     {
3114     _totalDataQuantity = _dataQuantity;
3115     }
3116 
3117   for(int dim=0;dim<m_NDims;dim++)
3118     {
3119     _indexMin[dim] *= subSamplingFactor;
3120     _indexMax[dim] *= subSamplingFactor;
3121     }
3122 
3123 
3124   if(META_DEBUG)
3125     {
3126     std::cout << "MetaImage: M_ReadElementsROI" << std::endl;
3127     }
3128 
3129   if(m_HeaderSize>(int)0)
3130     {
3131     _fstream->seekg(m_HeaderSize, std::ios::beg);
3132     if(!_fstream->good())
3133       {
3134       std::cerr << "MetaImage: M_ReadElementsROI: header not read correctly"
3135                           << std::endl;
3136       return false;
3137       }
3138     }
3139 
3140   int elementSize;
3141   MET_SizeOfType(m_ElementType, &elementSize);
3142   std::streamoff readSize = _dataQuantity*m_ElementNumberOfChannels*elementSize;
3143   int elementNumberOfBytes = elementSize*m_ElementNumberOfChannels;
3144 
3145   if(META_DEBUG)
3146     {
3147     std::cout << "MetaImage: M_ReadElementsROI: ReadSize = "
3148                         << readSize << std::endl;
3149     }
3150 
3151   if(m_HeaderSize == -1)
3152     {
3153     if(META_DEBUG)
3154       {
3155       std::cout << "MetaImage: M_ReadElementsROI: Skipping header"
3156                           << std::endl;
3157       }
3158     std::streamoff headSize = _totalDataQuantity*m_ElementNumberOfChannels*elementSize;
3159     _fstream->seekg(-headSize, std::ios::end);
3160     }
3161 
3162   std::streampos dataPos = _fstream->tellg();
3163   std::streamoff i;
3164 
3165   // If compressed we inflate
3166   if(m_BinaryData && m_CompressedData)
3167     {
3168     // if m_CompressedDataSize is not defined we assume the size of the
3169     // file is the size of the compressed data
3170     if(m_CompressedDataSize==0)
3171       {
3172       _fstream->seekg(0, std::ios::end);
3173       m_CompressedDataSize = _fstream->tellg();
3174       _fstream->seekg(0, std::ios::beg);
3175       }
3176 
3177       unsigned char* data = static_cast<unsigned char*>(_data);
3178       // Initialize the index
3179       int* currentIndex = new int[m_NDims];
3180       for(i=0;i<m_NDims;i++)
3181         {
3182         currentIndex[i] = _indexMin[i];
3183         }
3184 
3185       // Optimize the size of the buffer to read depending on the
3186       // region shape
3187       // This calculate the number of continuous bytes in the file
3188       // which can be read
3189       std::streamoff elementsToRead = 1;
3190       int movingDirection = 0;
3191       do
3192         {
3193         elementsToRead *= _indexMax[movingDirection] - _indexMin[movingDirection] + 1;
3194         ++movingDirection;
3195         }
3196       while(subSamplingFactor == 1
3197             && movingDirection < m_NDims
3198             && _indexMin[movingDirection-1] == 0
3199             && _indexMax[movingDirection-1] == m_DimSize[movingDirection-1]-1);
3200 
3201       std::streamoff bytesToRead = elementsToRead*elementNumberOfBytes;
3202       std::streamoff gc = 0;
3203 
3204       bool done = false;
3205       while(!done)
3206         {
3207         // Seek to the right position
3208         std::streamoff seekoff = 0;
3209         for(i=0; i<m_NDims; i++)
3210           {
3211           seekoff += m_SubQuantity[i]*elementNumberOfBytes*currentIndex[i];
3212           }
3213 
3214 
3215         if(subSamplingFactor > 1)
3216           {
3217           unsigned char* subdata = new unsigned char[static_cast<size_t>(bytesToRead)];
3218           std::streamoff rOff =
3219             MET_UncompressStream(_fstream, seekoff, subdata,
3220                                  bytesToRead, m_CompressedDataSize,
3221                                  m_CompressionTable);
3222           // if there was a read error
3223           if(rOff == -1)
3224             {
3225             delete [] currentIndex;
3226             return false;
3227             }
3228 
3229           for(std::streamoff p=0;
3230               p<bytesToRead;
3231               p+=(subSamplingFactor*m_ElementNumberOfChannels*elementSize))
3232             {
3233             for(int s=0; s<m_ElementNumberOfChannels*elementSize; s++)
3234               {
3235               *data = subdata[p+s];
3236               gc++;
3237               data++;
3238               }
3239             }
3240           delete [] subdata;
3241           }
3242         else
3243           {
3244           std::streamoff rOff =
3245             MET_UncompressStream(_fstream, seekoff, data,
3246                                  bytesToRead, m_CompressedDataSize,
3247                                  m_CompressionTable);
3248           if(rOff == -1)
3249             {
3250             delete [] currentIndex;
3251             return false;
3252             }
3253           data += bytesToRead;
3254           gc += rOff;
3255           }
3256 
3257         if(gc == readSize)
3258           {
3259           break;
3260           }
3261 
3262         // Go forward
3263         if(m_NDims == 1)
3264           {
3265           break;
3266           }
3267 
3268         currentIndex[movingDirection] += subSamplingFactor;
3269 
3270         // Check if we are still in the region
3271         for(i=1;i<m_NDims;i++)
3272           {
3273           if(currentIndex[i]>_indexMax[i])
3274             {
3275             if(i==m_NDims-1)
3276               {
3277               done = true;
3278               break;
3279               }
3280             else
3281               {
3282               currentIndex[i] = _indexMin[i];
3283               currentIndex[i+1] += subSamplingFactor;
3284               }
3285             }
3286           }
3287         }
3288 
3289       if(gc != readSize)
3290         {
3291         std::cerr
3292                   << "MetaImage: M_ReadElementsROI: data not read completely"
3293                   << std::endl;
3294         std::cerr << "   ideal = " << readSize << " : actual = " << gc
3295                   << std::endl;
3296         delete [] currentIndex;
3297         return false;
3298         }
3299 
3300       delete [] currentIndex;
3301     }
3302   else // if not compressed
3303     {
3304     double tf;
3305     MET_SizeOfType(m_ElementType, &elementSize);
3306 
3307     char* data = static_cast<char*>(_data);
3308     // Initialize the index
3309     int* currentIndex = new int[m_NDims];
3310     for(i=0;i<m_NDims;i++)
3311       {
3312       currentIndex[i] = _indexMin[i];
3313       }
3314 
3315     // Optimize the size of the buffer to read depending on the
3316     // region shape
3317     // This calculate the number of continuous bytes in the file
3318     // which can be read
3319     std::streamoff elementsToRead = 1;
3320     int movingDirection = 0;
3321     do
3322       {
3323       elementsToRead *= _indexMax[movingDirection] - _indexMin[movingDirection] + 1;
3324       ++movingDirection;
3325       }
3326     while(subSamplingFactor == 1
3327           && movingDirection < m_NDims
3328           && _indexMin[movingDirection-1] == 0
3329           && _indexMax[movingDirection-1] == m_DimSize[movingDirection-1]-1);
3330 
3331     //readLine *= m_ElementNumberOfChannels*elementSize;
3332     std::streamoff gc = 0;
3333 
3334     bool done = false;
3335     while(!done)
3336       {
3337       // Seek to the right position
3338       std::streamoff seekoff = 0;
3339       for(i=0;i<m_NDims;i++)
3340         {
3341         seekoff += m_SubQuantity[i]*m_ElementNumberOfChannels*elementSize*currentIndex[i];
3342         }
3343 
3344       _fstream->seekg(dataPos+seekoff, std::ios::beg);
3345 
3346       // Read a line
3347       if(subSamplingFactor > 1)
3348         {
3349         if(!m_BinaryData) // Not binary data
3350           {
3351           for(i=0; i<elementsToRead; i+=subSamplingFactor)
3352             {
3353             *_fstream >> tf;
3354             MET_DoubleToValue(tf, m_ElementType, _data, i);
3355 
3356             for(unsigned int j=0;j<subSamplingFactor;j++)
3357               {
3358               _fstream->get();
3359               }
3360             }
3361           }
3362         else // Binary data
3363           {
3364           char* subdata = new char[static_cast<size_t>(elementsToRead*elementNumberOfBytes)];
3365 
3366           _fstream->read(subdata, size_t(elementsToRead*elementNumberOfBytes));
3367 
3368           for(std::streamoff p=0;
3369               p<elementsToRead*elementNumberOfBytes;
3370               p+=(subSamplingFactor*elementNumberOfBytes))
3371             {
3372             for(int s=0;s<elementNumberOfBytes;s++)
3373               {
3374               *data = subdata[p+s];
3375               gc++;
3376               data++;
3377               }
3378             }
3379           delete [] subdata;
3380           }
3381         }
3382       else
3383         {
3384         if(!m_BinaryData) // Not binary data
3385           {
3386           // anyone using ROI reading of ASCII??
3387           // does this work? what about incrementing data?
3388           // what about data sizes and random access of file?
3389           std::streamoff blockSize = elementsToRead*m_ElementNumberOfChannels*elementSize;
3390           M_ReadElementData(  _fstream, data, (size_t)blockSize  );
3391           gc += blockSize;
3392 
3393           }
3394         else // binary data
3395           {
3396 
3397           M_ReadElementData(  _fstream, data, elementsToRead );
3398           gc += elementsToRead*elementNumberOfBytes;
3399           data += elementsToRead*elementNumberOfBytes;
3400           }
3401         }
3402 
3403       // I don't think this check is really needed -BCL
3404       if(gc == readSize)
3405         {
3406         break;
3407         }
3408 
3409       // check if there is only one read needed
3410       if ( movingDirection >= m_NDims )
3411         {
3412         break;
3413         }
3414 
3415       // Go forward
3416       currentIndex[movingDirection] += subSamplingFactor;
3417 
3418       // Check if we are still in the region
3419       for(i=movingDirection;i<m_NDims;i++)
3420         {
3421         if(currentIndex[i]>_indexMax[i])
3422           {
3423           if(i==m_NDims-1)
3424             {
3425             done = true;
3426             break;
3427             }
3428           else
3429             {
3430             currentIndex[i] = _indexMin[i];
3431             currentIndex[i+1] += subSamplingFactor;
3432             }
3433           }
3434         }
3435       }
3436 
3437     delete [] currentIndex;
3438 
3439     if(gc != readSize)
3440       {
3441       std::cerr
3442                 << "MetaImage: M_ReadElementsROI: data not read completely"
3443                 << std::endl;
3444       std::cerr << "   ideal = " << readSize << " : actual = " << gc
3445                 << std::endl;
3446       return false;
3447       }
3448     }
3449 
3450   return true;
3451 }
3452 
3453 
3454 bool MetaImage::
M_ReadElementData(std::ifstream * _fstream,void * _data,std::streamoff _dataQuantity)3455 M_ReadElementData(std::ifstream * _fstream,
3456                   void * _data,
3457                   std::streamoff _dataQuantity)
3458 {
3459   // NOTE: this method is different from WriteElementData
3460   std::streamoff gc = 0;
3461 
3462   if(!m_BinaryData)
3463     {
3464     double tf;
3465 
3466     for(int i=0; i<_dataQuantity; i++)
3467       {
3468       *_fstream >> tf;
3469       MET_DoubleToValue(tf, m_ElementType, _data, i);
3470       _fstream->get();
3471       ++gc;
3472       }
3473     }
3474   else
3475     {
3476     if(m_CompressedData)
3477       {
3478 
3479       // the data is read with calls no bigger then MaxIOChunk
3480       std::streamoff bytesRemaining = _dataQuantity;
3481       while ( bytesRemaining )
3482         {
3483         std::streamoff chunkToRead = bytesRemaining > MaxIOChunk ? MaxIOChunk : bytesRemaining;
3484         _fstream->read( (char *)_data, (size_t)chunkToRead );
3485         _data = (char *)(_data) + chunkToRead;
3486         bytesRemaining -= chunkToRead;
3487         std::streamsize numberOfBytesRead = _fstream->gcount();
3488         gc += numberOfBytesRead;
3489         }
3490 
3491       }
3492     else
3493       {
3494       int elementSize;
3495       MET_SizeOfType(m_ElementType, &elementSize);
3496       std::streamoff elementNumberOfBytes = elementSize*m_ElementNumberOfChannels;
3497 
3498       // the data is read with calls no bigger than MaxIOChunk
3499       std::streamoff bytesRemaining = _dataQuantity * elementNumberOfBytes;
3500       while ( bytesRemaining )
3501         {
3502         std::streamoff chunkToRead = bytesRemaining > MaxIOChunk ? MaxIOChunk : bytesRemaining;
3503         _fstream->read( (char *)_data, (size_t)chunkToRead );
3504         _data = (char *)(_data) + chunkToRead;
3505         bytesRemaining -= chunkToRead;
3506         std::streamsize numberOfBytesRead = _fstream->gcount();
3507         gc += numberOfBytesRead;
3508         }
3509       // convert to number of bytes so that it'll match gc's units
3510       _dataQuantity *= elementNumberOfBytes;
3511       }
3512     }
3513 
3514   // check that we actually read the correct number of bytes
3515   if( gc != _dataQuantity )
3516     {
3517     std::cerr
3518       << "MetaImage: M_ReadElementsData: data not read completely"
3519       << std::endl;
3520     std::cerr << "   ideal = " << _dataQuantity << " : actual = " << gc
3521                         << std::endl;
3522     return false;
3523     }
3524 
3525   // check if the io stream did not fail in the process of reading
3526   if ( _fstream->fail() )
3527     {
3528     std::cerr
3529       << "MetaImage: M_ReadElementsData: file stream is fail after read"
3530       << std::endl;
3531     return false;
3532     }
3533 
3534   return true;
3535 }
3536 
3537 
3538 #if (METAIO_USE_NAMESPACE)
3539 }
3540 #endif
3541