diff --git a/Modules/DICOMReader/include/mitkITKDICOMSeriesReaderHelper.h b/Modules/DICOMReader/include/mitkITKDICOMSeriesReaderHelper.h index 8792dce94f..4b40a4ffe6 100644 --- a/Modules/DICOMReader/include/mitkITKDICOMSeriesReaderHelper.h +++ b/Modules/DICOMReader/include/mitkITKDICOMSeriesReaderHelper.h @@ -1,101 +1,109 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkDICOMSeriesReaderHelper_h #define mitkDICOMSeriesReaderHelper_h #include "mitkImage.h" #include "mitkGantryTiltInformation.h" #include "mitkDICOMTag.h" #include /* Forward deceleration of an DCMTK class. Used in the txx but part of the interface.*/ class OFDateTime; namespace mitk { class ITKDICOMSeriesReaderHelper { public: static const DICOMTag AcquisitionDateTag; static const DICOMTag AcquisitionTimeTag; static const DICOMTag TriggerTimeTag; typedef std::vector StringContainer; typedef std::list StringContainerList; Image::Pointer Load( const StringContainer& filenames, bool correctTilt, const GantryTiltInformation& tiltInfo ); Image::Pointer Load3DnT( const StringContainerList& filenamesLists, bool correctTilt, const GantryTiltInformation& tiltInfo ); static bool CanHandleFile(const std::string& filename); private: typedef std::vector TimeBoundsList; typedef itk::FixedArray DateTimeBounds; /** Scans the given files for the acquisition time and returns the lowest and - highest acquisition time as time bounds via bounds. If no acquisition times can be found - the function return will be false. + highest acquisition date time as date time bounds via bounds. + @param bounds The acquisition date time bound extracted from the files. + @param triggerBounds Time bounds for trigger information extracted from the files. + If no trigger information was found than it returns trigger == [0.0, 0.0]. + @return If no acquisition date times can be found the function return will be false. Otherwise + it returns True. */ - static bool ExtractDateTimeBoundsOfTimeStep( const StringContainer& filenamesOfTimeStep, - DateTimeBounds& bounds ); + static bool ExtractDateTimeBoundsAndTriggerOfTimeStep( const StringContainer& filenamesOfTimeStep, + DateTimeBounds& bounds, TimeBounds& triggerBounds); - - static bool ExtractDateTimeBoundsOfTimeStep( const StringContainer& filenamesOfTimeStep, + /* Determine the time bounds in ms respective to the baselineDateTime for the passed + files. Additionaly it regards the trigger time tag if set and acquisition date time + carries not enough information.*/ + static bool ExtractTimeBoundsOfTimeStep(const StringContainer& filenamesOfTimeStep, TimeBounds& bounds, const OFDateTime& baselineDateTime ); /** Returns the list of time bounds of all passed time step containers. - (sa ExtracTimeBoundsOfTimeStep). Time steps where no time bounds could be extracted + (sa ExtractTimeBoundsOfTimeStep and ExtractDateTimeBoundsOfTimeStep). + Time steps where no time bounds could be extracted are indecated by "null" time bounds (both values "0"). The order of the returned - list equals of passed filenamesOfTimeSteps order.*/ + list equals of passed filenamesOfTimeSteps order. + @remark The function regards acquisition date time tags and trigger time tags.*/ static TimeBoundsList ExtractTimeBoundsOfTimeSteps (const StringContainerList& filenamesOfTimeSteps); /** Helper function that generates a time geometry using the template and the passed boundslist (which indicates the number of time steps). */ static TimeGeometry::Pointer GenerateTimeGeometry(const BaseGeometry* templateGeometry, const TimeBoundsList& boundsList); template typename ImageType::Pointer FixUpTiltedGeometry( ImageType* input, const GantryTiltInformation& tiltInfo ); template Image::Pointer LoadDICOMByITK( const StringContainer& filenames, bool correctTilt, const GantryTiltInformation& tiltInfo, itk::GDCMImageIO::Pointer& io); template Image::Pointer LoadDICOMByITK3DnT( const StringContainerList& filenames, bool correctTilt, const GantryTiltInformation& tiltInfo, itk::GDCMImageIO::Pointer& io); }; } #endif diff --git a/Modules/DICOMReader/src/mitkITKDICOMSeriesReaderHelper.cpp b/Modules/DICOMReader/src/mitkITKDICOMSeriesReaderHelper.cpp index fea2260e1a..2984ec7f19 100644 --- a/Modules/DICOMReader/src/mitkITKDICOMSeriesReaderHelper.cpp +++ b/Modules/DICOMReader/src/mitkITKDICOMSeriesReaderHelper.cpp @@ -1,492 +1,460 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ //#define MBILOG_ENABLE_DEBUG #include #define BOOST_DATE_TIME_NO_LIB //Prevent unnecessary/unwanted auto link in this compilation when activating boost libraries in the MITK superbuild //It is necessary because BOOST_ALL_DYN_LINK overwrites BOOST_DATE_TIME_NO_LIB #if defined(BOOST_ALL_DYN_LINK) #undef BOOST_ALL_DYN_LINK #endif #include #include "mitkITKDICOMSeriesReaderHelper.h" #include "mitkITKDICOMSeriesReaderHelper.txx" #include "mitkDICOMGDCMTagScanner.h" #include "mitkArbitraryTimeGeometry.h" #include "dcvrda.h" const mitk::DICOMTag mitk::ITKDICOMSeriesReaderHelper::AcquisitionDateTag = mitk::DICOMTag( 0x0008, 0x0022 ); const mitk::DICOMTag mitk::ITKDICOMSeriesReaderHelper::AcquisitionTimeTag = mitk::DICOMTag( 0x0008, 0x0032 ); const mitk::DICOMTag mitk::ITKDICOMSeriesReaderHelper::TriggerTimeTag = mitk::DICOMTag( 0x0018, 0x1060 ); #define switch3DCase( IOType, T ) \ case IOType: \ return LoadDICOMByITK( filenames, correctTilt, tiltInfo, io ); bool mitk::ITKDICOMSeriesReaderHelper::CanHandleFile( const std::string& filename ) { MITK_DEBUG << "ITKDICOMSeriesReaderHelper::CanHandleFile " << filename; itk::GDCMImageIO::Pointer tester = itk::GDCMImageIO::New(); return tester->CanReadFile( filename.c_str() ); } mitk::Image::Pointer mitk::ITKDICOMSeriesReaderHelper::Load( const StringContainer& filenames, bool correctTilt, const GantryTiltInformation& tiltInfo ) { if ( filenames.empty() ) { MITK_DEBUG << "Calling LoadDicomSeries with empty filename string container. Probably invalid application logic."; return nullptr; // this is not actually an error but the result is very simple } typedef itk::GDCMImageIO DcmIoType; DcmIoType::Pointer io = DcmIoType::New(); try { if ( io->CanReadFile( filenames.front().c_str() ) ) { io->SetFileName( filenames.front().c_str() ); io->ReadImageInformation(); if ( io->GetPixelType() == itk::ImageIOBase::SCALAR ) { switch ( io->GetComponentType() ) { switch3DCase(DcmIoType::UCHAR, unsigned char) switch3DCase(DcmIoType::CHAR, char) switch3DCase( DcmIoType::USHORT, unsigned short) switch3DCase(DcmIoType::SHORT, short) switch3DCase(DcmIoType::UINT, unsigned int) switch3DCase(DcmIoType::INT, int) switch3DCase( DcmIoType::ULONG, long unsigned int) switch3DCase(DcmIoType::LONG, long int) switch3DCase(DcmIoType::FLOAT, float) switch3DCase(DcmIoType::DOUBLE, double) default : MITK_ERROR << "Found unsupported DICOM scalar pixel type: (enum value) " << io->GetComponentType(); } } else if ( io->GetPixelType() == itk::ImageIOBase::RGB ) { switch ( io->GetComponentType() ) { switch3DCase(DcmIoType::UCHAR, itk::RGBPixel) switch3DCase( DcmIoType::CHAR, itk::RGBPixel) switch3DCase(DcmIoType::USHORT, itk::RGBPixel) switch3DCase(DcmIoType::SHORT, itk::RGBPixel) switch3DCase( DcmIoType::UINT, itk::RGBPixel) switch3DCase(DcmIoType::INT, itk::RGBPixel) switch3DCase(DcmIoType::ULONG, itk::RGBPixel) switch3DCase(DcmIoType::LONG, itk::RGBPixel) switch3DCase( DcmIoType::FLOAT, itk::RGBPixel) switch3DCase(DcmIoType::DOUBLE, itk::RGBPixel) default : MITK_ERROR << "Found unsupported DICOM scalar pixel type: (enum value) " << io->GetComponentType(); } } MITK_ERROR << "Unsupported DICOM pixel type"; return nullptr; } } catch ( const itk::MemoryAllocationError& e ) { MITK_ERROR << "Out of memory. Cannot load DICOM series: " << e.what(); } catch ( const std::exception& e ) { MITK_ERROR << "Error encountered when loading DICOM series:" << e.what(); } catch ( ... ) { MITK_ERROR << "Unspecified error encountered when loading DICOM series."; } return nullptr; } #define switch3DnTCase( IOType, T ) \ case IOType: \ return LoadDICOMByITK3DnT( filenamesLists, correctTilt, tiltInfo, io ); mitk::Image::Pointer mitk::ITKDICOMSeriesReaderHelper::Load3DnT( const StringContainerList& filenamesLists, bool correctTilt, const GantryTiltInformation& tiltInfo ) { if ( filenamesLists.empty() || filenamesLists.front().empty() ) { MITK_DEBUG << "Calling LoadDicomSeries with empty filename string container. Probably invalid application logic."; return nullptr; // this is not actually an error but the result is very simple } typedef itk::GDCMImageIO DcmIoType; DcmIoType::Pointer io = DcmIoType::New(); try { if ( io->CanReadFile( filenamesLists.front().front().c_str() ) ) { io->SetFileName( filenamesLists.front().front().c_str() ); io->ReadImageInformation(); if ( io->GetPixelType() == itk::ImageIOBase::SCALAR ) { switch ( io->GetComponentType() ) { switch3DnTCase(DcmIoType::UCHAR, unsigned char) switch3DnTCase(DcmIoType::CHAR, char) switch3DnTCase(DcmIoType::USHORT, unsigned short) switch3DnTCase( DcmIoType::SHORT, short) switch3DnTCase(DcmIoType::UINT, unsigned int) switch3DnTCase(DcmIoType::INT, int) switch3DnTCase(DcmIoType::ULONG, long unsigned int) switch3DnTCase(DcmIoType::LONG, long int) switch3DnTCase(DcmIoType::FLOAT, float) switch3DnTCase(DcmIoType::DOUBLE, double) default : MITK_ERROR << "Found unsupported DICOM scalar pixel type: (enum value) " << io->GetComponentType(); } } else if ( io->GetPixelType() == itk::ImageIOBase::RGB ) { switch ( io->GetComponentType() ) { switch3DnTCase(DcmIoType::UCHAR, itk::RGBPixel) switch3DnTCase(DcmIoType::CHAR, itk::RGBPixel) switch3DnTCase( DcmIoType::USHORT, itk::RGBPixel) switch3DnTCase(DcmIoType::SHORT, itk::RGBPixel) switch3DnTCase(DcmIoType::UINT, itk::RGBPixel) switch3DnTCase( DcmIoType::INT, itk::RGBPixel) switch3DnTCase(DcmIoType::ULONG, itk::RGBPixel) switch3DnTCase(DcmIoType::LONG, itk::RGBPixel) switch3DnTCase( DcmIoType::FLOAT, itk::RGBPixel) switch3DnTCase(DcmIoType::DOUBLE, itk::RGBPixel) default : MITK_ERROR << "Found unsupported DICOM scalar pixel type: (enum value) " << io->GetComponentType(); } } MITK_ERROR << "Unsupported DICOM pixel type"; return nullptr; } } catch ( const itk::MemoryAllocationError& e ) { MITK_ERROR << "Out of memory. Cannot load DICOM series: " << e.what(); } catch ( const std::exception& e ) { MITK_ERROR << "Error encountered when loading DICOM series:" << e.what(); } catch ( ... ) { MITK_ERROR << "Unspecified error encountered when loading DICOM series."; } return nullptr; } bool ConvertDICOMDateTimeString( const std::string& dateString, const std::string& timeString, OFDateTime& time ) { OFString content( timeString.c_str() ); if ( !dateString.empty() ) { content = OFString( dateString.c_str() ).append( content ); } else { // This is a workaround for DICOM data that has an AquisitionTime but no AquisitionDate. // In this case, we use the current date. That's not really nice, but is absolutely OK // as we're only interested in the time anyways... OFString currentDate; DcmDate::getCurrentDate( currentDate ); content = currentDate.append( content ); } const OFCondition result = DcmDateTime::getOFDateTimeFromString( content, time ); return result.good(); } boost::posix_time::ptime ConvertOFDateTimeToPTime( const OFDateTime& time ) { const boost::gregorian::date boostDate( time.getDate().getYear(), time.getDate().getMonth(), time.getDate().getDay() ); const boost::posix_time::time_duration boostTime = boost::posix_time::hours( time.getTime().getHour() ) + boost::posix_time::minutes( time.getTime().getMinute() ) + boost::posix_time::seconds( time.getTime().getSecond() ) + boost::posix_time::milliseconds( time.getTime().getMilliSecond() ); boost::posix_time::ptime result( boostDate, boostTime ); return result; } OFDateTime GetLowerDateTime( const OFDateTime& time1, const OFDateTime& time2 ) { OFDateTime result = time1; if ( ( time2.getDate() < time1.getDate() ) || ( ( time2.getDate() == time1.getDate() ) && ( time2.getTime() < time1.getTime() ) ) ) { result = time2; } return result; } OFDateTime GetUpperDateTime( const OFDateTime& time1, const OFDateTime& time2 ) { OFDateTime result = time1; if ( ( time2.getDate() > time1.getDate() ) || ( ( time2.getDate() == time1.getDate() ) && ( time2.getTime() > time1.getTime() ) ) ) { result = time2; } return result; } double ComputeMiliSecDuration( const OFDateTime& start, const OFDateTime& stop ) { const boost::posix_time::ptime startTime = ConvertOFDateTimeToPTime( start ); const boost::posix_time::ptime stopTime = ConvertOFDateTimeToPTime( stop ); ::boost::posix_time::time_duration duration = stopTime - startTime; return duration.total_milliseconds(); } -bool mitk::ITKDICOMSeriesReaderHelper::ExtractDateTimeBoundsOfTimeStep( - const StringContainer& filenamesOfTimeStep, DateTimeBounds& bounds ) +bool mitk::ITKDICOMSeriesReaderHelper::ExtractDateTimeBoundsAndTriggerOfTimeStep( + const StringContainer& filenamesOfTimeStep, DateTimeBounds& bounds, TimeBounds& triggerBounds) { DICOMGDCMTagScanner::Pointer filescanner = DICOMGDCMTagScanner::New(); - filescanner->SetInputFiles( filenamesOfTimeStep ); - filescanner->AddTag( AcquisitionDateTag ); - filescanner->AddTag( AcquisitionTimeTag ); - filescanner->Scan(); - - const DICOMDatasetAccessingIDICOMDatasetAccessingImageFrameListmageFrameList frameList = filescanner->GetFrameInfoList(); - - bool result = false; - bool first = true; - - for ( auto pos = frameList.cbegin(); pos != frameList.cend(); ++pos ) - { - const std::string dateStr = ( *pos )->GetTagValueAsString( AcquisitionDateTag ).value; - const std::string timeStr = ( *pos )->GetTagValueAsString( AcquisitionTimeTag ).value; - - OFDateTime time; - const bool convertResult = ConvertDICOMDateTimeString( dateStr, timeStr, time ); - - if ( convertResult ) - { - if ( first ) - { - bounds[0] = time; - bounds[1] = time; - first = false; - } - else - { - bounds[0] = GetLowerDateTime( bounds[0], time ); - bounds[1] = GetUpperDateTime( bounds[1], time ); - } - result = true; - } - } - - return result; -}; - -bool mitk::ITKDICOMSeriesReaderHelper::ExtractDateTimeBoundsOfTimeStep( - const StringContainer& filenamesOfTimeStep, TimeBounds& bounds, const OFDateTime& baselineDateTime ) -{ - DICOMGDCMTagScanner::Pointer filescanner = DICOMGDCMTagScanner::New(); - filescanner->SetInputFiles( filenamesOfTimeStep ); - filescanner->AddTag( AcquisitionDateTag ); - filescanner->AddTag( AcquisitionTimeTag ); - filescanner->AddTag( TriggerTimeTag ); + filescanner->SetInputFiles(filenamesOfTimeStep); + filescanner->AddTag(AcquisitionDateTag); + filescanner->AddTag(AcquisitionTimeTag); + filescanner->AddTag(TriggerTimeTag); filescanner->Scan(); const DICOMDatasetAccessingImageFrameList frameList = filescanner->GetFrameInfoList(); bool result = false; - bool firstAq = true; - bool firstTr = true; - - DateTimeBounds aqDTBounds; - TimeBounds triggerDTBounds( 0.0 ); + bool firstAq = true; + bool firstTr = true; + triggerBounds = TimeBounds(0.0); - for ( DICOMDatasetAccessingImageFrameList::const_iterator pos = frameList.cbegin(); pos != frameList.cend(); ++pos ) + for (DICOMDatasetAccessingImageFrameList::const_iterator pos = frameList.cbegin(); pos != frameList.cend(); ++pos) { - const std::string aqDateStr = ( *pos )->GetTagValueAsString( AcquisitionDateTag ).value; - const std::string aqTimeStr = ( *pos )->GetTagValueAsString( AcquisitionTimeTag ).value; - const std::string triggerTimeStr = ( *pos )->GetTagValueAsString( TriggerTimeTag ).value; + const std::string aqDateStr = (*pos)->GetTagValueAsString(AcquisitionDateTag).value; + const std::string aqTimeStr = (*pos)->GetTagValueAsString(AcquisitionTimeTag).value; + const std::string triggerTimeStr = (*pos)->GetTagValueAsString(TriggerTimeTag).value; OFDateTime aqDateTime; - const bool convertAqResult = ConvertDICOMDateTimeString( aqDateStr, aqTimeStr, aqDateTime ); + const bool convertAqResult = ConvertDICOMDateTimeString(aqDateStr, aqTimeStr, aqDateTime); OFBool convertTriggerResult; - mitk::ScalarType triggerTime = OFStandard::atof( triggerTimeStr.c_str(), &convertTriggerResult ); + mitk::ScalarType triggerTime = OFStandard::atof(triggerTimeStr.c_str(), &convertTriggerResult); - if ( convertAqResult ) + if (convertAqResult) { - if ( firstAq ) + if (firstAq) { - aqDTBounds[0] = aqDateTime; - aqDTBounds[1] = aqDateTime; - firstAq = false; + bounds[0] = aqDateTime; + bounds[1] = aqDateTime; + firstAq = false; } else { - aqDTBounds[0] = GetLowerDateTime( aqDTBounds[0], aqDateTime ); - aqDTBounds[1] = GetUpperDateTime( aqDTBounds[1], aqDateTime ); + bounds[0] = GetLowerDateTime(bounds[0], aqDateTime); + bounds[1] = GetUpperDateTime(bounds[1], aqDateTime); } result = true; } - if ( convertTriggerResult ) + if (convertTriggerResult) { - if ( firstTr ) + if (firstTr) { - triggerDTBounds[0] = triggerTime; - triggerDTBounds[1] = triggerTime; - firstTr = false; + triggerBounds[0] = triggerTime; + triggerBounds[1] = triggerTime; + firstTr = false; } else { - triggerDTBounds[0] = std::min( triggerDTBounds[0], triggerTime ); - triggerDTBounds[1] = std::max( triggerDTBounds[1], triggerTime ); + triggerBounds[0] = std::min(triggerBounds[0], triggerTime); + triggerBounds[1] = std::max(triggerBounds[1], triggerTime); } result = true; } } + return result; +}; + +bool mitk::ITKDICOMSeriesReaderHelper::ExtractTimeBoundsOfTimeStep( + const StringContainer& filenamesOfTimeStep, TimeBounds& bounds, const OFDateTime& baselineDateTime ) +{ + DateTimeBounds aqDTBounds; + TimeBounds triggerBounds; + + bool result = ExtractDateTimeBoundsAndTriggerOfTimeStep(filenamesOfTimeStep, aqDTBounds, triggerBounds); + mitk::ScalarType lowerBound = ComputeMiliSecDuration( baselineDateTime, aqDTBounds[0] ); mitk::ScalarType upperBound = ComputeMiliSecDuration( baselineDateTime, aqDTBounds[1] ); if ( lowerBound < mitk::eps || upperBound < mitk::eps ) { - lowerBound = triggerDTBounds[0]; - upperBound = triggerDTBounds[1]; + lowerBound = triggerBounds[0]; + upperBound = triggerBounds[1]; } bounds[0] = lowerBound; bounds[1] = upperBound; return result; }; mitk::ITKDICOMSeriesReaderHelper::TimeBoundsList mitk::ITKDICOMSeriesReaderHelper::ExtractTimeBoundsOfTimeSteps( const StringContainerList& filenamesOfTimeSteps ) { TimeBoundsList result; OFDateTime baseLine; bool baseLineSet = false; // extract the timebounds DateTimeBounds baselineDateTimeBounds; + TimeBounds triggerBounds; StringContainerList::const_iterator pos = filenamesOfTimeSteps.cbegin(); - ExtractDateTimeBoundsOfTimeStep( *pos, baselineDateTimeBounds ); + ExtractDateTimeBoundsAndTriggerOfTimeStep(*pos, baselineDateTimeBounds, triggerBounds); baseLine = baselineDateTimeBounds[0]; baseLineSet = true; // timebounds for baseline is 0 TimeBounds bounds( 0.0 ); result.push_back( bounds ); // iterate over the remaining timesteps for ( ++pos; pos != filenamesOfTimeSteps.cend(); ++pos ) { TimeBounds bounds( 0.0 ); TimeBounds dateTimeBounds; // extract the timebounds relative to the baseline - if ( ExtractDateTimeBoundsOfTimeStep( *pos, dateTimeBounds, baseLine ) ) + if ( ExtractTimeBoundsOfTimeStep( *pos, dateTimeBounds, baseLine ) ) { bounds[0] = dateTimeBounds[0]; bounds[1] = dateTimeBounds[1]; } result.push_back( bounds ); } return result; }; mitk::TimeGeometry::Pointer mitk::ITKDICOMSeriesReaderHelper::GenerateTimeGeometry( const BaseGeometry* templateGeometry, const TimeBoundsList& boundsList ) { TimeGeometry::Pointer timeGeometry; double check = 0.0; const auto boundListSize = boundsList.size(); for ( auto pos = 0; pos < boundListSize; ++pos ) { check += boundsList[pos][0]; check += boundsList[pos][1]; } if ( check < mitk::eps ) { // if all bounds are zero we assume that the bounds could not be correctly determined // and as a fallback generate a time geometry in the old mitk style ProportionalTimeGeometry::Pointer newTimeGeometry = ProportionalTimeGeometry::New(); newTimeGeometry->Initialize( templateGeometry, boundListSize ); timeGeometry = newTimeGeometry.GetPointer(); } else { ArbitraryTimeGeometry::Pointer newTimeGeometry = ArbitraryTimeGeometry::New(); newTimeGeometry->ClearAllGeometries(); newTimeGeometry->ReserveSpaceForGeometries( boundListSize ); for ( auto pos = 0; pos < boundListSize; ++pos ) { TimeBounds bounds = boundsList[pos]; if ( pos + 1 < boundListSize ) { //Currently we do not explicitly support "gaps" in the time coverage //thus we set the max time bound of a time step to the min time bound //of its successor. bounds[1] = boundsList[pos + 1][0]; } newTimeGeometry->AppendNewTimeStepClone(templateGeometry, bounds[0], bounds[1]); } timeGeometry = newTimeGeometry.GetPointer(); } return timeGeometry; };