diff --git a/Source/DataStructureAndEncodingDefinition/gdcmElement.h b/Source/DataStructureAndEncodingDefinition/gdcmElement.h index b49b093dc..15fb3a117 100644 --- a/Source/DataStructureAndEncodingDefinition/gdcmElement.h +++ b/Source/DataStructureAndEncodingDefinition/gdcmElement.h @@ -473,7 +473,7 @@ template<> class EncodingImplementation { assert( _is ); // Is stream valid ? _is.read( reinterpret_cast(data+0), type_size); for(unsigned long i=1; i(data+i), type_size ); } //ByteSwap::SwapRangeFromSwapCodeIntoSystem(data, @@ -489,7 +489,7 @@ template<> class EncodingImplementation { assert( _is ); // Is stream valid ? _is.read( reinterpret_cast(data+0), type_size); for(unsigned long i=1; i(data+i), type_size ); } //ByteSwap::SwapRangeFromSwapCodeIntoSystem(data, diff --git a/Source/MediaStorageAndFileFormat/gdcmDirectionCosines.cxx b/Source/MediaStorageAndFileFormat/gdcmDirectionCosines.cxx index 275d74855..119a7b4a8 100644 --- a/Source/MediaStorageAndFileFormat/gdcmDirectionCosines.cxx +++ b/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,7 +124,7 @@ 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++) { @@ -127,7 +132,7 @@ void DirectionCosines::Normalize() } } x = Values+3; - if ( (den = Norm(x)) != 0.0 ) + if ( (den = NormImpl(x)) != 0.0 ) { for (int i=0; i < 3; i++) { diff --git a/Source/MediaStorageAndFileFormat/gdcmDirectionCosines.h b/Source/MediaStorageAndFileFormat/gdcmDirectionCosines.h index 73717fdc0..1313146cf 100644 --- a/Source/MediaStorageAndFileFormat/gdcmDirectionCosines.h +++ b/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; } diff --git a/Source/MediaStorageAndFileFormat/gdcmImageChangeTransferSyntax.cxx b/Source/MediaStorageAndFileFormat/gdcmImageChangeTransferSyntax.cxx index fcb61e611..9457c5e9b 100644 --- a/Source/MediaStorageAndFileFormat/gdcmImageChangeTransferSyntax.cxx +++ b/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 ); diff --git a/Source/MediaStorageAndFileFormat/gdcmImageCodec.cxx b/Source/MediaStorageAndFileFormat/gdcmImageCodec.cxx index c2030c0a8..2c25a90f9 100644 --- a/Source/MediaStorageAndFileFormat/gdcmImageCodec.cxx +++ b/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 &) diff --git a/Source/MediaStorageAndFileFormat/gdcmImageHelper.cxx b/Source/MediaStorageAndFileFormat/gdcmImageHelper.cxx index e32e46b0e..3dd08ed97 100644 --- a/Source/MediaStorageAndFileFormat/gdcmImageHelper.cxx +++ b/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 ) { diff --git a/Source/MediaStorageAndFileFormat/gdcmJPEG2000Codec.cxx b/Source/MediaStorageAndFileFormat/gdcmJPEG2000Codec.cxx index 10ac23cca..430a24a87 100644 --- a/Source/MediaStorageAndFileFormat/gdcmJPEG2000Codec.cxx +++ b/Source/MediaStorageAndFileFormat/gdcmJPEG2000Codec.cxx @@ -826,8 +826,13 @@ std::pair 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(nullptr,0); + } if( comp->sgnd != PF.GetPixelRepresentation() ) { PF.SetPixelRepresentation( (uint16_t)comp->sgnd ); diff --git a/Source/MediaStorageAndFileFormat/gdcmLookupTable.cxx b/Source/MediaStorageAndFileFormat/gdcmLookupTable.cxx index 0d5a99c40..2c566923b 100644 --- a/Source/MediaStorageAndFileFormat/gdcmLookupTable.cxx +++ b/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 diff --git a/Source/MediaStorageAndFileFormat/gdcmPixmapReader.cxx b/Source/MediaStorageAndFileFormat/gdcmPixmapReader.cxx index 9c30ff8b9..258a23c1f 100644 --- a/Source/MediaStorageAndFileFormat/gdcmPixmapReader.cxx +++ b/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 ) ) { diff --git a/Source/MediaStorageAndFileFormat/gdcmRAWCodec.cxx b/Source/MediaStorageAndFileFormat/gdcmRAWCodec.cxx index 19f739399..46392461e 100644 --- a/Source/MediaStorageAndFileFormat/gdcmRAWCodec.cxx +++ b/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; diff --git a/Source/MediaStorageAndFileFormat/gdcmSplitMosaicFilter.cxx b/Source/MediaStorageAndFileFormat/gdcmSplitMosaicFilter.cxx index 1b1e6e6c9..6f368e0c5 100644 --- a/Source/MediaStorageAndFileFormat/gdcmSplitMosaicFilter.cxx +++ b/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" ); diff --git a/Source/MediaStorageAndFileFormat/gdcmSplitMosaicFilter.h b/Source/MediaStorageAndFileFormat/gdcmSplitMosaicFilter.h index 7f30b0a25..06f281a3d 100644 --- a/Source/MediaStorageAndFileFormat/gdcmSplitMosaicFilter.h +++ b/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; }