diff --git a/Modules/Pharmacokinetics/include/mitkConcentrationCurveGenerator.h b/Modules/Pharmacokinetics/include/mitkConcentrationCurveGenerator.h index 6536e5f630..9622837ac6 100644 --- a/Modules/Pharmacokinetics/include/mitkConcentrationCurveGenerator.h +++ b/Modules/Pharmacokinetics/include/mitkConcentrationCurveGenerator.h @@ -1,161 +1,161 @@ /*=================================================================== 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 CONCENTRATIONCURVEGENERATOR_H #define CONCENTRATIONCURVEGENERATOR_H #include #include #include "mitkConvertToConcentrationAbsoluteFunctor.h" #include "mitkConvertToConcentrationRelativeFunctor.h" #include "MitkPharmacokineticsExports.h" namespace mitk { /** \class ConcentrationCurveGenerator * \brief Converts a given 4D mitk::Image with MR signal values into a 4D mitk::Image with corresponding contrast agent concentration values * * From a given 4D image, the Generator takes the 3D image of the first time point as baseline image. It then loops over all time steps, casts * the current 3D image to itk and passes it to the ConvertToconcentrationFunctor. The returned 3D image has now values of concentration type and is stored at its timepoint * in the return image. */ class MITKPHARMACOKINETICS_EXPORT ConcentrationCurveGenerator : public itk::Object { public: mitkClassMacroItkParent(ConcentrationCurveGenerator, itk::Object); itkNewMacro(Self); //typedef itk::Image ImageType; typedef itk::Image ConvertedImageType; /** Getter and Setter for 4D mitk::Image*/ itkSetConstObjectMacro(DynamicImage,Image); itkGetConstObjectMacro(DynamicImage,Image); /** Parameters Relevant for conversion Calculation; Have to be Set externally (Sequence Dependend)*/ itkSetMacro(RelaxationTime, double); itkGetConstReferenceMacro(RelaxationTime, double); itkSetMacro(Relaxivity, double); itkGetConstReferenceMacro(Relaxivity, double); itkSetMacro(RecoveryTime, double); itkGetConstReferenceMacro(RecoveryTime, double); itkSetMacro(FlipAngle, double); itkGetConstReferenceMacro(FlipAngle, double); itkSetMacro(Factor, double); itkGetConstReferenceMacro(Factor, double); /** Getter and Setter for T10 Map image*/ itkSetConstObjectMacro(T10Image,Image); itkGetConstObjectMacro(T10Image,Image); itkSetMacro(T2Factor, double); itkGetConstReferenceMacro(T2Factor, double); itkSetMacro(T2EchoTime, double); itkGetConstReferenceMacro(T2EchoTime, double); /** @brief Calls Convert and returns the 4D mitk::image in Concentration units*/ - itkSetMacro(BaselineStartTimeStep, int); - itkGetConstReferenceMacro(BaselineStartTimeStep, int); + itkSetMacro(BaselineStartTimeStep, unsigned int); + itkGetConstReferenceMacro(BaselineStartTimeStep, unsigned int); - itkSetMacro(BaselineEndTimeStep, int); - itkGetConstReferenceMacro(BaselineEndTimeStep, int); + itkSetMacro(BaselineEndTimeStep, unsigned int); + itkGetConstReferenceMacro(BaselineEndTimeStep, unsigned int); itkSetMacro(isTurboFlashSequence,bool); itkGetConstReferenceMacro(isTurboFlashSequence,bool); itkSetMacro(AbsoluteSignalEnhancement,bool); itkGetConstReferenceMacro(AbsoluteSignalEnhancement,bool); itkSetMacro(RelativeSignalEnhancement,bool); itkGetConstReferenceMacro(RelativeSignalEnhancement,bool); itkSetMacro(UsingT1Map,bool); itkGetConstReferenceMacro(UsingT1Map,bool); itkSetMacro(isT2weightedImage,bool); itkGetConstReferenceMacro(isT2weightedImage,bool); Image::Pointer GetConvertedImage(); protected: ConcentrationCurveGenerator(); ~ConcentrationCurveGenerator() override; template mitk::Image::Pointer convertToConcentration(const itk::Image *itkInputImage, const itk::Image *itkBaselineImage); /** Calls ConvertToconcentrationFunctor for passed 3D itk::image*/ mitk::Image::Pointer ConvertSignalToConcentrationCurve(const mitk::Image* inputImage, const mitk::Image* baselineImage); /** @brief Takes the 3D image of the first timepoint to set as baseline image*/ void PrepareBaselineImage(); template void CalculateAverageBaselineImage(const itk::Image *itkBaselineImage); /** @brief loops over all timepoints, casts the current timepoint 3D mitk::image to itk and passes it to ConvertSignalToConcentrationCurve */ virtual void Convert(); private: Image::ConstPointer m_DynamicImage; Image::ConstPointer m_BaselineImage; Image::ConstPointer m_T10Image; Image::Pointer m_ConvertSignalToConcentrationCurve_OutputImage; Image::Pointer m_ConvertedImage; bool m_isT2weightedImage; bool m_isTurboFlashSequence; bool m_AbsoluteSignalEnhancement; bool m_RelativeSignalEnhancement; bool m_UsingT1Map; double m_Factor; //=Repetition Time TR double m_RecoveryTime; //= pre-CA T1 time double m_RelaxationTime; //= contrast agent relaxivity double m_Relaxivity; double m_FlipAngle; double m_T2Factor; double m_T2EchoTime; // The baseline image is averaged from the signal within time step range [m_BaselineStartTimeStep, m_BaselineEndTimeStep]. // m_BaselineStartTimeStep is the first time frame, that is included into the baseline averaging (starting with 0). - int m_BaselineStartTimeStep; + unsigned int m_BaselineStartTimeStep; // m_BaselinStopTimeStep is the last time frame, that is included into the baseline averaging. - int m_BaselineEndTimeStep; + unsigned int m_BaselineEndTimeStep; }; } #endif // CONCENTRATIONCURVEGENERATOR_H diff --git a/Modules/Pharmacokinetics/src/Common/mitkConcentrationCurveGenerator.cpp b/Modules/Pharmacokinetics/src/Common/mitkConcentrationCurveGenerator.cpp index c2c685f9f2..4265cd810d 100644 --- a/Modules/Pharmacokinetics/src/Common/mitkConcentrationCurveGenerator.cpp +++ b/Modules/Pharmacokinetics/src/Common/mitkConcentrationCurveGenerator.cpp @@ -1,327 +1,327 @@ #include "mitkConcentrationCurveGenerator.h" #include "mitkConvertToConcentrationTurboFlashFunctor.h" #include "mitkConvertT2ConcentrationFunctor.h" #include "mitkConvertToConcentrationViaT1Functor.h" #include "mitkImageTimeSelector.h" #include "mitkImageCast.h" #include "mitkITKImageImport.h" #include "mitkModelBase.h" #include "mitkExtractTimeGrid.h" #include "mitkArbitraryTimeGeometry.h" #include "itkNaryAddImageFilter.h" #include "mitkImageAccessByItk.h" #include "itkImageIOBase.h" #include "itkBinaryFunctorImageFilter.h" #include "itkTernaryFunctorImageFilter.h" #include #include "itkMeanProjectionImageFilter.h" mitk::ConcentrationCurveGenerator::ConcentrationCurveGenerator() : m_isT2weightedImage(false), m_isTurboFlashSequence(false), m_AbsoluteSignalEnhancement(false), m_RelativeSignalEnhancement(0.0), m_UsingT1Map(false), m_Factor(0.0), m_RecoveryTime(0.0), m_RelaxationTime(0.0), m_Relaxivity(0.0), m_FlipAngle(0.0), m_T2Factor(0.0), m_T2EchoTime(0.0) { } mitk::ConcentrationCurveGenerator::~ConcentrationCurveGenerator() { } mitk::Image::Pointer mitk::ConcentrationCurveGenerator::GetConvertedImage() { if(this->m_DynamicImage.IsNull()) { itkExceptionMacro( << "Dynamic Image not set!"); } else { Convert(); } return m_ConvertedImage; } void mitk::ConcentrationCurveGenerator::Convert() { mitk::Image::Pointer tempImage = mitk::Image::New(); mitk::PixelType pixeltype = mitk::MakeScalarPixelType(); tempImage->Initialize(pixeltype,*this->m_DynamicImage->GetTimeGeometry()); mitk::TimeGeometry::Pointer timeGeometry = (this->m_DynamicImage->GetTimeGeometry())->Clone(); tempImage->SetTimeGeometry(timeGeometry); PrepareBaselineImage(); mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); imageTimeSelector->SetInput(this->m_DynamicImage); for(unsigned int i = 0; i< this->m_DynamicImage->GetTimeSteps(); ++i) { imageTimeSelector->SetTimeNr(i); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image::Pointer mitkInputImage = imageTimeSelector->GetOutput(); mitk::Image::Pointer outputImage = mitk::Image::New(); outputImage = ConvertSignalToConcentrationCurve(mitkInputImage,this->m_BaselineImage); mitk::ImageReadAccessor accessor(outputImage); tempImage->SetVolume(accessor.GetData(), i); } this->m_ConvertedImage = tempImage; } void mitk::ConcentrationCurveGenerator::PrepareBaselineImage() { mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); imageTimeSelector->SetInput(this->m_DynamicImage); imageTimeSelector->SetTimeNr(0); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image::Pointer baselineImage; baselineImage = imageTimeSelector->GetOutput(); if (m_BaselineStartTimeStep == m_BaselineEndTimeStep) { this->m_BaselineImage = imageTimeSelector->GetOutput(); } else { try { AccessFixedDimensionByItk(this->m_DynamicImage, mitk::ConcentrationCurveGenerator::CalculateAverageBaselineImage, 4); } catch (itk::ExceptionObject & err) { std::cerr << "ExceptionObject in ConcentrationCurveGenerator::CalculateAverageBaselineImage caught!" << std::endl; std::cerr << err << std::endl; } } } template void mitk::ConcentrationCurveGenerator::CalculateAverageBaselineImage(const itk::Image *itkBaselineImage) { if (itkBaselineImage == NULL) { mitkThrow() << "Error in ConcentrationCurveGenerator::CalculateAverageBaselineImage. Input image is NULL."; } if (m_BaselineStartTimeStep > m_BaselineEndTimeStep) { mitkThrow() << "Error in ConcentrationCurveGenerator::CalculateAverageBaselineImage. End time point is before start time point."; } typedef itk::Image TPixel4DImageType; typedef itk::Image TPixel3DImageType; typedef itk::Image Double3DImageType; typedef itk::Image Double4DImageType; typedef itk::ExtractImageFilter ExtractImageFilterType; typedef itk::ExtractImageFilter Extract3DImageFilterType; typedef itk::MeanProjectionImageFilter MeanProjectionImageFilterType; - MeanProjectionImageFilterType::Pointer MeanProjectionImageFilter = MeanProjectionImageFilterType::New(); - Extract3DImageFilterType::Pointer Extract3DImageFilter = Extract3DImageFilterType::New(); - TPixel4DImageType::RegionType region_input = itkBaselineImage->GetLargestPossibleRegion(); + typename MeanProjectionImageFilterType::Pointer MeanProjectionImageFilter = MeanProjectionImageFilterType::New(); + typename Extract3DImageFilterType::Pointer Extract3DImageFilter = Extract3DImageFilterType::New(); + typename TPixel4DImageType::RegionType region_input = itkBaselineImage->GetLargestPossibleRegion(); if (m_BaselineEndTimeStep > region_input.GetSize()[3]) { mitkThrow() << "Error in ConcentrationCurveGenerator::CalculateAverageBaselineImage. End time point is larger than total number of time points."; } ExtractImageFilterType::Pointer ExtractFilter = ExtractImageFilterType::New(); - TPixel4DImageType::Pointer baselineTimeFrameImage = TPixel4DImageType::New(); - TPixel4DImageType::RegionType extractionRegion; - TPixel4DImageType::SizeType size_input_aux = region_input.GetSize(); + typename TPixel4DImageType::Pointer baselineTimeFrameImage = TPixel4DImageType::New(); + typename TPixel4DImageType::RegionType extractionRegion; + typename TPixel4DImageType::SizeType size_input_aux = region_input.GetSize(); size_input_aux[3] = m_BaselineEndTimeStep - m_BaselineStartTimeStep+1; - TPixel4DImageType::IndexType start_input_aux = region_input.GetIndex(); + typename TPixel4DImageType::IndexType start_input_aux = region_input.GetIndex(); start_input_aux[3] = m_BaselineStartTimeStep; extractionRegion.SetSize(size_input_aux); extractionRegion.SetIndex(start_input_aux); ExtractFilter->SetExtractionRegion(extractionRegion); ExtractFilter->SetInput(itkBaselineImage); ExtractFilter->SetDirectionCollapseToSubmatrix(); try { ExtractFilter->Update(); } catch (itk::ExceptionObject & err) { std::cerr << "ExceptionObject caught!" << std::endl; std::cerr << err << std::endl; } baselineTimeFrameImage = ExtractFilter->GetOutput(); MeanProjectionImageFilter->SetProjectionDimension(3); MeanProjectionImageFilter->SetInput(baselineTimeFrameImage); try { MeanProjectionImageFilter->Update(); } catch (itk::ExceptionObject & err) { std::cerr << "ExceptionObject caught!" << std::endl; std::cerr << err << std::endl; } Extract3DImageFilter->SetInput(MeanProjectionImageFilter->GetOutput()); size_input_aux[3] = 0; start_input_aux[3] = 0; extractionRegion.SetSize(size_input_aux); extractionRegion.SetIndex(start_input_aux); Extract3DImageFilter->SetExtractionRegion(extractionRegion); Extract3DImageFilter->SetDirectionCollapseToSubmatrix(); try { Extract3DImageFilter->Update(); } catch (itk::ExceptionObject & err) { std::cerr << "ExceptionObject caught!" << std::endl; std::cerr << err << std::endl; } Image::Pointer mitkBaselineImage = Image::New(); CastToMitkImage(Extract3DImageFilter->GetOutput(), mitkBaselineImage); this->m_BaselineImage = mitkBaselineImage; } mitk::Image::Pointer mitk::ConcentrationCurveGenerator::ConvertSignalToConcentrationCurve(const mitk::Image* inputImage,const mitk::Image* baselineImage) { mitk::PixelType m_PixelType = inputImage->GetPixelType(); AccessTwoImagesFixedDimensionByItk(inputImage, baselineImage, mitk::ConcentrationCurveGenerator::convertToConcentration, 3); return m_ConvertSignalToConcentrationCurve_OutputImage; } template mitk::Image::Pointer mitk::ConcentrationCurveGenerator::convertToConcentration(const itk::Image *itkInputImage, const itk::Image *itkBaselineImage) { typedef itk::Image InputImageType; typedef itk::Image BaselineImageType; typedef itk::Image DoubleImageType; typename DoubleImageType::Pointer itkDoubleInputImage = DoubleImageType::New(); typename DoubleImageType::Pointer itkDoubleBaselineImage = DoubleImageType::New(); typedef itk::CastImageFilter CastInputImageToDoubleImageFilterType; - CastInputImageToDoubleImageFilterType::Pointer CastInputImageToDoubleImageFilter = CastInputImageToDoubleImageFilterType::New(); + typename CastInputImageToDoubleImageFilterType::Pointer CastInputImageToDoubleImageFilter = CastInputImageToDoubleImageFilterType::New(); CastInputImageToDoubleImageFilter->SetInput(itkInputImage); CastInputImageToDoubleImageFilter->Update(); itkDoubleInputImage = CastInputImageToDoubleImageFilter->GetOutput(); typedef itk::CastImageFilter CastBaselineImageToDoubleImageFilterType; - CastBaselineImageToDoubleImageFilterType::Pointer CastBaselineImageToDoubleImageFilter = CastBaselineImageToDoubleImageFilterType::New(); + typename CastBaselineImageToDoubleImageFilterType::Pointer CastBaselineImageToDoubleImageFilter = CastBaselineImageToDoubleImageFilterType::New(); CastBaselineImageToDoubleImageFilter->SetInput(itkBaselineImage); CastBaselineImageToDoubleImageFilter->Update(); itkDoubleBaselineImage = CastBaselineImageToDoubleImageFilter->GetOutput(); if (this->m_isT2weightedImage) { typedef mitk::ConvertT2ConcentrationFunctor ConversionFunctorT2Type; typedef itk::BinaryFunctorImageFilter FilterT2Type; ConversionFunctorT2Type ConversionT2Functor; ConversionT2Functor.initialize(this->m_T2Factor, this->m_T2EchoTime); typename FilterT2Type::Pointer ConversionT2Filter = FilterT2Type::New(); ConversionT2Filter->SetFunctor(ConversionT2Functor); ConversionT2Filter->SetInput1(itkDoubleInputImage); ConversionT2Filter->SetInput2(itkDoubleBaselineImage); ConversionT2Filter->Update(); m_ConvertSignalToConcentrationCurve_OutputImage = mitk::ImportItkImage(ConversionT2Filter->GetOutput())->Clone(); } else { if(this->m_isTurboFlashSequence) { typedef mitk::ConvertToConcentrationTurboFlashFunctor ConversionFunctorTurboFlashType; typedef itk::BinaryFunctorImageFilter FilterTurboFlashType; ConversionFunctorTurboFlashType ConversionTurboFlashFunctor; ConversionTurboFlashFunctor.initialize(this->m_RelaxationTime, this->m_Relaxivity, this->m_RecoveryTime); typename FilterTurboFlashType::Pointer ConversionTurboFlashFilter = FilterTurboFlashType::New(); ConversionTurboFlashFilter->SetFunctor(ConversionTurboFlashFunctor); ConversionTurboFlashFilter->SetInput1(itkDoubleInputImage); ConversionTurboFlashFilter->SetInput2(itkDoubleBaselineImage); ConversionTurboFlashFilter->Update(); m_ConvertSignalToConcentrationCurve_OutputImage = mitk::ImportItkImage(ConversionTurboFlashFilter->GetOutput())->Clone(); } else if(this->m_UsingT1Map) { typename DoubleImageType::Pointer itkT10Image = DoubleImageType::New(); mitk::CastToItkImage(m_T10Image, itkT10Image); typedef mitk::ConvertToConcentrationViaT1CalcFunctor ConvertToConcentrationViaT1CalcFunctorType; typedef itk::TernaryFunctorImageFilter FilterT1MapType; ConvertToConcentrationViaT1CalcFunctorType ConversionT1MapFunctor; ConversionT1MapFunctor.initialize(this->m_Relaxivity, this->m_RecoveryTime, this->m_FlipAngle); typename FilterT1MapType::Pointer ConversionT1MapFilter = FilterT1MapType::New(); ConversionT1MapFilter->SetFunctor(ConversionT1MapFunctor); ConversionT1MapFilter->SetInput1(itkDoubleInputImage); ConversionT1MapFilter->SetInput2(itkDoubleBaselineImage); ConversionT1MapFilter->SetInput3(itkT10Image); ConversionT1MapFilter->Update(); m_ConvertSignalToConcentrationCurve_OutputImage = mitk::ImportItkImage(ConversionT1MapFilter->GetOutput())->Clone(); } else if(this->m_AbsoluteSignalEnhancement) { typedef mitk::ConvertToConcentrationAbsoluteFunctor ConversionFunctorAbsoluteType; typedef itk::BinaryFunctorImageFilter FilterAbsoluteType; ConversionFunctorAbsoluteType ConversionAbsoluteFunctor; ConversionAbsoluteFunctor.initialize(this->m_Factor); typename FilterAbsoluteType::Pointer ConversionAbsoluteFilter = FilterAbsoluteType::New(); ConversionAbsoluteFilter->SetFunctor(ConversionAbsoluteFunctor); ConversionAbsoluteFilter->SetInput1(itkDoubleInputImage); ConversionAbsoluteFilter->SetInput2(itkDoubleBaselineImage); ConversionAbsoluteFilter->Update(); m_ConvertSignalToConcentrationCurve_OutputImage = mitk::ImportItkImage(ConversionAbsoluteFilter->GetOutput())->Clone(); } else if(this->m_RelativeSignalEnhancement) { typedef mitk::ConvertToConcentrationRelativeFunctor ConversionFunctorRelativeType; typedef itk::BinaryFunctorImageFilter FilterRelativeType; ConversionFunctorRelativeType ConversionRelativeFunctor; ConversionRelativeFunctor.initialize(this->m_Factor); typename FilterRelativeType::Pointer ConversionRelativeFilter = FilterRelativeType::New(); ConversionRelativeFilter->SetFunctor(ConversionRelativeFunctor); ConversionRelativeFilter->SetInput1(itkDoubleInputImage); ConversionRelativeFilter->SetInput2(itkDoubleBaselineImage); ConversionRelativeFilter->Update(); m_ConvertSignalToConcentrationCurve_OutputImage = mitk::ImportItkImage(ConversionRelativeFilter->GetOutput())->Clone(); } } return m_ConvertSignalToConcentrationCurve_OutputImage; }