diff --git a/core/CHANGELOG b/core/CHANGELOG
index b56a72af..03974615 100644
--- a/core/CHANGELOG
+++ b/core/CHANGELOG
@@ -1,3 +1,11 @@
+Release CMTK-3.4.0 (2023-09-08):
+
+(r5471) Updated test snapshots to accomodate updated Siemens support.
+
+(r5468) Added support for extracting values from:
+ * Siemens MOSAIC format
+ * Siemens CSA DICOM headers
+
Release CMTK-3.3.2 (2022-11-06):
(r5456) Merged 'ndar' branch:
diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt
index 7c46226c..bc875587 100644
--- a/core/CMakeLists.txt
+++ b/core/CMakeLists.txt
@@ -45,8 +45,8 @@ SET(CMAKE_OVERRIDE_COMPILER_MISMATCH 1)
PROJECT(CMTK)
SET(CMTK_VERSION_MAJOR "3")
-SET(CMTK_VERSION_MINOR "3")
-SET(CMTK_VERSION_PATCH "2")
+SET(CMTK_VERSION_MINOR "4")
+SET(CMTK_VERSION_PATCH "0")
SET(CMTK_VERSION_STRING "${CMTK_VERSION_MAJOR}.${CMTK_VERSION_MINOR}.${CMTK_VERSION_PATCH}")
diff --git a/core/libs/IO/cmtkDICOM.cxx b/core/libs/IO/cmtkDICOM.cxx
index 2aa75c27..0891301a 100644
--- a/core/libs/IO/cmtkDICOM.cxx
+++ b/core/libs/IO/cmtkDICOM.cxx
@@ -1,535 +1,515 @@
-/*
-//
-// Copyright 2014-2015 Google Inc.
-//
-// Copyright 2004-2014 SRI International
-//
-// Copyright 1997-2009 Torsten Rohlfing
-//
-// This file is part of the Computational Morphometry Toolkit.
-//
-// http://www.nitrc.org/projects/cmtk/
-//
-// The Computational Morphometry Toolkit is free software: you can
-// redistribute it and/or modify it under the terms of the GNU General Public
-// License as published by the Free Software Foundation, either version 3 of
-// the License, or (at your option) any later version.
-//
-// The Computational Morphometry Toolkit is distributed in the hope that it
-// will be useful, but WITHOUT ANY WARRANTY; without even the implied
-// warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License along
-// with the Computational Morphometry Toolkit. If not, see
-// .
-//
-// $Revision$
-//
-// $LastChangedDate$
-//
-// $LastChangedBy$
-//
-*/
-
-#include "cmtkDICOM.h"
-
-#include
-#include
-
-#include
-#include
-
-#include
-
-
-#ifdef CMTK_USE_DCMTK_JPEG
-# include
-#endif
-
-#ifdef HAVE_SYS_TYPES_H
-# include
-#endif
-
-#include
-#include
-#include
-
-#ifndef DCM_ACR_NEMA_ImagePosition
-# define DCM_ACR_NEMA_ImagePosition DcmTagKey(0x0020, 0x0030)
-#endif
-
-#ifndef DCM_ACR_NEMA_ImageOrientation
-# define DCM_ACR_NEMA_ImageOrientation DcmTagKey(0x0020, 0x0035)
-#endif
-
-#ifndef DCM_ACR_NEMA_Location
-# define DCM_ACR_NEMA_Location DcmTagKey(0x0020, 0x0050)
-#endif
-
-#ifndef DCM_ACR_NEMA_2C_VariablePixelData
-# define DCM_ACR_NEMA_2C_VariablePixelData DcmTagKey(0x7f00, 0x0010)
-#endif
-
-namespace
-cmtk
-{
-
-/** \addtogroup IO */
-//@{
-
-void
-DICOM::InitFromFile( const std::string& path )
-{
- this->m_Path = path;
-
-#ifdef CMTK_USE_DCMTK_JPEG
- // register global decompression codecs
- static bool decodersRegistered = false;
- if ( ! decodersRegistered )
- {
- DJDecoderRegistration::registerCodecs( EDC_photometricInterpretation, EUC_default, EPC_default, 1 );
- decodersRegistered = true;
- }
-#endif
-
- std::auto_ptr fileformat( new DcmFileFormat );
- if (!fileformat.get())
- {
- throw Exception( "Could not create DICOM file format object." );
- }
-
- OFCondition status = fileformat->loadFile( path.c_str() );
-
- if ( !status.good() )
- {
- throw Exception( "Cannot read DICOM file.." );
-// StdErr << "Error: cannot read DICOM file " << path << " (" << status.text() << ")\n";
- }
-
- this->m_Dataset = fileformat->getAndRemoveDataset();
-
- if ( !this->m_Dataset )
- {
- throw Exception( "File format has NULL dataset." );
- }
-
- this->m_Document = std::auto_ptr( new DiDocument( this->m_Dataset, this->m_Dataset->getOriginalXfer(), CIF_AcrNemaCompatibility ) );
- if ( ! this->m_Document.get() || ! this->m_Document->good() )
- {
- throw Exception( "Could not create document representation." );
- }
-}
-
-const FixedVector<3,int>
-DICOM::GetDims() const
-{
- FixedVector<3,int> dims( 0 );
-
- Uint16 tempUint16 = 1;
- if ( this->Document().getValue( DCM_Columns, tempUint16 ) )
- {
- dims[0] = static_cast( tempUint16 );
- }
-
- if ( this->Document().getValue( DCM_Rows, tempUint16 ) )
- {
- dims[1] = static_cast( tempUint16 );
- }
-
- // detect and treat multi-frame files
- if ( ! this->Document().getValue( DCM_NumberOfFrames, tempUint16 ) )
- {
- // unlike Rows/Columns, NumberofFrames defaults to 1
- tempUint16 = 1;
- }
- dims[2] = tempUint16;
-
- return dims;
-}
-
-const FixedVector<3,double>
-DICOM::GetPixelSize() const
-{
- FixedVector<3,double> pixelSize( 0.0 );
-
- // get calibration from image
- const bool hasPixelSpacing = (this->Document().getValue(DCM_PixelSpacing, pixelSize[0], 0) > 0);
- if ( hasPixelSpacing )
- {
- if (this->Document().getValue(DCM_PixelSpacing, pixelSize[1], 1) < 2)
- {
- throw Exception( "DICOM file does not have two elements in pixel size tag" );
- }
- }
- else
- throw Exception( "DICOM file does not specify pixel size" );
-
- // get slice spacing from multi-slice images.
- if ( ! this->Document().getValue( DCM_SpacingBetweenSlices, pixelSize[2] ) )
- {
- pixelSize[2] = 0;
- }
-
- return pixelSize;
-}
-
-const FixedVector<3,double>
-DICOM::GetImageOrigin() const
-{
- FixedVector<3,double> imageOrigin( 0.0 );
-
- const char *image_position_s = NULL;
- if ( ! this->Document().getValue( DCM_ImagePositionPatient, image_position_s ) )
- {
- // ImagePositionPatient tag not present, try ImagePosition instead
-#ifdef DCM_ImagePosition
- if ( ! this->Document().getValue( DCM_ImagePosition, image_position_s ) )
- image_position_s = NULL;
-#else
- if ( ! this->Document().getValue( DCM_ACR_NEMA_ImagePosition, image_position_s ) )
- image_position_s = NULL;
-#endif
- }
-
- if ( image_position_s )
- {
- double xyz[3];
- if ( 3 == sscanf( image_position_s,"%20lf%*c%20lf%*c%20lf", xyz, xyz+1, xyz+2 ) )
- {
- imageOrigin = FixedVector<3,double>::FromPointer( xyz );
- }
- }
-
- return imageOrigin;
-}
-
-const FixedArray< 2, FixedVector<3,double> >
-DICOM::GetImageOrientation() const
-{
- FixedArray< 2, FixedVector<3,double> > orientation;
-
- orientation[0] = FixedVector<3,double>( 0.0 );
- orientation[1] = FixedVector<3,double>( 0.0);
-
- orientation[0][0] = 1;
- orientation[1][1] = 1;
-
- const char *image_orientation_s = NULL;
-#ifdef DCM_ImageOrientation
- if ( ! this->Document().getValue( DCM_ImageOrientation, image_orientation_s ) )
-#else
- if ( ! this->Document().getValue( DCM_ACR_NEMA_ImageOrientation, image_orientation_s ) )
-#endif
- {
- // ImageOrientation tag not present, try ImageOrientationPatient instead
- if ( ! this->Document().getValue( DCM_ImageOrientationPatient, image_orientation_s ) )
- image_orientation_s = NULL;
- }
-
- if ( image_orientation_s )
- {
- double dx[3], dy[3];
- if ( 6 == sscanf( image_orientation_s, "%20lf%*c%20lf%*c%20lf%*c%20lf%*c%20lf%*c%20lf", dx, dx+1, dx+2, dy, dy+1, dy+2 ) )
- {
- orientation[0] = ( FixedVector<3,double>::FromPointer( dx ) );
- orientation[1] = ( FixedVector<3,double>::FromPointer( dy ) );
- }
- }
-
- return orientation;
-}
-
-TypedArray::SmartPtr
-DICOM::GetPixelDataArray( const size_t pixelDataLength )
-{
- DcmElement *delem = NULL;
-
- unsigned short bitsAllocated = 0;
- if ( ( delem = this->m_Document->search( DCM_BitsAllocated ) ) )
- {
- delem->getUint16( bitsAllocated );
- }
- else
- {
- // No "BitsAllocated" tag; use "BitsStored" instead.
- if ( ( delem = this->m_Document->search( DCM_BitsStored ) ) )
- {
- delem->getUint16( bitsAllocated );
- }
- }
-
- bool pixelDataSigned = false;
- Uint16 pixelRepresentation = 0;
- if ( this->m_Document->getValue( DCM_PixelRepresentation, pixelRepresentation ) > 0)
- pixelDataSigned = (pixelRepresentation == 1);
-
- double rescaleIntercept, rescaleSlope;
- const bool haveRescaleIntercept = (0 != this->m_Document->getValue( DCM_RescaleIntercept, rescaleIntercept ));
- if ( ! haveRescaleIntercept )
- rescaleIntercept = 0;
-
- const bool haveRescaleSlope = (0 != this->m_Document->getValue( DCM_RescaleSlope, rescaleSlope ));
- if ( ! haveRescaleSlope )
- rescaleSlope = 1;
-
- pixelDataSigned = pixelDataSigned || (rescaleIntercept < 0);
-
- Uint16 paddingValue = 0;
- const bool paddingFlag = (this->m_Dataset->findAndGetUint16( DCM_PixelPaddingValue, paddingValue )).good();
-
- TypedArray::SmartPtr pixelDataArray;
-
-#ifdef DCM_VariablePixelData
- delem = this->m_Document->search( DCM_VariablePixelData );
-#else
- delem = this->m_Document->search( DCM_ACR_NEMA_2C_VariablePixelData );
-#endif
- if (!delem)
- delem = this->m_Document->search( DCM_PixelData );
-
- if (delem)
- {
- if ( (delem->getTag().getEVR() == EVR_OW) || (bitsAllocated > 8) )
- {
- Uint16 *pdata = NULL;
- delem->getUint16Array(pdata);
- if ( pixelDataSigned )
- {
- const short paddingShort = static_cast( paddingValue );
- pixelDataArray = TypedArray::Create( TYPE_SHORT, pdata, pixelDataLength, paddingFlag, &paddingShort, Memory::ArrayCXX::DeleteWrapper );
- }
- else
- {
- const unsigned short paddingUShort = static_cast( paddingValue );
- pixelDataArray = TypedArray::Create( TYPE_USHORT, pdata, pixelDataLength, paddingFlag, &paddingUShort, Memory::ArrayCXX::DeleteWrapper );
- }
- }
- else
- {
- Uint8 *pdata = NULL;
- delem->getUint8Array(pdata);
- if ( pixelDataSigned )
- {
- const char paddingChar = static_cast( paddingValue );
- pixelDataArray = TypedArray::Create( TYPE_CHAR, pdata, pixelDataLength, paddingFlag, &paddingChar, Memory::ArrayCXX::DeleteWrapper );
- }
- else
- {
- const char paddingByte = static_cast( paddingValue );
- pixelDataArray = TypedArray::Create( TYPE_BYTE, pdata, pixelDataLength, paddingFlag, &paddingByte, Memory::ArrayCXX::DeleteWrapper );
- }
- }
-
- delem->detachValueField();
- }
-
- if ( ! pixelDataArray )
- {
- throw Exception( "Could not read pixel data from DICOM file" );
- }
-
- if ( haveRescaleIntercept || haveRescaleSlope )
- {
- double intpart = 0;
- if ( fabs(modf(rescaleSlope, &intpart) / rescaleSlope) > 1e-5 )
- {
- pixelDataArray = pixelDataArray->Convert(TYPE_FLOAT);
- }
-
- pixelDataArray->Rescale( rescaleSlope, rescaleIntercept );
- }
-
- return pixelDataArray;
-}
-
-
-const FixedVector<3,double>
-DICOM::DemosaicAndGetNormal
-( const FixedArray< 2, FixedVector<3,double> >& imageOrientation, const FixedVector<3,Types::Coordinate>& deltas, FixedVector<3,int>& dims, TypedArray::SmartPtr& pixelDataArray, FixedVector<3,double>& imageOrigin )
-{
- // without further information, we "guess" the image normal vector
- FixedVector<3,double> sliceNormal = SurfaceNormal( imageOrientation[0], imageOrientation[1] ).Get();
-
- // detect and treat Siemens multi-slice mosaics
- const char* tmpStr = NULL;
- if ( this->Document().getValue( DCM_Manufacturer, tmpStr ) )
- {
- if ( !strncmp( tmpStr, "SIEMENS", 7 ) )
- {
- Uint16 tempUint16 = 0;
- const DcmTagKey nSlicesTag(0x0019,0x100a);
- if ( this->Document().getValue( nSlicesTag, tempUint16 ) )
- {
- dims[2] = tempUint16;
- }
-
- // check for mosaic
- if ( dims[2] || (this->Document().getValue( DCM_ImageType, tmpStr ) && strstr( tmpStr, "MOSAIC" )) )
- {
- int unmosaicImageRows;
- int unmosaicImageCols;
-
- const DcmTagKey mosaicTag(0x0051,0x100b);
- if ( this->Document().getValue( mosaicTag, tmpStr ) )
- {
- if ( 2 != sscanf( tmpStr, "%6dp*%6ds", &unmosaicImageRows, &unmosaicImageCols) )
- {
- if ( 2 != sscanf( tmpStr, "%6d*%6ds", &unmosaicImageRows, &unmosaicImageCols) )
- {
- StdErr << "ERROR: unable to parse mosaic size from (0x0051,0x100b): " << tmpStr << "\n";
- }
- }
- }
-
-// For the following, see here: http://nipy.sourceforge.net/nibabel/dicom/siemens_csa.html#csa-header
- this->ParseSiemensCSA( DcmTagKey(0x0029,0x1020), unmosaicImageCols, unmosaicImageRows, dims[2], sliceNormal, imageOrigin ); // series information
- this->ParseSiemensCSA( DcmTagKey(0x0029,0x1010), unmosaicImageCols, unmosaicImageRows, dims[2], sliceNormal, imageOrigin ); // image information
-
- // hopefully we have figured out the mosaic dimensions by now.
- if ( (unmosaicImageCols > 0) && (unmosaicImageRows > 0 ) )
- {
- const int xMosaic = dims[0] / unmosaicImageCols;
-
- dims[0] = unmosaicImageCols;
- dims[1] = unmosaicImageRows;
-
- // de-mosaic the data array
- const unsigned long imageSizePixels = dims[0] * dims[1] * dims[2];
- TypedArray::SmartPtr newDataArray( TypedArray::Create( pixelDataArray->GetType(), imageSizePixels ) );
-
- const size_t pixelsPerSlice = unmosaicImageCols * unmosaicImageRows;
- size_t toOffset = 0;
- for ( int slice = 0; slice < dims[2]; ++slice )
- {
- for ( int j = 0; j < unmosaicImageRows; ++j, toOffset += dims[0] )
- {
- const size_t iPatch = slice % xMosaic;
- const size_t jPatch = slice / xMosaic;
-
- const size_t fromOffset = jPatch * xMosaic * pixelsPerSlice + j * xMosaic * unmosaicImageCols + iPatch * unmosaicImageCols;
- pixelDataArray->BlockCopy( *newDataArray, toOffset, fromOffset, unmosaicImageCols );
- }
- }
-
- pixelDataArray = newDataArray;
-
- // convert CSA-header center-of-image origin to corner-of-image standard origin (Issue #6754)
- imageOrigin -= (0.5 * ((dims[0]-1) * deltas[0] * imageOrientation[0] + (dims[1]-1) * deltas[1] * imageOrientation[1]) );
- }
- }
- }
- }
- return sliceNormal;
-}
-
-void
-DICOM::ParseSiemensCSA( const DcmTagKey& tagKey, int& unmosaicImageCols, int& unmosaicImageRows, int& slices, FixedVector<3,double>& sliceNormal, FixedVector<3,double>& imageOrigin )
-{
- const Uint8* csaHeaderInfo = NULL;
- unsigned long csaHeaderLength = 0;
- if ( this->Dataset().findAndGetUint8Array ( tagKey, csaHeaderInfo, &csaHeaderLength ).status() == OF_ok )
- {
- SiemensCSAHeader csaHeader( (const char*)csaHeaderInfo, csaHeaderLength );
-
- SiemensCSAHeader::const_iterator it = csaHeader.find( "AcquisitionMatrixText" );
- if ( (it != csaHeader.end()) && !it->second.empty() )
- {
- if ( 2 != sscanf( it->second[0].c_str(), "%6dp*%6ds", &unmosaicImageRows, &unmosaicImageCols) )
- {
- if ( 2 != sscanf( it->second[0].c_str(), "%6d*%6ds", &unmosaicImageRows, &unmosaicImageCols) )
- {
- StdErr << "ERROR: unable to parse mosaic size from CSA field AcquisitionMatrixText: " << it->second[0] << " in file " << this->m_Path << "\n";
- }
- }
- }
-
- it = csaHeader.find( "NumberOfImagesInMosaic" );
- if ( (it != csaHeader.end()) && !it->second.empty() )
- slices = atof( it->second[0].c_str() );
-
- it = csaHeader.find( "SliceNormalVector" );
- if ( (it != csaHeader.end()) && (it->second.size() >= 3) )
- {
- for ( size_t i = 0; i < 3; ++i )
- sliceNormal[i] = atof( it->second[i].c_str() );
- }
-
- // get true slice0 location from CSA header
- it = csaHeader.find( "MrPhoenixProtocol" );
- if ( (it != csaHeader.end()) && !it->second.empty() )
- {
- // ID strings for three axes in CSA header
- const std::string sliceOrientationString[] = { "dSag", "dCor", "dTra" };
-
- for ( int i = 0; i < 3; ++i )
- {
- const size_t sliceTagPos = it->second[0].find( std::string( "sSliceArray.asSlice[0].sPosition." ) + sliceOrientationString[i] );
- if ( sliceTagPos != std::string::npos )
- {
- const size_t equalPos = it->second[0].find( '=', sliceTagPos );
- if ( equalPos != std::string::npos )
- {
- imageOrigin[i] = atof( it->second[0].substr( equalPos + 1 ).c_str() );
- }
- else
- {
- StdErr << "ERROR: unable to get image origin component from: " << it->second[0] << " in file " << this->m_Path << "\nAssuming zero.\n";
- imageOrigin[i] = 0;
- }
- }
- else
- {
- StdErr << "ERROR: unable to get image origin tag for component " << sliceOrientationString[i] << " from CSA header in file " << this->m_Path << "\nAssuming zero.\n";
- imageOrigin[i] = 0;
- }
- }
- }
- }
-}
-
-ScalarImage*
-DICOM::Read
-( const char *path )
-{
- ScalarImage* image = NULL;
-
- Self dicom( path );
-
- FixedVector<3,int> dims = dicom.GetDims();
- FixedVector<3,double> pixelSize = dicom.GetPixelSize();
- ScalarImage::SpaceVectorType imageOrigin = dicom.GetImageOrigin();
-
- image = new ScalarImage( dims[0], dims[1], dims[2] );
- image->SetPixelSize( pixelSize[0], pixelSize[1] );
- image->SetFrameToFrameSpacing( pixelSize[2] );
-
- TypedArray::SmartPtr pixelDataArray = dicom.GetPixelDataArray( dims[0] * dims[1] * dims[2] );
- image->SetPixelData( pixelDataArray );
-
- // now some more manual readings...
-
- // get original table position from image.
- double sliceLocation = 0;
- if ( ! dicom.Document().getValue( DCM_SliceLocation, sliceLocation ) )
- {
-#ifdef DCM_Location
- dicom.Document().getValue( DCM_Location, sliceLocation );
-#else
- dicom.Document().getValue( DCM_ACR_NEMA_Location, sliceLocation );
-#endif
- }
- image->SetImageSlicePosition( sliceLocation );
- image->SetImageOrigin( imageOrigin );
-
- // get original image direction from file.
- FixedArray< 2, FixedVector<3,double> > imageOrientation = dicom.GetImageOrientation();
- image->SetImageDirectionX( imageOrientation[0] );
- image->SetImageDirectionY( imageOrientation[1] );
-
- return image;
-}
-
-//@}
-
-} // namespace cmtk
+/*
+//
+// Copyright 2014-2015 Google Inc.
+//
+// Copyright 2004-2014 SRI International
+//
+// Copyright 1997-2009 Torsten Rohlfing
+//
+// This file is part of the Computational Morphometry Toolkit.
+//
+// http://www.nitrc.org/projects/cmtk/
+//
+// The Computational Morphometry Toolkit is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// The Computational Morphometry Toolkit is distributed in the hope that it
+// will be useful, but WITHOUT ANY WARRANTY; without even the implied
+// warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with the Computational Morphometry Toolkit. If not, see
+// .
+//
+// $Revision$
+//
+// $LastChangedDate$
+//
+// $LastChangedBy$
+//
+*/
+
+#include "cmtkDICOM.h"
+
+#include
+#include
+
+#include
+#include
+
+#include
+
+
+#ifdef CMTK_USE_DCMTK_JPEG
+# include
+#endif
+
+#ifdef HAVE_SYS_TYPES_H
+# include
+#endif
+
+#include
+#include
+#include
+
+#ifndef DCM_ACR_NEMA_ImagePosition
+# define DCM_ACR_NEMA_ImagePosition DcmTagKey(0x0020, 0x0030)
+#endif
+
+#ifndef DCM_ACR_NEMA_ImageOrientation
+# define DCM_ACR_NEMA_ImageOrientation DcmTagKey(0x0020, 0x0035)
+#endif
+
+#ifndef DCM_ACR_NEMA_Location
+# define DCM_ACR_NEMA_Location DcmTagKey(0x0020, 0x0050)
+#endif
+
+#ifndef DCM_ACR_NEMA_2C_VariablePixelData
+# define DCM_ACR_NEMA_2C_VariablePixelData DcmTagKey(0x7f00, 0x0010)
+#endif
+
+namespace
+cmtk
+{
+
+/** \addtogroup IO */
+//@{
+
+void
+DICOM::InitFromFile( const std::string& path )
+{
+ this->m_Path = path;
+
+#ifdef CMTK_USE_DCMTK_JPEG
+ // register global decompression codecs
+ static bool decodersRegistered = false;
+ if ( ! decodersRegistered )
+ {
+ DJDecoderRegistration::registerCodecs( EDC_photometricInterpretation, EUC_default, EPC_default, 1 );
+ decodersRegistered = true;
+ }
+#endif
+
+ std::auto_ptr fileformat( new DcmFileFormat );
+ if (!fileformat.get())
+ {
+ throw Exception( "Could not create DICOM file format object." );
+ }
+
+ OFCondition status = fileformat->loadFile( path.c_str() );
+
+ if ( !status.good() )
+ {
+ throw Exception( "Cannot read DICOM file.." );
+// StdErr << "Error: cannot read DICOM file " << path << " (" << status.text() << ")\n";
+ }
+
+ this->m_Dataset = fileformat->getAndRemoveDataset();
+
+ if ( !this->m_Dataset )
+ {
+ throw Exception( "File format has NULL dataset." );
+ }
+
+ this->m_Document = std::auto_ptr( new DiDocument( this->m_Dataset, this->m_Dataset->getOriginalXfer(), CIF_AcrNemaCompatibility ) );
+ if ( ! this->m_Document.get() || ! this->m_Document->good() )
+ {
+ throw Exception( "Could not create document representation." );
+ }
+}
+
+const FixedVector<3,int>
+DICOM::GetDims() const
+{
+ FixedVector<3,int> dims( 0 );
+
+ Uint16 tempUint16 = 1;
+ if ( this->Document().getValue( DCM_Columns, tempUint16 ) )
+ {
+ dims[0] = static_cast( tempUint16 );
+ }
+
+ if ( this->Document().getValue( DCM_Rows, tempUint16 ) )
+ {
+ dims[1] = static_cast( tempUint16 );
+ }
+
+ // detect and treat multi-frame files
+ if ( ! this->Document().getValue( DCM_NumberOfFrames, tempUint16 ) )
+ {
+ // unlike Rows/Columns, NumberofFrames defaults to 1
+ tempUint16 = 1;
+ }
+ dims[2] = tempUint16;
+
+ return dims;
+}
+
+const FixedVector<3,double>
+DICOM::GetPixelSize() const
+{
+ FixedVector<3,double> pixelSize( 0.0 );
+
+ // get calibration from image
+ const bool hasPixelSpacing = (this->Document().getValue(DCM_PixelSpacing, pixelSize[0], 0) > 0);
+ if ( hasPixelSpacing )
+ {
+ if (this->Document().getValue(DCM_PixelSpacing, pixelSize[1], 1) < 2)
+ {
+ throw Exception( "DICOM file does not have two elements in pixel size tag" );
+ }
+ }
+ else
+ throw Exception( "DICOM file does not specify pixel size" );
+
+ // get slice spacing from multi-slice images.
+ if ( ! this->Document().getValue( DCM_SpacingBetweenSlices, pixelSize[2] ) )
+ {
+ pixelSize[2] = 0;
+ }
+
+ return pixelSize;
+}
+
+const FixedVector<3,double>
+DICOM::GetImageOrigin() const
+{
+ FixedVector<3,double> imageOrigin( 0.0 );
+
+ const char *image_position_s = NULL;
+ if ( ! this->Document().getValue( DCM_ImagePositionPatient, image_position_s ) )
+ {
+ // ImagePositionPatient tag not present, try ImagePosition instead
+#ifdef DCM_ImagePosition
+ if ( ! this->Document().getValue( DCM_ImagePosition, image_position_s ) )
+ image_position_s = NULL;
+#else
+ if ( ! this->Document().getValue( DCM_ACR_NEMA_ImagePosition, image_position_s ) )
+ image_position_s = NULL;
+#endif
+ }
+
+ if ( image_position_s )
+ {
+ double xyz[3];
+ if ( 3 == sscanf( image_position_s,"%20lf%*c%20lf%*c%20lf", xyz, xyz+1, xyz+2 ) )
+ {
+ imageOrigin = FixedVector<3,double>::FromPointer( xyz );
+ }
+ }
+
+ return imageOrigin;
+}
+
+const FixedArray< 2, FixedVector<3,double> >
+DICOM::GetImageOrientation() const
+{
+ FixedArray< 2, FixedVector<3,double> > orientation;
+
+ orientation[0] = FixedVector<3,double>( 0.0 );
+ orientation[1] = FixedVector<3,double>( 0.0);
+
+ orientation[0][0] = 1;
+ orientation[1][1] = 1;
+
+ const char *image_orientation_s = NULL;
+#ifdef DCM_ImageOrientation
+ if ( ! this->Document().getValue( DCM_ImageOrientation, image_orientation_s ) )
+#else
+ if ( ! this->Document().getValue( DCM_ACR_NEMA_ImageOrientation, image_orientation_s ) )
+#endif
+ {
+ // ImageOrientation tag not present, try ImageOrientationPatient instead
+ if ( ! this->Document().getValue( DCM_ImageOrientationPatient, image_orientation_s ) )
+ image_orientation_s = NULL;
+ }
+
+ if ( image_orientation_s )
+ {
+ double dx[3], dy[3];
+ if ( 6 == sscanf( image_orientation_s, "%20lf%*c%20lf%*c%20lf%*c%20lf%*c%20lf%*c%20lf", dx, dx+1, dx+2, dy, dy+1, dy+2 ) )
+ {
+ orientation[0] = ( FixedVector<3,double>::FromPointer( dx ) );
+ orientation[1] = ( FixedVector<3,double>::FromPointer( dy ) );
+ }
+ }
+
+ return orientation;
+}
+
+TypedArray::SmartPtr
+DICOM::GetPixelDataArray( const size_t pixelDataLength )
+{
+ DcmElement *delem = NULL;
+
+ unsigned short bitsAllocated = 0;
+ if ( ( delem = this->m_Document->search( DCM_BitsAllocated ) ) )
+ {
+ delem->getUint16( bitsAllocated );
+ }
+ else
+ {
+ // No "BitsAllocated" tag; use "BitsStored" instead.
+ if ( ( delem = this->m_Document->search( DCM_BitsStored ) ) )
+ {
+ delem->getUint16( bitsAllocated );
+ }
+ }
+
+ bool pixelDataSigned = false;
+ Uint16 pixelRepresentation = 0;
+ if ( this->m_Document->getValue( DCM_PixelRepresentation, pixelRepresentation ) > 0)
+ pixelDataSigned = (pixelRepresentation == 1);
+
+ double rescaleIntercept, rescaleSlope;
+ const bool haveRescaleIntercept = (0 != this->m_Document->getValue( DCM_RescaleIntercept, rescaleIntercept ));
+ if ( ! haveRescaleIntercept )
+ rescaleIntercept = 0;
+
+ const bool haveRescaleSlope = (0 != this->m_Document->getValue( DCM_RescaleSlope, rescaleSlope ));
+ if ( ! haveRescaleSlope )
+ rescaleSlope = 1;
+
+ pixelDataSigned = pixelDataSigned || (rescaleIntercept < 0);
+
+ Uint16 paddingValue = 0;
+ const bool paddingFlag = (this->m_Dataset->findAndGetUint16( DCM_PixelPaddingValue, paddingValue )).good();
+
+ TypedArray::SmartPtr pixelDataArray;
+
+#ifdef DCM_VariablePixelData
+ delem = this->m_Document->search( DCM_VariablePixelData );
+#else
+ delem = this->m_Document->search( DCM_ACR_NEMA_2C_VariablePixelData );
+#endif
+ if (!delem)
+ delem = this->m_Document->search( DCM_PixelData );
+
+ if (delem)
+ {
+ if ( (delem->getTag().getEVR() == EVR_OW) || (bitsAllocated > 8) )
+ {
+ Uint16 *pdata = NULL;
+ delem->getUint16Array(pdata);
+ if ( pixelDataSigned )
+ {
+ const short paddingShort = static_cast( paddingValue );
+ pixelDataArray = TypedArray::Create( TYPE_SHORT, pdata, pixelDataLength, paddingFlag, &paddingShort, Memory::ArrayCXX::DeleteWrapper );
+ }
+ else
+ {
+ const unsigned short paddingUShort = static_cast( paddingValue );
+ pixelDataArray = TypedArray::Create( TYPE_USHORT, pdata, pixelDataLength, paddingFlag, &paddingUShort, Memory::ArrayCXX::DeleteWrapper );
+ }
+ }
+ else
+ {
+ Uint8 *pdata = NULL;
+ delem->getUint8Array(pdata);
+ if ( pixelDataSigned )
+ {
+ const char paddingChar = static_cast( paddingValue );
+ pixelDataArray = TypedArray::Create( TYPE_CHAR, pdata, pixelDataLength, paddingFlag, &paddingChar, Memory::ArrayCXX::DeleteWrapper );
+ }
+ else
+ {
+ const char paddingByte = static_cast( paddingValue );
+ pixelDataArray = TypedArray::Create( TYPE_BYTE, pdata, pixelDataLength, paddingFlag, &paddingByte, Memory::ArrayCXX::DeleteWrapper );
+ }
+ }
+
+ delem->detachValueField();
+ }
+
+ if ( ! pixelDataArray )
+ {
+ throw Exception( "Could not read pixel data from DICOM file" );
+ }
+
+ if ( haveRescaleIntercept || haveRescaleSlope )
+ {
+ double intpart = 0;
+ if ( fabs(modf(rescaleSlope, &intpart) / rescaleSlope) > 1e-5 )
+ {
+ pixelDataArray = pixelDataArray->Convert(TYPE_FLOAT);
+ }
+
+ pixelDataArray->Rescale( rescaleSlope, rescaleIntercept );
+ }
+
+ return pixelDataArray;
+}
+
+
+inline const FixedVector<3,double>
+CornerToCenterTranslation
+(const FixedArray< 2, FixedVector<3,double> >& imageOrientation, const FixedVector<3,Types::Coordinate>& deltas, const FixedVector<3,int>& dims)
+{
+ return 0.5 * ((dims[0]-1) * deltas[0] * imageOrientation[0] + (dims[1]-1) * deltas[1] * imageOrientation[1]);
+}
+
+
+const FixedVector<3,double>
+DICOM::DemosaicAndGetNormal
+( const FixedArray< 2, FixedVector<3,double> >& imageOrientation, const FixedVector<3,Types::Coordinate>& deltas, FixedVector<3,int>& dims, TypedArray::SmartPtr& pixelDataArray, FixedVector<3,double>& imageOrigin )
+{
+ // without further information, we "guess" the image normal vector
+ FixedVector<3,double> sliceNormal = SurfaceNormal( imageOrientation[0], imageOrientation[1] ).Get();
+
+ // detect and treat Siemens multi-slice mosaics
+ const char* tmpStr = NULL;
+ if ( this->Document().getValue( DCM_Manufacturer, tmpStr ) )
+ {
+ if ( !strncmp( tmpStr, "SIEMENS", 7 ) )
+ {
+ Uint16 tempUint16 = 0;
+ const DcmTagKey nSlicesTag(0x0019,0x100a);
+ if ( this->Document().getValue( nSlicesTag, tempUint16 ) )
+ {
+ dims[2] = tempUint16;
+ }
+
+ // check for mosaic
+ if ( dims[2] || (this->Document().getValue( DCM_ImageType, tmpStr ) && strstr( tmpStr, "MOSAIC" )) )
+ {
+ int unmosaicImageRows;
+ int unmosaicImageCols;
+
+ const DcmTagKey mosaicTag(0x0051,0x100b);
+ if ( this->Document().getValue( mosaicTag, tmpStr ) )
+ {
+ if ( 2 != sscanf( tmpStr, "%6dp*%6ds", &unmosaicImageRows, &unmosaicImageCols) )
+ {
+ if ( 2 != sscanf( tmpStr, "%6d*%6ds", &unmosaicImageRows, &unmosaicImageCols) )
+ {
+ StdErr << "ERROR: unable to parse mosaic size from (0x0051,0x100b): " << tmpStr << "\n";
+ }
+ }
+ }
+
+ // For the following, see here: https://github.com/nipy/nibabel/blob/master/doc/source/dicom/siemens_csa.rst
+ this->ParseSiemensCSA( DcmTagKey(0x0029,0x1020), unmosaicImageCols, unmosaicImageRows, dims[2], sliceNormal); // series information
+ this->ParseSiemensCSA( DcmTagKey(0x0029,0x1010), unmosaicImageCols, unmosaicImageRows, dims[2], sliceNormal); // image information
+
+ // hopefully we have figured out the mosaic dimensions by now.
+ if ( (unmosaicImageCols > 0) && (unmosaicImageRows > 0 ) )
+ {
+ const int xMosaic = dims[0] / unmosaicImageCols;
+
+ // For reasoning, see the folllowing https://github.com/nipy/nibabel/blob/master/doc/source/dicom/dicom_mosaic.rst
+ imageOrigin += CornerToCenterTranslation(imageOrientation, deltas, dims);
+
+ dims[0] = unmosaicImageCols;
+ dims[1] = unmosaicImageRows;
+
+ // de-mosaic the data array
+ const unsigned long imageSizePixels = dims[0] * dims[1] * dims[2];
+ TypedArray::SmartPtr newDataArray( TypedArray::Create( pixelDataArray->GetType(), imageSizePixels ) );
+
+ const size_t pixelsPerSlice = unmosaicImageCols * unmosaicImageRows;
+ size_t toOffset = 0;
+ for ( int slice = 0; slice < dims[2]; ++slice )
+ {
+ for ( int j = 0; j < unmosaicImageRows; ++j, toOffset += dims[0] )
+ {
+ const size_t iPatch = slice % xMosaic;
+ const size_t jPatch = slice / xMosaic;
+
+ const size_t fromOffset = jPatch * xMosaic * pixelsPerSlice + j * xMosaic * unmosaicImageCols + iPatch * unmosaicImageCols;
+ pixelDataArray->BlockCopy( *newDataArray, toOffset, fromOffset, unmosaicImageCols );
+ }
+ }
+
+ pixelDataArray = newDataArray;
+
+ // convert CSA-header center-of-image origin to corner-of-image standard origin (Issue #6754)
+ imageOrigin -= CornerToCenterTranslation(imageOrientation, deltas, dims);
+ }
+ }
+ }
+ }
+ return sliceNormal;
+}
+
+void
+DICOM::ParseSiemensCSA( const DcmTagKey& tagKey, int& unmosaicImageCols, int& unmosaicImageRows, int& slices, FixedVector<3,double>& sliceNormal)
+{
+ const Uint8* csaHeaderInfo = NULL;
+ unsigned long csaHeaderLength = 0;
+ if ( this->Dataset().findAndGetUint8Array ( tagKey, csaHeaderInfo, &csaHeaderLength ).status() == OF_ok )
+ {
+ SiemensCSAHeader csaHeader( (const char*)csaHeaderInfo, csaHeaderLength );
+
+ SiemensCSAHeader::const_iterator it = csaHeader.find( "AcquisitionMatrixText" );
+ if ( (it != csaHeader.end()) && !it->second.empty() )
+ {
+ if ( 2 != sscanf( it->second[0].c_str(), "%6dp*%6ds", &unmosaicImageRows, &unmosaicImageCols) )
+ {
+ if ( 2 != sscanf( it->second[0].c_str(), "%6d*%6ds", &unmosaicImageRows, &unmosaicImageCols) )
+ {
+ StdErr << "ERROR: unable to parse mosaic size from CSA field AcquisitionMatrixText: " << it->second[0] << " in file " << this->m_Path << "\n";
+ }
+ }
+ }
+
+ it = csaHeader.find( "NumberOfImagesInMosaic" );
+ if ( (it != csaHeader.end()) && !it->second.empty() )
+ slices = atof( it->second[0].c_str() );
+
+ it = csaHeader.find( "SliceNormalVector" );
+ if ( (it != csaHeader.end()) && (it->second.size() >= 3) )
+ {
+ for ( size_t i = 0; i < 3; ++i )
+ sliceNormal[i] = atof( it->second[i].c_str() );
+ }
+ }
+}
+
+ScalarImage*
+DICOM::Read
+( const char *path )
+{
+ ScalarImage* image = NULL;
+
+ Self dicom( path );
+
+ FixedVector<3,int> dims = dicom.GetDims();
+ FixedVector<3,double> pixelSize = dicom.GetPixelSize();
+ ScalarImage::SpaceVectorType imageOrigin = dicom.GetImageOrigin();
+
+ image = new ScalarImage( dims[0], dims[1], dims[2] );
+ image->SetPixelSize( pixelSize[0], pixelSize[1] );
+ image->SetFrameToFrameSpacing( pixelSize[2] );
+
+ TypedArray::SmartPtr pixelDataArray = dicom.GetPixelDataArray( dims[0] * dims[1] * dims[2] );
+ image->SetPixelData( pixelDataArray );
+
+ // now some more manual readings...
+
+ // get original table position from image.
+ double sliceLocation = 0;
+ if ( ! dicom.Document().getValue( DCM_SliceLocation, sliceLocation ) )
+ {
+#ifdef DCM_Location
+ dicom.Document().getValue( DCM_Location, sliceLocation );
+#else
+ dicom.Document().getValue( DCM_ACR_NEMA_Location, sliceLocation );
+#endif
+ }
+ image->SetImageSlicePosition( sliceLocation );
+ image->SetImageOrigin( imageOrigin );
+
+ // get original image direction from file.
+ FixedArray< 2, FixedVector<3,double> > imageOrientation = dicom.GetImageOrientation();
+ image->SetImageDirectionX( imageOrientation[0] );
+ image->SetImageDirectionY( imageOrientation[1] );
+
+ return image;
+}
+
+//@}
+
+} // namespace cmtk
diff --git a/core/libs/IO/cmtkDICOM.h b/core/libs/IO/cmtkDICOM.h
index 9cb87a64..b61a89ac 100644
--- a/core/libs/IO/cmtkDICOM.h
+++ b/core/libs/IO/cmtkDICOM.h
@@ -1,8 +1,8 @@
-/*
-//
-// Copyright 2004-2012, 2013 SRI International
-//
-// Copyright 1997-2009 Torsten Rohlfing
+/*
+//
+// Copyright 2004-2012, 2013 SRI International
+//
+// Copyright 1997-2009 Torsten Rohlfing
//
// This file is part of the Computational Morphometry Toolkit.
//
@@ -10,13 +10,13 @@
//
// The Computational Morphometry Toolkit is free software: you can
// redistribute it and/or modify it under the terms of the GNU General Public
-// License as published by the Free Software Foundation, either version 3 of
-// the License, or (at your option) any later version.
-//
-// The Computational Morphometry Toolkit is distributed in the hope that it
-// will be useful, but WITHOUT ANY WARRANTY; without even the implied
-// warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-// GNU General Public License for more details.
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// The Computational Morphometry Toolkit is distributed in the hope that it
+// will be useful, but WITHOUT ANY WARRANTY; without even the implied
+// warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with the Computational Morphometry Toolkit. If not, see
@@ -30,123 +30,122 @@
//
*/
-#ifndef __cmtkDICOM_h_included_
-#define __cmtkDICOM_h_included_
-
-#include
-
-#include
-
-#include
-#include
-#include
-
-#include
-
-namespace
-cmtk
-{
-
-/** \addtogroup IO */
-//@{
-
-/** Reader/writer class for DICOM images.
- */
-class DICOM
-{
-public:
- /// This class.
- typedef DICOM Self;
-
- /// Default constructor.
- DICOM() : m_Dataset( NULL ) {}
-
- /// Read file constructor.
- DICOM( const std::string& path ) { this->InitFromFile( path ); }
-
- /// Read object from DICOM file.
- void InitFromFile( const std::string& path );
-
- /// Get image dimensions (number of pixels per axis).
- const FixedVector<3,int> GetDims() const;
-
- /// Get image pixel size.
- const FixedVector<3,double> GetPixelSize() const;
-
- /// Get image origin in scanner coordinates.
- const FixedVector<3,double> GetImageOrigin() const;
-
- /// Get image orientation of the i and j grid axes.
- const FixedArray< 2, FixedVector<3,double> > GetImageOrientation() const;
-
- /** Get pixel data array.
- * The pixel data type is determined automatically based on bits allocated and signed vs. unsigned representation.
- *
- * If the RescaleSlope or RescaleIntercept tags are present, intensity rescaling is applied.
- *
- * If a padding value is defined in the DICOM file, this value is also set as padding in the output array.
- *
- *\warning As a side effect, this function releases the pixel data array pointer from the DICOM object, i.e.,
- * this function can only be called ONCE for each object.
- */
- TypedArray::SmartPtr GetPixelDataArray( const size_t pixelDataLength );
-
- /// Demosaic (if necessary) the image and return slice normal vector.
- const FixedVector<3,double> DemosaicAndGetNormal( const FixedArray< 2, FixedVector<3,double> >& imageOrientation /*!< Image orientation vectors (read-only, for default normal computation)*/,
- const FixedVector<3,Types::Coordinate>& deltas /*!< Pixel size - this is needed for computing image origin coordinates from image center location */,
- FixedVector<3,int>& dims /*!< Image dimensions - these may be modified if the image is a mosaic */,
- TypedArray::SmartPtr& pixelDataArray /*!< Pixel data array - a new array will be returned if the image is a mosaic.*/,
- FixedVector<3,double>& imageOrigin /*!< True image origin from CSA header if this is a mosaic file. */);
-
- /// Get const DICOM dataset.
- const DcmDataset& Dataset() const
- {
- return *(this->m_Dataset);
- }
-
- /// Get DICOM dataset.
- DcmDataset& Dataset()
- {
- return *(this->m_Dataset);
- }
-
- /// Get const DICOM document.
- const DiDocument& Document() const
- {
- return *(this->m_Document);
- }
-
- /// Get DICOM document.
- DiDocument& Document()
- {
- return *(this->m_Document);
- }
-
- /** Read ScalarImage from DICOM file.
- */
- static ScalarImage* Read( const char *path );
-
-private:
- /// Path of this DICOM file.
- std::string m_Path;
-
- /// Pointer to the DICOM dataset object
- DcmDataset* m_Dataset;
-
- /// Pointer to the DICOM document object
- std::auto_ptr m_Document;
-
- /// Parse private Siemens CSA data.
- void ParseSiemensCSA( const DcmTagKey& tagKey /*!< DCM tag with the Siemens CSA header data */,
- int& unmosaicImageCols /*!< Return demosaiced image columns */,
- int& unmosaicImageRows /*!< Return demosaiced image rows */,
- int& slices /*!< Return demosaiced image slices */,
- FixedVector<3,double>& sliceNormal /*!< Return demosaiced image slice normal */,
- FixedVector<3,double>& imageOrigin /*!< Return demosaiced image origin */ );
-};
-
-//@}
-
-} // namespace cmtk
-
-#endif // #ifndef __cmtkDICOM_h_included_
+#ifndef __cmtkDICOM_h_included_
+#define __cmtkDICOM_h_included_
+
+#include
+
+#include
+
+#include
+#include
+#include
+
+#include
+
+namespace
+cmtk
+{
+
+/** \addtogroup IO */
+//@{
+
+/** Reader/writer class for DICOM images.
+ */
+class DICOM
+{
+public:
+ /// This class.
+ typedef DICOM Self;
+
+ /// Default constructor.
+ DICOM() : m_Dataset( NULL ) {}
+
+ /// Read file constructor.
+ DICOM( const std::string& path ) { this->InitFromFile( path ); }
+
+ /// Read object from DICOM file.
+ void InitFromFile( const std::string& path );
+
+ /// Get image dimensions (number of pixels per axis).
+ const FixedVector<3,int> GetDims() const;
+
+ /// Get image pixel size.
+ const FixedVector<3,double> GetPixelSize() const;
+
+ /// Get image origin in scanner coordinates.
+ const FixedVector<3,double> GetImageOrigin() const;
+
+ /// Get image orientation of the i and j grid axes.
+ const FixedArray< 2, FixedVector<3,double> > GetImageOrientation() const;
+
+ /** Get pixel data array.
+ * The pixel data type is determined automatically based on bits allocated and signed vs. unsigned representation.
+ *
+ * If the RescaleSlope or RescaleIntercept tags are present, intensity rescaling is applied.
+ *
+ * If a padding value is defined in the DICOM file, this value is also set as padding in the output array.
+ *
+ *\warning As a side effect, this function releases the pixel data array pointer from the DICOM object, i.e.,
+ * this function can only be called ONCE for each object.
+ */
+ TypedArray::SmartPtr GetPixelDataArray( const size_t pixelDataLength );
+
+ /// Demosaic (if necessary) the image and return slice normal vector.
+ const FixedVector<3,double> DemosaicAndGetNormal( const FixedArray< 2, FixedVector<3,double> >& imageOrientation /*!< Image orientation vectors (read-only, for default normal computation)*/,
+ const FixedVector<3,Types::Coordinate>& deltas /*!< Pixel size - this is needed for computing image origin coordinates from image center location */,
+ FixedVector<3,int>& dims /*!< Image dimensions - these may be modified if the image is a mosaic */,
+ TypedArray::SmartPtr& pixelDataArray /*!< Pixel data array - a new array will be returned if the image is a mosaic.*/,
+ FixedVector<3,double>& imageOrigin /*!< True image origin from CSA header if this is a mosaic file. */);
+
+ /// Get const DICOM dataset.
+ const DcmDataset& Dataset() const
+ {
+ return *(this->m_Dataset);
+ }
+
+ /// Get DICOM dataset.
+ DcmDataset& Dataset()
+ {
+ return *(this->m_Dataset);
+ }
+
+ /// Get const DICOM document.
+ const DiDocument& Document() const
+ {
+ return *(this->m_Document);
+ }
+
+ /// Get DICOM document.
+ DiDocument& Document()
+ {
+ return *(this->m_Document);
+ }
+
+ /** Read ScalarImage from DICOM file.
+ */
+ static ScalarImage* Read( const char *path );
+
+private:
+ /// Path of this DICOM file.
+ std::string m_Path;
+
+ /// Pointer to the DICOM dataset object
+ DcmDataset* m_Dataset;
+
+ /// Pointer to the DICOM document object
+ std::auto_ptr m_Document;
+
+ /// Parse private Siemens CSA data.
+ void ParseSiemensCSA( const DcmTagKey& tagKey /*!< DCM tag with the Siemens CSA header data */,
+ int& unmosaicImageCols /*!< Return demosaiced image columns */,
+ int& unmosaicImageRows /*!< Return demosaiced image rows */,
+ int& slices /*!< Return demosaiced image slices */,
+ FixedVector<3,double>& sliceNormal /*!< Return demosaiced image slice normal */);
+};
+
+//@}
+
+} // namespace cmtk
+
+#endif // #ifndef __cmtkDICOM_h_included_
diff --git a/core/libs/IO/cmtkImageFileDICOM.cxx b/core/libs/IO/cmtkImageFileDICOM.cxx
index 2ce810d4..8d4b143e 100644
--- a/core/libs/IO/cmtkImageFileDICOM.cxx
+++ b/core/libs/IO/cmtkImageFileDICOM.cxx
@@ -35,13 +35,13 @@
#include
#include
-#include
#include
#include
#include
#include
+#include
namespace cmtk
{
@@ -219,19 +219,19 @@ namespace cmtk
this->m_AcquisitionNumber = 0;
// check for which vendor and deal with specifics elsewhere
- if (this->m_TagToStringMap[DCM_Manufacturer].substr(0, 7) == "SIEMENS")
+ std::string vendor( this->m_TagToStringMap[DCM_Manufacturer] );
+ std::transform(vendor.begin(), vendor.end(), vendor.begin(), ::toupper);
+ if ( vendor.compare( 0, 7, "SIEMENS" ) == 0 )
{
this->DoVendorTagsSiemens();
}
-
- if (this->m_TagToStringMap[DCM_Manufacturer].substr(0, 2) == "GE")
- {
- this->DoVendorTagsGE();
- }
-
- if (this->m_TagToStringMap[DCM_Manufacturer].substr(0, 7) == "Philips")
+ else if ( vendor.compare( 0, 7, "PHILIPS" ) == 0 )
{
this->DoVendorTagsPhilips();
+ }
+ else if ( vendor.compare( 0, 2, "GE" ) == 0 )
+ {
+ this->DoVendorTagsGE();
}
}
@@ -289,12 +289,186 @@ namespace cmtk
this->m_TagToStringMap[DCM_AcquisitionMatrix] = join(vAcquisitionMatrix.begin(), vAcquisitionMatrix.end(), "\\");
}
+ bool
+ ImageFileDICOM::DoVendorTagsSiemensDwellTime(const SiemensCSAHeader* csaImageHeader)
+ {
+ // Get the dwell time from the dedicated "RealDwellTime" private
+ // SIEMENS DICOM attribute if it exists, or from the CSA header if it
+ // exists and contains it. The units of the SIEMENS dwell times are
+ // nano-seconds.
+ bool dwell_time_found = false;
+ double dwell_time_value;
+ unsigned long dwell_time_length = 0;
+ const Uint8* dwell_time_buffer = NULL;
+ OFString dwell_time_str;
+
+ if ( this->m_Dataset->findAndGetUint8Array( DCM_Siemens_RealDwellTime, dwell_time_buffer, &dwell_time_length ).good() )
+ {
+ if( dwell_time_found = ((dwell_time_buffer != NULL) && (dwell_time_length > 0)) )
+ {
+ dwell_time_str.assign( OFstatic_cast(const char*, dwell_time_buffer), dwell_time_length);
+ dwell_time_value = atof( dwell_time_str.c_str() );
+ }
+ }
+
+ if ( !dwell_time_found && this->m_Document->getValue(DCM_Siemens_RealDwellTime, dwell_time_str) > 0 )
+ {
+ if ( dwell_time_found = !dwell_time_str.empty() )
+ dwell_time_value = atof( dwell_time_str.c_str() );
+ }
+
+ if ( !dwell_time_found && csaImageHeader != NULL )
+ {
+ SiemensCSAHeader::const_iterator csa_it = csaImageHeader->find( "RealDwellTime" );
+ if ( dwell_time_found = ((csa_it != csaImageHeader->end()) && !csa_it->second.empty()) )
+ dwell_time_value = atof( csa_it->second[0].c_str() );
+ }
+ // I do not know why the reciprocal of the reported dwell time is being
+ // stored in the dwell time member. Does this not give us Bandwidth
+ // instead of dwell time?
+ if ( dwell_time_found )
+ {
+ if ( dwell_time_value != 0.0 )
+ this->m_DwellTime = 1.0 / dwell_time_value;
+ else
+ this->m_DwellTime = 0.0;
+ }
+ return dwell_time_found;
+ }
+
+ bool
+ ImageFileDICOM::DoVendorTagsSiemensPhaseEncodeDirection(const SiemensCSAHeader* csaImageHeader)
+ {
+ // Get the polarity of the phase encoding from the dedicated
+ // "PhaseEncodingDirectionPositive" private SIEMENS DICOM attribute if
+ // it exists, or from the CSA header if it exists and contains it
+ bool polarity_found = false;
+ char polarity_value = '\0';
+ unsigned long polarity_length = 0;
+ const Uint8* polarity_buffer = NULL;
+ OFString polarity_str;
+
+ if ( this->m_Dataset->findAndGetUint8Array( DCM_Siemens_PhaseEncodingDirectionPositive, polarity_buffer, &polarity_length ).good() )
+ {
+ if( polarity_found = (polarity_buffer != NULL && polarity_length > 0) )
+ polarity_value = *( OFstatic_cast(const char*, polarity_buffer) );
+ }
+
+ if ( !polarity_found && this->m_Document->getValue(DCM_Siemens_PhaseEncodingDirectionPositive, polarity_str) )
+ {
+ if ( polarity_found = !polarity_str.empty() )
+ polarity_value = polarity_str[0];
+ }
+
+ if ( !polarity_found && csaImageHeader != NULL )
+ {
+ SiemensCSAHeader::const_iterator csa_it = csaImageHeader->find( "PhaseEncodingDirectionPositive");
+ if ( polarity_found = (csa_it != csaImageHeader->end() && !csa_it->second.empty()) )
+ polarity_value = csa_it->second[0][0];
+
+ }
+ // use of the polarity_found flag here is to preserve the previous
+ // functionality where the m_PhaseEncodeDirectionSign member was
+ // only set when polarity_value was found.
+ if ( polarity_found )
+ {
+ switch ( polarity_value )
+ {
+ case '0':
+ this->m_PhaseEncodeDirectionSign = "NEG";
+ break;
+ case '1':
+ this->m_PhaseEncodeDirectionSign = "POS";
+ break;
+ default:
+ this->m_PhaseEncodeDirectionSign = "UNKNOWN";
+ break;
+ }
+ }
+ return polarity_found;
+ }
+
+ bool
+ ImageFileDICOM::DoVendorTagsSiemensDiffusionBValue(const SiemensCSAHeader* csaImageHeader)
+ {
+ unsigned long bvalue_length = 0;
+ const Uint8* bvalue_buffer = NULL;
+ OFString bvalue_str;
+
+ if ( this->m_Dataset->findAndGetUint8Array( DCM_Siemens_DiffusionBValue, bvalue_buffer, &bvalue_length ).good() )
+ {
+ if( bvalue_buffer != NULL && bvalue_length > 0 )
+ {
+ bvalue_str.assign( OFstatic_cast(const char*, bvalue_buffer), bvalue_length);
+ this->m_BValue = atof( bvalue_str.c_str() );
+ return true;
+ }
+ }
+
+ if ( this->m_Document->getValue(DCM_Siemens_DiffusionBValue, bvalue_str) > 0 )
+ {
+ if( !bvalue_str.empty() )
+ {
+ this->m_BValue = atof( bvalue_str.c_str() );
+ return true;
+ }
+ }
+
+ if ( csaImageHeader != NULL )
+ {
+ SiemensCSAHeader::const_iterator csa_it = csaImageHeader->find( "B_value" );
+ if ( (csa_it != csaImageHeader->end()) && !csa_it->second.empty() )
+ {
+ this->m_BValue = atof( csa_it->second[0].c_str() );
+ return true;
+ }
+ }
+ return false;
+ }
+
+ bool
+ ImageFileDICOM::DoVendorTagsSiemensDiffusionBVector(const SiemensCSAHeader* csaImageHeader)
+ {
+ unsigned long bvector_length = 0;
+ const Uint8* bvector_buffer = NULL;
+
+ if ( this->m_Dataset->findAndGetUint8Array( DCM_Siemens_DiffusionGradientOrientation, bvector_buffer, &bvector_length ).good() )
+ {
+ if ( bvector_buffer != NULL && (bvector_length * sizeof(Uint8)) >= (this->m_BVector.Size() * sizeof(double)) )
+ {
+ const double* input_itr = OFstatic_cast(const double*, bvector_buffer);
+ std::copy(input_itr, input_itr + this->m_BVector.Size(), this->m_BVector.begin());
+ return true;
+ }
+ }
+
+ if ( this->m_Document->getValue(DCM_Siemens_DiffusionGradientOrientation, this->m_BVector[0], 0) >= this->m_BVector.Size() )
+ {
+ for ( size_t idx = 1; idx < this->m_BVector.Size(); ++idx )
+ this->m_Document->getValue(DCM_Siemens_DiffusionGradientOrientation, this->m_BVector[idx], idx);
+ return true;
+ }
+
+ if ( csaImageHeader != NULL )
+ {
+ SiemensCSAHeader::const_iterator csa_it = csaImageHeader->find( "DiffusionGradientDirection" );
+ if ( (csa_it != csaImageHeader->end()) && (csa_it->second.size() >= this->m_BVector.Size()) )
+ {
+ for ( size_t idx = 0; idx < this->m_BVector.Size(); ++idx )
+ this->m_BVector[idx] = atof( csa_it->second[idx].c_str() );
+ return true;
+ }
+ }
+ return false;
+ }
+
void
ImageFileDICOM::DoVendorTagsSiemens()
{
Uint16 nFrames = 0;
double tmpDbl = 0;
const char *tmpStr = NULL;
+ SiemensCSAHeader::const_iterator csa_it;
this->m_IsMultislice = (0 != this->m_Document->getValue(DcmTagKey(0x0019, 0x100a), nFrames)); // Number of Slices tag
this->m_IsMultislice |= (this->m_Document->getValue(DCM_ImageType, tmpStr) && strstr(tmpStr, "MOSAIC")); // mosaics are always multi-slice
@@ -319,88 +493,60 @@ namespace cmtk
this->m_RawDataType = "real";
}
- // Get dwell time if cas header exists and contains it
- const Uint8 *csaImageHeaderInfo = NULL;
+ // If the CSA Header exists, instantiate SiemensCSAHeader to use in
+ // the retrieval of subsequent imagining quantities
+ const Uint8* csaImageHeaderInfo = NULL;
unsigned long csaImageHeaderLength = 0;
- if (this->m_Dataset->findAndGetUint8Array(DcmTagKey(0x0029, 0x1010), csaImageHeaderInfo, &csaImageHeaderLength).status() == OF_ok) // the "Image" CSA header, not the "Series" header.
+ std::auto_ptr csaImageHeader;
+ if ( this->m_Dataset->findAndGetUint8Array( DcmTagKey(0x0029,0x1010), csaImageHeaderInfo, &csaImageHeaderLength ).status() == OF_ok ) // the "Image" CSA header, not the "Series" header.
{
- SiemensCSAHeader csaImageHeader((const char *)csaImageHeaderInfo, csaImageHeaderLength);
- SiemensCSAHeader::const_iterator it = csaImageHeader.find("RealDwellTime");
- if ((it != csaImageHeader.end()) && !it->second.empty())
- {
- this->m_DwellTime = 1.0 / atof(it->second[0].c_str());
- }
- else
- {
- this->m_DwellTime = 0.0;
- }
+ csaImageHeader.reset( new SiemensCSAHeader( OFstatic_cast(const char*, csaImageHeaderInfo), csaImageHeaderLength ) );
+ }
- it = csaImageHeader.find("PhaseEncodingDirectionPositive");
- if ((it != csaImageHeader.end()) && !it->second.empty())
- {
- switch (it->second[0][0])
- {
- case '0':
- this->m_PhaseEncodeDirectionSign = "NEG";
- break;
- case '1':
- this->m_PhaseEncodeDirectionSign = "POS";
- break;
- default:
- this->m_PhaseEncodeDirectionSign = "UNKNOWN";
- break;
- }
- }
+ DoVendorTagsSiemensDwellTime(csaImageHeader.get());
+ DoVendorTagsSiemensPhaseEncodeDirection(csaImageHeader.get());
+
+ bool bvalue_found = DoVendorTagsSiemensDiffusionBValue(csaImageHeader.get());
+ // 20161028 djk: make m_IsDWI is always true if the B_value field is defined
+ this->m_IsDWI |= bvalue_found;
+
+ this->m_HasBVector = DoVendorTagsSiemensDiffusionBVector(csaImageHeader.get());
+ this->m_HasBVector |= bvalue_found;
+ this->m_IsDWI |= this->m_HasBVector;
+
+ bool found_directionality = false;
+ OFString directionality;
+ unsigned long directionality_length = 0;
+ const Uint8* directionality_buffer = NULL;
+ if ( this->m_Dataset->findAndGetUint8Array( DCM_Siemens_DiffusionDirectionality, directionality_buffer, &directionality_length ).good() )
+ {
+ if( found_directionality = (directionality_buffer != NULL && directionality_length > 0) )
+ directionality.assign( OFstatic_cast(const char *, directionality_buffer), directionality_length );
}
- // for DWI, first check standard DICOM vendor tags
- if ((this->m_IsDWI = (this->m_Document->getValue(DcmTagKey(0x0019, 0x100d), tmpStr) != 0))) // "Directionality" tag
+ if ( !found_directionality )
+ found_directionality = (this->m_Document->getValue(DCM_Siemens_DiffusionDirectionality, directionality) > 0);
+
+ if ( !found_directionality && csaImageHeader.get() != NULL )
{
- if (this->m_Document->getValue(DcmTagKey(0x0019, 0x100c), tmpStr) != 0) // bValue tag
- {
- this->m_BValue = atof(tmpStr);
- this->m_IsDWI |= (this->m_BValue > 0);
- }
+ csa_it = csaImageHeader->find( "DiffusionDirectionality" );
+ if ( found_directionality = (csa_it != csaImageHeader->end() && !csa_it->second.empty()) )
+ directionality = csa_it->second[0].c_str();
+ }
- if (this->m_BValue > 0)
+ if ( found_directionality )
+ {
+ if ( 0 == directionality.compare( 0, 11, "DIRECTIONAL" ) )
{
- for (int idx = 0; idx < 3; ++idx)
- {
- this->m_IsDWI |= (this->m_Document->getValue(DcmTagKey(0x0019, 0x100e), this->m_BVector[idx], idx) != 0);
- }
+ this->m_IsDWI = true;
+ this->m_HasBVector = true;
}
- }
- else
- {
- // no hint of DWI in standard tags, look into CSA header to confirm
- if (csaImageHeaderInfo)
+ else if ( 0 == directionality.compare( 0, 9, "ISOTROPIC" ) )
{
- SiemensCSAHeader csaImageHeader((const char *)csaImageHeaderInfo, csaImageHeaderLength);
- SiemensCSAHeader::const_iterator it = csaImageHeader.find("DiffusionDirectionality");
- if ((it != csaImageHeader.end()) && !it->second.empty())
- {
- this->m_IsDWI = (0 == it->second[0].compare(0, 11, "DIRECTIONAL"));
- }
-
- it = csaImageHeader.find("B_value");
- if ((it != csaImageHeader.end()) && !it->second.empty())
- {
- this->m_BValue = atof(it->second[0].c_str());
- // 20161028 djk: make m_IsDWI is always true if the B_value field is defined
- this->m_IsDWI = true;
- }
-
- it = csaImageHeader.find("DiffusionGradientDirection");
- if ((it != csaImageHeader.end()) && (it->second.size() >= 3))
- {
- for (int idx = 0; idx < 3; ++idx)
- {
- this->m_BVector[idx] = atof(it->second[idx].c_str());
- }
- }
+ this->m_IsDWI = true;
+ this->m_HasBVector = false;
}
}
- this->m_HasBVector = this->m_IsDWI;
}
}
diff --git a/core/libs/IO/cmtkImageFileDICOM.h b/core/libs/IO/cmtkImageFileDICOM.h
index 24699694..e3483831 100644
--- a/core/libs/IO/cmtkImageFileDICOM.h
+++ b/core/libs/IO/cmtkImageFileDICOM.h
@@ -40,6 +40,8 @@
#include
+#include
+
#include
#include
@@ -64,6 +66,13 @@
#define DCM_Siemens_MosaicRefAcqTimes DcmTagKey(0x0019,0x1029)
#endif
+// Siemens private DICOM attributes for Diffusion MR.
+#define DCM_Siemens_DiffusionDirectionality DcmTagKey(0x0019, 0x100d)
+#define DCM_Siemens_DiffusionBValue DcmTagKey(0x0019,0x100c)
+#define DCM_Siemens_DiffusionGradientOrientation DcmTagKey(0x0019,0x100e)
+#define DCM_Siemens_PhaseEncodingDirectionPositive DcmTagKey(0x0021,0x111c)
+#define DCM_Siemens_RealDwellTime DcmTagKey(0x0021,0x1142)
+
#include
#include