Merge branch 'release'

Source/DataStructureAndEncodingDefinition/gdcmElement.h

@@ -473,7 +473,7 @@ template<> class EncodingImplementation<VR::VRBINARY> {
     assert( _is ); // Is stream valid ?
     _is.read( reinterpret_cast<char*>(data+0), type_size);
     for(unsigned long i=1; i<length; ++i) {
-      assert( _is );
+      if( _is )
       _is.read( reinterpret_cast<char*>(data+i), type_size );
     }
     //ByteSwap<T>::SwapRangeFromSwapCodeIntoSystem(data,
@@ -489,7 +489,7 @@ template<> class EncodingImplementation<VR::VRBINARY> {
     assert( _is ); // Is stream valid ?
     _is.read( reinterpret_cast<char*>(data+0), type_size);
     for(unsigned long i=1; i<length; ++i) {
-      assert( _is );
+      if( _is )
       _is.read( reinterpret_cast<char*>(data+i), type_size );
     }
     //ByteSwap<T>::SwapRangeFromSwapCodeIntoSystem(data,

Source/MediaStorageAndFileFormat/gdcmDirectionCosines.cxx

@@ -98,15 +98,20 @@ double DirectionCosines::Dot() const
 }
 
 // static function is within gdcm:: namespace, so should not pollute too much on UNIX
-static inline double Norm(const double x[3])
+static inline double NormImpl(const double x[3])
 {
   return sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
 }
 
+double DirectionCosines::Norm(const double v[3])
+{
+  return NormImpl(v);
+}
+
 void DirectionCosines::Normalize(double v[3])
 {
   double den;
-  if ( (den = Norm(v)) != 0.0 )
+  if ( (den = NormImpl(v)) != 0.0 )
     {
     for (int i=0; i < 3; i++)
       {
@@ -119,15 +124,15 @@ void DirectionCosines::Normalize()
 {
   double *x = Values;
   double den;
-  if ( (den = Norm(x)) != 0.0 )
+  if ( (den = NormImpl(x)) != 0.0 )
     {
     for (int i=0; i < 3; i++)
       {
       x[i] /= den;
       }
     }
   x = Values+3;
-  if ( (den = Norm(x)) != 0.0 )
+  if ( (den = NormImpl(x)) != 0.0 )
     {
     for (int i=0; i < 3; i++)
       {

Source/MediaStorageAndFileFormat/gdcmDirectionCosines.h

@@ -49,6 +49,9 @@ class GDCM_EXPORT DirectionCosines
   /// Normalize in-place
   static void Normalize(double v[3]);
 
+  /// Return norm of the vector
+  static double Norm(const double v[3]);
+
   /// Make the class behave like a const double *
   operator const double* () const { return Values; }
 

Source/MediaStorageAndFileFormat/gdcmImageChangeTransferSyntax.cxx

@@ -421,6 +421,7 @@ bool ImageChangeTransferSyntax::Change()
     if( !b )
       {
       gdcmErrorMacro( "Error in getting buffer from input image." );
+      delete bv0;
       return false;
       }
     pixeldata.SetValue( *bv0 );

Source/MediaStorageAndFileFormat/gdcmImageCodec.cxx

@@ -675,6 +675,7 @@ bool ImageCodec::DecodeByStreams(std::istream &is, std::ostream &os)
 
   // Do the overlay cleanup (cleanup the unused bits)
   // must be the last operation (duh!)
+  bool copySuccess = false;
   if ( PF.GetBitsAllocated() != PF.GetBitsStored()
     && PF.GetBitsAllocated() != 8 )
     {
@@ -685,21 +686,23 @@ bool ImageCodec::DecodeByStreams(std::istream &is, std::ostream &os)
     // Sigh, I finally found someone not declaring that unused bits where not zero:
     // gdcmConformanceTests/dcm4chee_unusedbits_not_zero.dcm
     if( NeedOverlayCleanup )
-      DoOverlayCleanup(*cur_is,os);
+      {
+      copySuccess = DoOverlayCleanup(*cur_is, os);
+      }
     else
       {
       // Once the issue with IMAGES/JPLY/RG3_JPLY aka gdcmData/D_CLUNIE_RG3_JPLY.dcm is solved the previous
       // code will be replace with a simple call to:
-      DoSimpleCopy(*cur_is,os);
+      copySuccess = DoSimpleCopy(*cur_is, os);
       }
     }
   else
     {
     assert( PF.GetBitsAllocated() == PF.GetBitsStored() );
-    DoSimpleCopy(*cur_is,os);
+    copySuccess = DoSimpleCopy(*cur_is, os);
     }
 
-  return true;
+  return copySuccess;
 }
 
 bool ImageCodec::IsValid(PhotometricInterpretation const &)

Source/MediaStorageAndFileFormat/gdcmImageHelper.cxx

@@ -205,6 +205,7 @@ static bool ComputeZSpacingFromIPP(const DataSet &ds, double &zspacing)
   double normal[3];
   DirectionCosines dc( cosines.data() );
   dc.Cross( normal );
+  DirectionCosines::Normalize(normal);
 
   // For each item
   SequenceOfItems::SizeType nitems = sqi->GetNumberOfItems();
@@ -2027,6 +2028,7 @@ void ImageHelper::SetOriginValue(DataSet & ds, const Image & image)
 
     double normal[3];
     dc.Cross( normal );
+    DirectionCosines::Normalize(normal);
 
     for(unsigned int i = 0; i < dimz; ++i )
       {

Source/MediaStorageAndFileFormat/gdcmJPEG2000Codec.cxx

@@ -826,8 +826,13 @@ std::pair<char *, size_t> JPEG2000Codec::DecodeByStreamsCommon(char *dummy_buffe
 
     // ELSCINT1_JP2vsJ2K.dcm
     // -> prec = 12, bpp = 0, sgnd = 0
-    //assert( wr == Dimensions[0] );
-    //assert( hr == Dimensions[1] );
+    if( wr != Dimensions[0] || hr != Dimensions[1] ) {
+	    gdcmErrorMacro("Invalid dimension");
+	    delete[] raw;
+    opj_destroy_codec(dinfo);
+  opj_image_destroy(image);
+    return std::pair<char*,size_t>(nullptr,0);
+    }
     if( comp->sgnd != PF.GetPixelRepresentation() )
       {
       PF.SetPixelRepresentation( (uint16_t)comp->sgnd );

Source/MediaStorageAndFileFormat/gdcmLookupTable.cxx

@@ -130,7 +130,10 @@ void LookupTable::SetLUT(LookupTableType type, const unsigned char *array,
 
   if( !IncompleteLUT )
     {
-    assert( Internal->RGB.size() == 3*Internal->Length[type]*(BitSample/8) );
+    if( Internal->RGB.size() != 3*Internal->Length[type]*(BitSample/8) ) {
+      gdcmErrorMacro( "Invalid length for LUT data" );
+      return;
+    }
     }
   // Too funny: 05115014-mr-siemens-avanto-syngo-with-palette-icone.dcm
   // There is pseudo PALETTE_COLOR LUT in the Icon, if one look carefully the LUT values

Source/MediaStorageAndFileFormat/gdcmPixmapReader.cxx

@@ -306,8 +306,12 @@ static void DoIconImage(const DataSet& rootds, Pixmap& image)
           unsigned long check =
             (el_us3.GetValue(0) ? el_us3.GetValue(0) : 65536)
             * el_us3.GetValue(2) / 8;
-          assert( check == lut_raw->GetLength() || 2 * check == lut_raw->GetLength()
-            || check + 1 == lut_raw->GetLength() ); (void)check;
+          if(!( check == lut_raw->GetLength() || 2 * check == lut_raw->GetLength()
+            || check + 1 == lut_raw->GetLength() )) {
+          gdcmErrorMacro( "Icon Sequence is invalid. Giving up" );
+          pixeldata.Clear();
+          return;
+	  }
           }
         else if( ds.FindDataElement( seglut ) )
           {

Source/MediaStorageAndFileFormat/gdcmRAWCodec.cxx

@@ -112,9 +112,7 @@ bool RAWCodec::DecodeBytes(const char* inBytes, size_t inBufferLength,
   if(!r) return false;
 
   std::string str = os.str();
-  //std::string::size_type check = str.size();//unused
 
-
   if( this->GetPixelFormat() == PixelFormat::UINT12 ||
       this->GetPixelFormat() == PixelFormat::INT12 )
     {
@@ -135,7 +133,14 @@ bool RAWCodec::DecodeBytes(const char* inBytes, size_t inBufferLength,
     // DermaColorLossLess.dcm
     //assert (check == inOutBufferLength || check == inOutBufferLength + 1);
     // problem with: SIEMENS_GBS_III-16-ACR_NEMA_1.acr
-    memcpy(outBytes, str.c_str(), inOutBufferLength);
+    size_t len = str.size();
+    if( inOutBufferLength <= len )
+      memcpy(outBytes, str.c_str(), inOutBufferLength);
+    else
+    {
+      gdcmWarningMacro( "Requesting too much data. Truncating result" );
+      memcpy(outBytes, str.c_str(), len);
+    }
     }
 
   return r;

Source/MediaStorageAndFileFormat/gdcmSplitMosaicFilter.cxx

@@ -287,6 +287,96 @@ bool SplitMosaicFilter::ComputeMOSAICSliceNormal( double slicenormalvector[3], b
   return snvfound;
 }
 
+bool SplitMosaicFilter::ComputeMOSAICImagePositionPatient( double ret[3], 
+    const double ipp[6],
+    const double dircos[6],
+    const double pixelspacing[3],
+    const unsigned int image_dims[3] ,
+    const unsigned int mosaic_dims[3] , bool inverted)
+{
+  CSAHeader csa;
+  DataSet& ds = GetFile().GetDataSet();
+  DirectionCosines dc( dircos );
+  dc.Normalize();
+  const double *dircos_normalized = dc;
+  const double *x = dircos_normalized;
+  const double *y = dircos_normalized + 3;
+
+  double ipp_csa[3] = {};
+  bool hasIppCsa = false;
+  MrProtocol mrprot;
+  // https://www.nmr.mgh.harvard.edu/~greve/dicom-unpack
+  if( csa.GetMrProtocol(ds, mrprot) ) 
+  {
+    MrProtocol::SliceArray sa;
+    bool b = mrprot.GetSliceArray(sa);
+    if( b ) {
+      size_t size = sa.Slices.size();
+      if( size ) {
+        // two cases:
+        if( size == mosaic_dims[2] ) {
+          // all mosaic have there own slice position, simply need to pick the right one
+          // Handle inverted case:
+          size_t index = inverted ? size - 1 : 0;
+          MrProtocol::Slice & slice = sa.Slices[index];
+          MrProtocol::Vector3 & p = slice.Position;
+          double pos[3];
+          pos[0] = p.dSag;
+          pos[1] = p.dCor;
+          pos[2] = p.dTra;
+          for(int i = 0; i < 3; ++i ) {
+            ipp_csa[i] = pos[i] - mosaic_dims[0] / 2. * pixelspacing[0] * x[i] - mosaic_dims[1] / 2. * pixelspacing[1] * y[i];
+          }
+          hasIppCsa = true;
+        } else if( size == 1 /*&& mosaic_dims[2] % 2 == 0*/) {
+          // there is a single SliceArray but multiple mosaics, assume this is exactly the center one
+          double z[3]={};
+          dc.Cross(z);
+          DirectionCosines::Normalize(z);
+          size_t index = 0;
+          MrProtocol::Slice & slice = sa.Slices[index];
+          MrProtocol::Vector3 & p = slice.Position;
+          double pos[3];
+          pos[0] = p.dSag;
+          pos[1] = p.dCor;
+          pos[2] = p.dTra;
+          for(int i = 0; i < 3; ++i ) {
+            ipp_csa[i] = pos[i] - mosaic_dims[0] / 2. * pixelspacing[0] * x[i] - mosaic_dims[1] / 2. * pixelspacing[1] * y[i]
+              - (mosaic_dims[2] - 1) / 2. * pixelspacing[2] * z[i];
+          }
+          hasIppCsa = true;
+        } else {
+          gdcmWarningMacro( "Inconsistent SliceArray: " << size << " vs expected: " << mosaic_dims[2] );
+        }
+      }
+    }
+  }
+
+  // https://nipy.org/nibabel/dicom/dicom_mosaic.html#dicom-mosaic
+  double ipp_dcm[3] = {};
+  {
+    for(int i = 0; i < 3; ++i ) {
+      ipp_dcm[i] = ipp[i] + (image_dims[0]  - mosaic_dims[0]) / 2. * pixelspacing[0] * x[i] +
+        (image_dims[1]  - mosaic_dims[1]) / 2. * pixelspacing[1] * y[i] ;
+    }
+  }
+  if(hasIppCsa ) {
+    double diff[3];
+    for(int i = 0; i < 3; ++i ) {
+      diff[i] = ipp_dcm[i] - ipp_csa[i];
+    }
+    double n = DirectionCosines::Norm(diff);
+    if( n > 1e-4 ) {
+      gdcmWarningMacro( "Inconsistent values for IPP/CSA: (" << ipp_dcm[0] << "," << ipp_dcm[1] << "," << ipp_dcm[2] << ") vs (" << ipp_csa[0] << "," << ipp_csa[1] << "," << ipp_csa[2] << ")" );
+    }
+  }
+  // no error set origin:
+  for(int i = 0; i < 3; ++i )
+    ret[i] = ipp_dcm[i];
+
+  return true;
+}
+
 bool SplitMosaicFilter::ComputeMOSAICSlicePosition( double pos[3], bool )
 {
   CSAHeader csa;
@@ -301,25 +391,6 @@ bool SplitMosaicFilter::ComputeMOSAICSlicePosition( double pos[3], bool )
 
   size_t size = sa.Slices.size();
   if( !size ) return false;
-#if 0
-  {
-    double z[3];
-    for( int i = 0; i < size; ++i )
-    {
-      MrProtocol::Slice & slice = sa.Slices[i];
-      MrProtocol::Vector3 & p = slice.Position;
-      z[0] = p.dSag;
-      z[1] = p.dCor;
-      z[2] = p.dTra;
-      const double snv_dot = DirectionCosines::Dot( slicenormalvector, z );
-      if( (1. - snv_dot) < 1e-6 )
-      {
-        gdcmErrorMacro("Invalid direction found");
-        return false;
-      }
-    }
-  }
-#endif
 
   size_t index = 0;
   MrProtocol::Slice & slice = sa.Slices[index];
@@ -353,13 +424,21 @@ bool SplitMosaicFilter::Split()
   }
   (void)hasNormalCSA;
   double origin[3];
+  const Image &inputimage = GetImage();
+#if 0
   if( !ComputeMOSAICSlicePosition( origin, inverted ) )
   {
     gdcmWarningMacro( "Origin will not be accurate" );
     hasOriginCSA = false;
   }
+#else
+  if(!ComputeMOSAICImagePositionPatient( origin, inputimage.GetOrigin(), inputimage.GetDirectionCosines(), inputimage.GetSpacing(), inputimage.GetDimensions(), dims, inverted ))
+  {
+    gdcmWarningMacro( "Origin will not be accurate" );
+    hasOriginCSA = false;
+  }
+#endif
 
-  const Image &inputimage = GetImage();
   if( inputimage.GetPixelFormat() != PixelFormat::UINT16 )
     {
     gdcmErrorMacro( "Expecting UINT16 PixelFormat" );

Source/MediaStorageAndFileFormat/gdcmSplitMosaicFilter.h

@@ -57,8 +57,17 @@ class GDCM_EXPORT SplitMosaicFilter
   bool ComputeMOSAICSliceNormal( double dims[3], bool & inverted );
 
   /// Extract the value for ImagePositionPatient (requires inverted flag)
+  /// Deprecated
   bool ComputeMOSAICSlicePosition( double pos[3], bool inverted );
 
+  /// Extract the value for ImagePositionPatient
+  bool ComputeMOSAICImagePositionPatient( double pos[3],
+    const double ipp[6],
+    const double dircos[6],
+    const double pixelspacing[3],
+    const unsigned int image_dims[3] ,
+    const unsigned int mosaic_dims[3], bool inverted );
+
   void SetImage(const Image& image);
   const Image &GetImage() const { return *I; }
   Image &GetImage() { return *I; }