diff --git a/Modules/DICOM/include/mitkITKDICOMSeriesReaderHelper.h b/Modules/DICOM/include/mitkITKDICOMSeriesReaderHelper.h index 4bcae594c7..68816b35ac 100644 --- a/Modules/DICOM/include/mitkITKDICOMSeriesReaderHelper.h +++ b/Modules/DICOM/include/mitkITKDICOMSeriesReaderHelper.h @@ -1,105 +1,109 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkITKDICOMSeriesReaderHelper_h #define mitkITKDICOMSeriesReaderHelper_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 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 ExtractDateTimeBoundsAndTriggerOfTimeStep( const StringContainer& filenamesOfTimeStep, DateTimeBounds& bounds, TimeBounds& triggerBounds); /* Determine the time bounds in ms respective to the baselineDateTime for the passed files. Additionally it regards the trigger time tag if set and acquisition date time - carries not enough information.*/ + carries not enough information. + @param [in, out] usedTriggerBounds Variable used to communicate if preceding time bounds extractions + used the trigger time (true). If true is passed to the method, it will use always trigger time + to extract the next bounds. Default is false. After the method call is returning the variable indicates + the last call used trigger time or not to extract the time bounds.*/ static bool ExtractTimeBoundsOfTimeStep(const StringContainer& filenamesOfTimeStep, TimeBounds& bounds, - const OFDateTime& baselineDateTime ); - + const OFDateTime& baselineDateTime, + bool& usedTriggerBounds); /** Returns the list of time bounds of all passed time step containers. (sa ExtractTimeBoundsOfTimeStep and ExtractDateTimeBoundsOfTimeStep). Time steps where no time bounds could be extracted are indicated by "null" time bounds (both values "0"). The order of the returned 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 + /** Helper function that generates a time geometry using the template and the passed bounds list (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/DICOM/src/mitkITKDICOMSeriesReaderHelper.cpp b/Modules/DICOM/src/mitkITKDICOMSeriesReaderHelper.cpp index 5ac21da5d4..ff93c00ab8 100644 --- a/Modules/DICOM/src/mitkITKDICOMSeriesReaderHelper.cpp +++ b/Modules/DICOM/src/mitkITKDICOMSeriesReaderHelper.cpp @@ -1,455 +1,464 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ //#define MBILOG_ENABLE_DEBUG #include #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 "dcmtk/dcmdata/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::IOPixelEnum::SCALAR ) { switch ( io->GetComponentType() ) { switch3DCase(itk::IOComponentEnum::UCHAR, unsigned char) switch3DCase(itk::IOComponentEnum::CHAR, char) switch3DCase( itk::IOComponentEnum::USHORT, unsigned short) switch3DCase(itk::IOComponentEnum::SHORT, short) switch3DCase(itk::IOComponentEnum::UINT, unsigned int) switch3DCase(itk::IOComponentEnum::INT, int) switch3DCase( itk::IOComponentEnum::ULONG, long unsigned int) switch3DCase(itk::IOComponentEnum::LONG, long int) switch3DCase(itk::IOComponentEnum::FLOAT, float) switch3DCase(itk::IOComponentEnum::DOUBLE, double) default : MITK_ERROR << "Found unsupported DICOM scalar pixel type: (enum value) " << io->GetComponentType(); } } else if ( io->GetPixelType() == itk::IOPixelEnum::RGB ) { switch ( io->GetComponentType() ) { switch3DCase(itk::IOComponentEnum::UCHAR, itk::RGBPixel) switch3DCase( itk::IOComponentEnum::CHAR, itk::RGBPixel) switch3DCase(itk::IOComponentEnum::USHORT, itk::RGBPixel) switch3DCase(itk::IOComponentEnum::SHORT, itk::RGBPixel) switch3DCase( itk::IOComponentEnum::UINT, itk::RGBPixel) switch3DCase(itk::IOComponentEnum::INT, itk::RGBPixel) switch3DCase(itk::IOComponentEnum::ULONG, itk::RGBPixel) switch3DCase(itk::IOComponentEnum::LONG, itk::RGBPixel) switch3DCase( itk::IOComponentEnum::FLOAT, itk::RGBPixel) switch3DCase(itk::IOComponentEnum::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::IOPixelEnum::SCALAR ) { switch ( io->GetComponentType() ) { switch3DnTCase(itk::IOComponentEnum::UCHAR, unsigned char) switch3DnTCase(itk::IOComponentEnum::CHAR, char) switch3DnTCase(itk::IOComponentEnum::USHORT, unsigned short) switch3DnTCase( itk::IOComponentEnum::SHORT, short) switch3DnTCase(itk::IOComponentEnum::UINT, unsigned int) switch3DnTCase(itk::IOComponentEnum::INT, int) switch3DnTCase(itk::IOComponentEnum::ULONG, long unsigned int) switch3DnTCase(itk::IOComponentEnum::LONG, long int) switch3DnTCase(itk::IOComponentEnum::FLOAT, float) switch3DnTCase(itk::IOComponentEnum::DOUBLE, double) default : MITK_ERROR << "Found unsupported DICOM scalar pixel type: (enum value) " << io->GetComponentType(); } } else if ( io->GetPixelType() == itk::IOPixelEnum::RGB ) { switch ( io->GetComponentType() ) { switch3DnTCase(itk::IOComponentEnum::UCHAR, itk::RGBPixel) switch3DnTCase(itk::IOComponentEnum::CHAR, itk::RGBPixel) switch3DnTCase( itk::IOComponentEnum::USHORT, itk::RGBPixel) switch3DnTCase(itk::IOComponentEnum::SHORT, itk::RGBPixel) switch3DnTCase(itk::IOComponentEnum::UINT, itk::RGBPixel) switch3DnTCase( itk::IOComponentEnum::INT, itk::RGBPixel) switch3DnTCase(itk::IOComponentEnum::ULONG, itk::RGBPixel) switch3DnTCase(itk::IOComponentEnum::LONG, itk::RGBPixel) switch3DnTCase( itk::IOComponentEnum::FLOAT, itk::RGBPixel) switch3DnTCase(itk::IOComponentEnum::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( static_cast(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::ExtractDateTimeBoundsAndTriggerOfTimeStep( const StringContainer& filenamesOfTimeStep, DateTimeBounds& bounds, TimeBounds& triggerBounds) { DICOMGDCMTagScanner::Pointer filescanner = DICOMGDCMTagScanner::New(); 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; triggerBounds = TimeBounds(0.0); for (auto 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; OFDateTime aqDateTime; const bool convertAqResult = ConvertDICOMDateTimeString(aqDateStr, aqTimeStr, aqDateTime); OFBool convertTriggerResult; mitk::ScalarType triggerTime = OFStandard::atof(triggerTimeStr.c_str(), &convertTriggerResult); if (convertAqResult) { if (firstAq) { bounds[0] = aqDateTime; bounds[1] = aqDateTime; firstAq = false; } else { bounds[0] = GetLowerDateTime(bounds[0], aqDateTime); bounds[1] = GetUpperDateTime(bounds[1], aqDateTime); } result = true; } if (convertTriggerResult) { if (firstTr) { triggerBounds[0] = triggerTime; triggerBounds[1] = triggerTime; firstTr = false; } else { 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 ) + const StringContainer& filenamesOfTimeStep, TimeBounds& bounds, const OFDateTime& baselineDateTime, bool& usedTriggerBounds ) { 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 ) + if ( lowerBound < mitk::eps || upperBound < mitk::eps || usedTriggerBounds) { lowerBound = triggerBounds[0]; upperBound = triggerBounds[1]; + usedTriggerBounds = true; + } + else + { + usedTriggerBounds = false; } bounds[0] = lowerBound; bounds[1] = upperBound; return result; }; mitk::ITKDICOMSeriesReaderHelper::TimeBoundsList mitk::ITKDICOMSeriesReaderHelper::ExtractTimeBoundsOfTimeSteps( const StringContainerList& filenamesOfTimeSteps ) { TimeBoundsList result; OFDateTime baseLine; // extract the timebounds DateTimeBounds baselineDateTimeBounds; TimeBounds triggerBounds; auto pos = filenamesOfTimeSteps.cbegin(); ExtractDateTimeBoundsAndTriggerOfTimeStep(*pos, baselineDateTimeBounds, triggerBounds); baseLine = baselineDateTimeBounds[0]; // timebounds for baseline is 0 TimeBounds bounds( 0.0 ); result.push_back( bounds ); + //start with not using trigger time. Will be changed by ExtractTimeBoundsOfTimeStep if needed. + bool usedTriggerTime = false; // 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 ( ExtractTimeBoundsOfTimeStep( *pos, dateTimeBounds, baseLine ) ) + if ( ExtractTimeBoundsOfTimeStep( *pos, dateTimeBounds, baseLine, usedTriggerTime) ) { - bounds[0] = dateTimeBounds[0]; bounds[1] = dateTimeBounds[1]; } result.push_back( bounds ); } + if (usedTriggerTime) + MITK_DEBUG << "Used trigger time to extract time bounds of passed files"; + 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 ( std::size_t 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 ( std::size_t 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; };