diff --git a/Modules/MatchPointRegistration/CMakeLists.txt b/Modules/MatchPointRegistration/CMakeLists.txt index 97bb55792b..acbcd26ff3 100644 --- a/Modules/MatchPointRegistration/CMakeLists.txt +++ b/Modules/MatchPointRegistration/CMakeLists.txt @@ -1,25 +1,28 @@ MITK_CREATE_MODULE( - INCLUDE_DIRS PUBLIC Rendering Helper algorithms + INCLUDE_DIRS + PUBLIC algorithms + PRIVATE src/Helper src/Rendering DEPENDS MitkCore MitkSceneSerializationBase PACKAGE_DEPENDS PUBLIC MatchPoint ) if(TARGET ${MODULE_TARGET}) set(ALG_PROFILE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/algorithms) include(${MatchPoint_SOURCE_DIR}/CMake/mapFunctionCreateAlgorithmProfile.cmake) file(GLOB ALG_PROFILE_FILES LIST_DIRECTORIES false RELATIVE ${ALG_PROFILE_DIR} "${ALG_PROFILE_DIR}/*.profile") foreach(profile_file ${ALG_PROFILE_FILES}) get_filename_component(profile_name ${profile_file} NAME_WE) MESSAGE(STATUS "... generate MDRA profile ${profile_name} (from ${profile_file})...") CREATE_ALGORITHM_PROFILE(${profile_name} ${ALG_PROFILE_DIR}/${profile_file}) endforeach(profile_file) ADD_SUBDIRECTORY(autoload/IO) ADD_SUBDIRECTORY(deployment) if(BUILD_TESTING) ADD_SUBDIRECTORY(Testing) endif(BUILD_TESTING) + ADD_SUBDIRECTORY(cmdapps) endif() diff --git a/Modules/MatchPointRegistration/Testing/files.cmake b/Modules/MatchPointRegistration/Testing/files.cmake index 5efdc15458..b1268f0481 100644 --- a/Modules/MatchPointRegistration/Testing/files.cmake +++ b/Modules/MatchPointRegistration/Testing/files.cmake @@ -1,3 +1,4 @@ SET(MODULE_TESTS mitkTimeFramesRegistrationHelperTest.cpp + itkStitchImageFilterTest.cpp ) \ No newline at end of file diff --git a/Modules/MatchPointRegistration/Testing/itkStitchImageFilterTest.cpp b/Modules/MatchPointRegistration/Testing/itkStitchImageFilterTest.cpp new file mode 100644 index 0000000000..10c37646e5 --- /dev/null +++ b/Modules/MatchPointRegistration/Testing/itkStitchImageFilterTest.cpp @@ -0,0 +1,160 @@ +/*============================================================================ + +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. + +============================================================================*/ + +#include "mitkTestingMacros.h" +#include "mitkTestFixture.h" +#include "mitkTestDynamicImageGenerator.h" + +#include "itkStitchImageFilter.h" +#include "itkTranslationTransform.h" +#include "itkNearestNeighborInterpolateImageFunction.h" + +class itkStitchImageFilterTestSuite : public mitk::TestFixture +{ + CPPUNIT_TEST_SUITE(itkStitchImageFilterTestSuite); + // Test the append method + MITK_TEST(StitchWithNoTransformAndNoInterp); + MITK_TEST(StitchWithNoInterp); + MITK_TEST(Stitch); + CPPUNIT_TEST_SUITE_END(); + + using InputImageType = mitk::TestImageType; + using OutputImageType = itk::Image; +private: + using FilterType = itk::StitchImageFilter; + FilterType::Pointer m_Filter; + + InputImageType::Pointer m_Input1; + InputImageType::Pointer m_Input2; + InputImageType::Pointer m_Input3; + +public: + void setUp() override + { + InputImageType::PointType origin; + origin.Fill(0.); + InputImageType::SpacingType spacing; + spacing.Fill(1.); + spacing[1] = 5.; + + m_Filter = FilterType::New(); + m_Input1 = mitk::GenerateTestImage(1); + m_Input1->SetSpacing(spacing); + + m_Input2 = mitk::GenerateTestImage(10); + origin[1] = 10.; + m_Input2->SetOrigin(origin); + m_Input2->SetSpacing(spacing); + + m_Input3 = mitk::GenerateTestImage(100); + origin[1] = 20.; + m_Input3->SetOrigin(origin); + m_Input3->SetSpacing(spacing); + + FilterType::SizeType size = { 3, 9 }; + m_Filter->SetDefaultPixelValue(1000); + m_Filter->SetSize(size); + m_Filter->SetOutputSpacing(spacing); + } + + void tearDown() override + { + } + + bool CheckPixels(const OutputImageType* image, const std::vector& pixels) + { + bool result = true; + + itk::ImageRegionConstIteratorWithIndex iter(image, image->GetLargestPossibleRegion()); + auto refIter = pixels.begin(); + iter.GoToBegin(); + while (!iter.IsAtEnd()) + { + if (refIter == pixels.end()) + { + std::cerr << "Error image to check has a different pixel count then the reference pixel value vector."<SetInput(0, m_Input1); + m_Filter->SetInput(1, m_Input2); + m_Filter->SetInput(2, m_Input3); + + m_Filter->Update(); + auto output = m_Filter->GetOutput(); + + CPPUNIT_ASSERT(CheckPixels(output, {1, 2, 3, 4, 5, 6, 8.5, 14, 19.5, 40, 50, 60, 85, 140, 195, 400, 500, 600, 700, 800, 900, 1000, 1000, 1000, 1000, 1000, 1000})); + } + + void StitchWithNoInterp() + { + m_Filter->SetInput(0, m_Input1); + + using TranslationType = itk::TranslationTransform; + TranslationType::OutputVectorType offset; + offset[0] = 0.; + offset[1] = -5.; + auto translation1 = TranslationType::New(); + translation1->SetOffset(offset); + m_Filter->SetInput(1, m_Input2, translation1); + + offset[1] = -10.; + auto translation2 = TranslationType::New(); + translation2->SetOffset(offset); + m_Filter->SetInput(2, m_Input3, translation2); + + m_Filter->Update(); + auto output = m_Filter->GetOutput(); + + CPPUNIT_ASSERT(CheckPixels(output, { 1,2,3,4,5,6,7,8,9,10,20,30,40,50,60,70,80,90,100,200,300,400,500,600,700,800,900})); + } + + void Stitch() + { + using TranslationType = itk::TranslationTransform; + TranslationType::OutputVectorType offset; + offset[0] = 0; + offset[1] = -7.5; + auto translation1 = TranslationType::New(); + translation1->SetOffset(offset); + + offset[1] = -12.5; + auto translation2 = TranslationType::New(); + translation2->SetOffset(offset); + + m_Filter->SetInput(0, m_Input1); + m_Filter->SetInput(1, m_Input2, translation1, itk::NearestNeighborInterpolateImageFunction::New()); + m_Filter->SetInput(2, m_Input3, translation2); + + m_Filter->Update(); + auto output = m_Filter->GetOutput(); + + CPPUNIT_ASSERT(CheckPixels(output, { 1,2,3,4,5,6,7,8,9,10,20,30,40,50,60,70,80,90,100,200,300,250,350,450,550,650,750 })); + } + +}; + +MITK_TEST_SUITE_REGISTRATION(itkStitchImageFilter) diff --git a/Modules/MatchPointRegistration/cmdapps/CMakeLists.txt b/Modules/MatchPointRegistration/cmdapps/CMakeLists.txt new file mode 100644 index 0000000000..f02ab70bc2 --- /dev/null +++ b/Modules/MatchPointRegistration/cmdapps/CMakeLists.txt @@ -0,0 +1,33 @@ +option(BUILD_MatchPointCmdApps "Build commandline tools for the MatchPoint module" OFF) + +if(BUILD_MatchPointCmdApps OR MITK_BUILD_ALL_APPS) + + # needed include directories + include_directories( + ${CMAKE_CURRENT_SOURCE_DIR} + ${CMAKE_CURRENT_BINARY_DIR} + ) + # list of CmdApps + # if an app requires additional dependencies + # they are added after a "^^" and separated by "_" + set( cmdapps + StitchImagesMiniApp^^ + ) + + foreach(cmdapp ${cmdapps}) + # extract cmd name and dependencies + string(REPLACE "^^" "\\;" cmdapp_info ${cmdapp}) + set(cmdapp_info_list ${cmdapp_info}) + list(GET cmdapp_info_list 0 appname) + list(GET cmdapp_info_list 1 raw_dependencies) + string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") + set(dependencies_list ${dependencies}) + + mitkFunctionCreateCommandLineApp( + NAME ${appname} + DEPENDS MitkCore MitkMatchPointRegistration ${dependencies_list} + PACKAGE_DEPENDS ITK + ) + endforeach() + +endif(BUILD_MatchPointCmdApps OR MITK_BUILD_ALL_APPS) diff --git a/Modules/MatchPointRegistration/cmdapps/StitchImagesMiniApp.cpp b/Modules/MatchPointRegistration/cmdapps/StitchImagesMiniApp.cpp new file mode 100644 index 0000000000..ef8fc82bd2 --- /dev/null +++ b/Modules/MatchPointRegistration/cmdapps/StitchImagesMiniApp.cpp @@ -0,0 +1,226 @@ +/*============================================================================ + +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. + +============================================================================*/ + +// std includes +#include +#include + +// itk includes +#include "itksys/SystemTools.hxx" + +// CTK includes +#include "mitkCommandLineParser.h" + +// MITK includes +#include +#include +#include +#include +#include + +mitkCommandLineParser::StringContainerType inFilenames; +mitkCommandLineParser::StringContainerType regFilenames; +std::string outFileName; +std::string refGeometryFileName; + +std::vector images; +std::vector registrations; +mitk::BaseGeometry::Pointer refGeometry; +mitk::ImageMappingInterpolator::Type interpolatorType = mitk::ImageMappingInterpolator::Linear; +double paddingValue = 0; +itk::StitchStrategy stitchStratgy = itk::StitchStrategy::Mean; + +void setupParser(mitkCommandLineParser& parser) +{ + // set general information about your MiniApp + parser.setCategory("Mapping Tools"); + parser.setTitle("Stitch 3D Images"); + parser.setDescription("MiniApp that allows to map and stitch 3D images into a given output geometry."); + parser.setContributor("DKFZ MIC"); + //! [create parser] + + //! [add arguments] + // how should arguments be prefixed + parser.setArgumentPrefix("--", "-"); + // add each argument, unless specified otherwise each argument is optional + // see mitkCommandLineParser::addArgument for more information + parser.beginGroup("Required I/O parameters"); + parser.addArgument( + "inputs", "i", mitkCommandLineParser::StringList, "Input files", "Pathes to the input images that should be mapped and stitched", us::Any(), false, false, false, mitkCommandLineParser::Input); + parser.addArgument("output", + "o", + mitkCommandLineParser::File, + "Output file path", + "Path to the fused 3D+t image.", + us::Any(), + false, false, false, mitkCommandLineParser::Output); + parser.endGroup(); + + parser.beginGroup("Optional parameters"); + parser.addArgument("template", + "t", + mitkCommandLineParser::File, + "Output template image.", + "File path to an image that serves as template for the output geometry.", + us::Any(), + false, false, false, mitkCommandLineParser::Input); + parser.addArgument( + "registrations", "r", mitkCommandLineParser::StringList, "Registration files", "Pathes to the registrations that should be used to map the input images. If this parameter is not set, identity transforms are assumed. If this parameter is set, it must have the same number of entries then the parameter inputs. If you want to use and identity transform for a specific input, specify an empty string. The application assumes that inputs and registrations have the same order, so the n-th input should use thr n-th registration.", us::Any(), true, false, false, mitkCommandLineParser::Input); + parser.addArgument("interpolator", "n", mitkCommandLineParser::Int, "Interpolator type", "Interpolator used for mapping the images. Default: 2; allowed values: 1: Nearest Neighbour, 2: Linear, 3: BSpline 3, 4: WSinc Hamming, 5: WSinc Welch", us::Any(2), true); + parser.addArgument("strategy", "s", mitkCommandLineParser::Int, "Stitch strategy", "Strategy used for stitching the images. 0: Mean -> computes the mean value of all input images that cover an output pixel (default strategy). 1: BorderDistance -> Uses the input pixel that has the largest minimal distance to its image borders", us::Any(2), true); + parser.addArgument("padding", "p", mitkCommandLineParser::Float, "Padding value", "Value used for output voxels that are not covered by any input image.", us::Any(0.), true); + parser.addArgument("help", "h", mitkCommandLineParser::Bool, "Help:", "Show this help text"); + parser.endGroup(); + //! [add arguments] +} + +bool configureApplicationSettings(std::map parsedArgs) +{ + try + { + if (parsedArgs.size() == 0) + return false; + + inFilenames = us::any_cast(parsedArgs["inputs"]); + outFileName = us::any_cast(parsedArgs["output"]); + + if (parsedArgs.count("template")) + { + refGeometryFileName = us::any_cast(parsedArgs["template"]); + } + + if (parsedArgs.count("registrations")) + { + regFilenames = us::any_cast(parsedArgs["registrations"]); + } + else + { + regFilenames.resize(inFilenames.size()); + std::fill(regFilenames.begin(), regFilenames.end(), ""); + } + + if (parsedArgs.count("interpolator")) + { + auto interpolator = us::any_cast(parsedArgs["interpolator"]); + interpolatorType = static_cast(interpolator); + } + + if (parsedArgs.count("padding")) + { + paddingValue = us::any_cast(parsedArgs["padding"]); + } + + if (parsedArgs.count("strategy")) + { + auto temp = us::any_cast(parsedArgs["strategy"]); + stitchStratgy = static_cast(temp); + } + } + catch (...) + { + return false; + } + + return true; +} + +int main(int argc, char* argv[]) +{ + mitkCommandLineParser parser; + setupParser(parser); + + mitk::PreferenceListReaderOptionsFunctor readerFilterFunctor = mitk::PreferenceListReaderOptionsFunctor({ "MITK DICOM Reader v2 (autoselect)" }, { "" }); + + const std::map& parsedArgs = parser.parseArguments(argc, argv); + if (!configureApplicationSettings(parsedArgs)) + { + MITK_ERROR << "Command line arguments are invalid. To see the correct usage please call with -h or --help to show the help information."; + return EXIT_FAILURE; + }; + + // Show a help message + if (parsedArgs.count("help") || parsedArgs.count("h")) + { + std::cout << parser.helpText(); + return EXIT_SUCCESS; + } + + if(regFilenames.size() != inFilenames.size()) + { + MITK_ERROR << "Cannot stitch inputs. The number of specified registrations does not match the number of inputs."; + return EXIT_FAILURE; + } + + //! [do processing] + try + { + std::cout << "Load images:" << std::endl; + + unsigned int index = 0; + for (auto path : inFilenames) + { + std::cout << "#"<(path, &readerFilterFunctor); + images.push_back(image.GetPointer()); + if (regFilenames[index].empty()) + { + std::cout << " associated registration: identity" << std::endl; + registrations.push_back(mitk::GenerateIdentityRegistration3D().GetPointer()); + } + else + { + std::cout << " associated registration: " << regFilenames[index] << std::endl; + auto reg = mitk::IOUtil::Load(regFilenames[index]); + registrations.push_back(reg.GetPointer()); + } + ++index; + } + std::cout << "Reference image: " << refGeometryFileName << std::endl << std::endl; + auto refImage = mitk::IOUtil::Load(refGeometryFileName, &readerFilterFunctor); + if (refImage.IsNotNull()) + { + refGeometry = refImage->GetGeometry(); + } + std::cout << "Padding value: " << paddingValue << std::endl; + std::cout << "Stitch strategy: "; + if (itk::StitchStrategy::Mean == stitchStratgy) + { + std::cout << "Mean " << std::endl; + } + else + { + std::cout << "BorderDistance" << std::endl; + } + + std::cout << "Stitch the images ..." << std::endl; + + auto output = mitk::StitchImages(images, registrations, refGeometry,paddingValue,stitchStratgy,interpolatorType); + + std::cout << "Save output image: " << outFileName << std::endl; + + mitk::IOUtil::Save(output, outFileName); + + std::cout << "Processing finished." << std::endl; + + return EXIT_SUCCESS; + } + catch (const std::exception& e) + { + MITK_ERROR << e.what(); + return EXIT_FAILURE; + } + catch (...) + { + MITK_ERROR << "Unexpected error encountered."; + return EXIT_FAILURE; + } +} diff --git a/Modules/MatchPointRegistration/files.cmake b/Modules/MatchPointRegistration/files.cmake index 1b9c8ff7e1..0bf01d264c 100644 --- a/Modules/MatchPointRegistration/files.cmake +++ b/Modules/MatchPointRegistration/files.cmake @@ -1,65 +1,32 @@ +file(GLOB_RECURSE H_FILES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}/include/*.h") +file(GLOB_RECURSE TPP_FILES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}/include/*.tpp") + set(CPP_FILES mitkMAPRegistrationWrapper.cpp mitkMAPRegistrationWrapperObjectFactory.cpp mitkRegEvaluationObjectFactory.cpp mitkRegEvaluationObject.cpp Helper/mitkUIDHelper.cpp Helper/mitkMAPAlgorithmHelper.cpp Helper/mitkMaskedAlgorithmHelper.cpp Helper/mitkRegistrationHelper.cpp Helper/mitkImageMappingHelper.cpp + Helper/mitkImageStitchingHelper.cpp Helper/mitkPointSetMappingHelper.cpp Helper/mitkResultNodeGenerationHelper.cpp Helper/mitkTimeFramesRegistrationHelper.cpp Rendering/mitkRegistrationWrapperMapper2D.cpp Rendering/mitkRegistrationWrapperMapper3D.cpp Rendering/mitkRegistrationWrapperMapperBase.cpp Rendering/mitkRegEvaluationMapper2D.cpp Rendering/mitkRegVisStyleProperty.cpp Rendering/mitkRegVisDirectionProperty.cpp Rendering/mitkRegVisColorStyleProperty.cpp Rendering/mitkRegVisPropertyTags.cpp Rendering/mitkRegVisHelper.cpp Rendering/mitkRegEvalStyleProperty.cpp Rendering/mitkRegEvalWipeStyleProperty.cpp ) -set(H_FILES - mitkMatchPointPropertyTags.h - mitkMAPRegistrationWrapper.h - mitkMAPRegistrationWrapperObjectFactory.h - mitkRegEvaluationObjectFactory.h - mitkRegEvaluationObject.h - algorithms/mitkMultiModalAffineDefaultRegistrationAlgorithm.h - algorithms/mitkMultiModalRigidDefaultRegistrationAlgorithm.h - algorithms/mitkMultiModalTransDefaultRegistrationAlgorithm.h - algorithms/mitkFastSymmetricForcesDemonsMultiResDefaultRegistrationAlgorithm.h - algorithms/mitkLevelSetMotionMultiResDefaultRegistrationAlgorithm.h - algorithms/mitkRigidClosedFormPointsDefaultRegistrationAlgorithm.h - algorithms/mitkRigidICPDefaultRegistrationAlgorithm.h - Helper/mitkUIDHelper.h - Helper/mitkMAPAlgorithmHelper.h - Helper/mitkMaskedAlgorithmHelper.h - Helper/mitkRegistrationHelper.h - Helper/mitkImageMappingHelper.h - Helper/mitkPointSetMappingHelper.h - Helper/mitkResultNodeGenerationHelper.h - Helper/mitkTimeFramesRegistrationHelper.h - Rendering/mitkRegistrationWrapperMapper2D.h - Rendering/mitkRegistrationWrapperMapper3D.h - Rendering/mitkRegistrationWrapperMapperBase.h - Rendering/mitkRegVisStyleProperty.h - Rendering/mitkRegVisDirectionProperty.h - Rendering/mitkRegVisColorStyleProperty.h - Rendering/mitkRegVisPropertyTags.h - Rendering/mitkRegVisHelper.h - Rendering/mitkRegEvaluationMapper2D.h - Rendering/mitkRegEvalStyleProperty.h - Rendering/mitkRegEvalWipeStyleProperty.h -) - -set(TPP_FILES -) - set(MOC_H_FILES ) diff --git a/Modules/MatchPointRegistration/Helper/QmitkAlgorithmListModel.h b/Modules/MatchPointRegistration/include/QmitkAlgorithmListModel.h similarity index 100% rename from Modules/MatchPointRegistration/Helper/QmitkAlgorithmListModel.h rename to Modules/MatchPointRegistration/include/QmitkAlgorithmListModel.h diff --git a/Modules/MatchPointRegistration/Helper/QmitkMapPropertyDelegate.h b/Modules/MatchPointRegistration/include/QmitkMapPropertyDelegate.h similarity index 100% rename from Modules/MatchPointRegistration/Helper/QmitkMapPropertyDelegate.h rename to Modules/MatchPointRegistration/include/QmitkMapPropertyDelegate.h diff --git a/Modules/MatchPointRegistration/include/itkStitchImageFilter.h b/Modules/MatchPointRegistration/include/itkStitchImageFilter.h new file mode 100644 index 0000000000..5d5edce5d0 --- /dev/null +++ b/Modules/MatchPointRegistration/include/itkStitchImageFilter.h @@ -0,0 +1,330 @@ +/*============================================================================ + +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 itkStitchImageFilter_h +#define itkStitchImageFilter_h + +#include "itkFixedArray.h" +#include "itkTransform.h" +#include "itkImageRegionIterator.h" +#include "itkImageToImageFilter.h" +#include "itkLinearInterpolateImageFunction.h" +#include "itkSize.h" +#include "itkDefaultConvertPixelTraits.h" +#include "itkDataObjectDecorator.h" + + +namespace itk +{ + enum class StitchStrategy + { + Mean = 0, //use the mean value of all inputs that can provide a pixel vaule + BorderDistance = 1 //use the value that is largest minimal distance to its image borders (use e.g. if vaules tend to be not reliable at borders) + }; + + std::ostream& operator<< (std::ostream& os, const itk::StitchStrategy& strategy) + { + if (itk::StitchStrategy::Mean == strategy) + os << "Mean"; + else if (itk::StitchStrategy::BorderDistance == strategy) + os << "BorderDistance"; + else + os << "unkown"; + + return os; + }; + + /** \class StitchImageFilter + * \brief ITK filter that resamples/stitches multiple images into a given reference geometry. + * + * StitchImageFilter is similar to itk's ResampleImageFilter, but in difference to the last + * mentioned StitchImageFilter is able to resample multiple input images at once (with a transform + * for each input image). If multiple input images cover the output region the behavior depends on + * the StitchStragy: + * - Mean: a weighted sum of all voxels mapped input pixel values will be calculated. + * - BorderDistance: the voxels will be choosen that have the largest minimal distance to its own image borders. + * + * All other behaviors are similar to itk::ResampleImageFilter. See the filter's description for + * more details. + */ +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType = double, + typename TTransformPrecisionType = TInterpolatorPrecisionType> +class StitchImageFilter : + public ImageToImageFilter< TInputImage, TOutputImage > +{ +public: + /** Standard class typedefs. */ + typedef StitchImageFilter Self; + typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass; + typedef SmartPointer< Self > Pointer; + typedef SmartPointer< const Self > ConstPointer; + + typedef TInputImage InputImageType; + typedef TOutputImage OutputImageType; + typedef typename InputImageType::Pointer InputImagePointer; + typedef typename InputImageType::ConstPointer InputImageConstPointer; + typedef typename OutputImageType::Pointer OutputImagePointer; + typedef typename InputImageType::RegionType InputImageRegionType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(StitchImageFilter, ImageToImageFilter); + + /** Number of dimensions. */ + itkStaticConstMacro(ImageDimension, unsigned int, + TOutputImage::ImageDimension); + itkStaticConstMacro(InputImageDimension, unsigned int, + TInputImage::ImageDimension); + + /** base type for images of the current ImageDimension */ + typedef ImageBase< itkGetStaticConstMacro(ImageDimension) > ImageBaseType; + + /** + * Transform typedef. + */ + typedef Transform< TTransformPrecisionType, + itkGetStaticConstMacro(ImageDimension), + itkGetStaticConstMacro(ImageDimension) > TransformType; + typedef typename TransformType::ConstPointer TransformPointerType; + typedef DataObjectDecorator DecoratedTransformType; + typedef typename DecoratedTransformType::Pointer DecoratedTransformPointer; + + + /** Interpolator typedef. */ + typedef InterpolateImageFunction< InputImageType, + TInterpolatorPrecisionType > InterpolatorType; + typedef typename InterpolatorType::Pointer InterpolatorPointerType; + + typedef typename InterpolatorType::OutputType InterpolatorOutputType; + + typedef DefaultConvertPixelTraits< InterpolatorOutputType > InterpolatorConvertType; + + typedef typename InterpolatorConvertType::ComponentType ComponentType; + + typedef LinearInterpolateImageFunction< InputImageType, + TInterpolatorPrecisionType > LinearInterpolatorType; + typedef typename LinearInterpolatorType::Pointer + LinearInterpolatorPointerType; + + /** Image size typedef. */ + typedef Size< itkGetStaticConstMacro(ImageDimension) > SizeType; + + /** Image index typedef. */ + typedef typename TOutputImage::IndexType IndexType; + + /** Image point typedef. */ + typedef typename InterpolatorType::PointType PointType; + //typedef typename TOutputImage::PointType PointType; + + /** Image pixel value typedef. */ + typedef typename TOutputImage::PixelType PixelType; + typedef typename TInputImage::PixelType InputPixelType; + + typedef DefaultConvertPixelTraits PixelConvertType; + + typedef typename PixelConvertType::ComponentType PixelComponentType; + + /** Input pixel continuous index typdef */ + typedef ContinuousIndex< TTransformPrecisionType, ImageDimension > + ContinuousInputIndexType; + + /** Typedef to describe the output image region type. */ + typedef typename TOutputImage::RegionType OutputImageRegionType; + + /** Image spacing,origin and direction typedef */ + typedef typename TOutputImage::SpacingType SpacingType; + typedef typename TOutputImage::PointType OriginPointType; + typedef typename TOutputImage::DirectionType DirectionType; + + using Superclass::GetInput; + + /** Typedef the reference image type to be the ImageBase of the OutputImageType */ + typedef ImageBase ReferenceImageBaseType; + + using Superclass::SetInput; + void SetInput(const InputImageType* image) override; + void SetInput(unsigned int index, const InputImageType* image) override; + /** Convinience methods that allows setting of input image and its transform in + one call.*/ + virtual void SetInput(unsigned int index, const InputImageType* image, const TransformType* transform); + virtual void SetInput(unsigned int index, const InputImageType* image, const TransformType* transform, InterpolatorType* interpolator); + + const TransformType* GetTransform(unsigned int index) const; + + const InterpolatorType* GetInterpolator(unsigned int index) const; + + /** Get/Set the size of the output image. */ + itkSetMacro(Size, SizeType); + itkGetConstReferenceMacro(Size, SizeType); + + /** Get/Set the pixel value when a transformed pixel is outside of the + * image. The default default pixel value is 0. */ + itkSetMacro(DefaultPixelValue, PixelType); + itkGetConstReferenceMacro(DefaultPixelValue, PixelType); + + /** Set the output image spacing. */ + itkSetMacro(OutputSpacing, SpacingType); + virtual void SetOutputSpacing(const double *values); + + /** Get the output image spacing. */ + itkGetConstReferenceMacro(OutputSpacing, SpacingType); + + /** Set the output image origin. */ + itkSetMacro(OutputOrigin, OriginPointType); + virtual void SetOutputOrigin(const double *values); + + /** Get the output image origin. */ + itkGetConstReferenceMacro(OutputOrigin, OriginPointType); + + /** Set the output direciton cosine matrix. */ + itkSetMacro(OutputDirection, DirectionType); + itkGetConstReferenceMacro(OutputDirection, DirectionType); + + /** Helper method to set the output parameters based on this image. */ + void SetOutputParametersFromImage(const ImageBaseType *image); + + /** Set the start index of the output largest possible region. + * The default is an index of all zeros. */ + itkSetMacro(OutputStartIndex, IndexType); + + /** Get the start index of the output largest possible region. */ + itkGetConstReferenceMacro(OutputStartIndex, IndexType); + + /** Set a reference image to use to define the output information. + * By default, output information is specificed through the + * SetOutputSpacing, Origin, and Direction methods. Alternatively, + * this method can be used to specify an image from which to + * copy the information. UseReferenceImageOn must be set to utilize the + * reference image. */ + itkSetInputMacro(ReferenceImage, ReferenceImageBaseType); + + /** Get the reference image that is defining the output information. */ + itkGetInputMacro(ReferenceImage, ReferenceImageBaseType); + + /** Turn on/off whether a specified reference image should be used to define + * the output information. */ + itkSetMacro(UseReferenceImage, bool); + itkBooleanMacro(UseReferenceImage); + itkGetConstMacro(UseReferenceImage, bool); + + itkSetMacro(StitchStrategy, StitchStrategy); + itkGetConstMacro(StitchStrategy, StitchStrategy); + + /** StitchImageFilter produces an image which is a different size + * than its input. As such, it needs to provide an implementation + * for GenerateOutputInformation() in order to inform the pipeline + * execution model. The original documentation of this method is + * below. \sa ProcessObject::GenerateOutputInformaton() */ + virtual void GenerateOutputInformation() ITK_OVERRIDE; + + /** StitchImageFilter needs a different input requested region than + * the output requested region. As such, StitchImageFilter needs + * to provide an implementation for GenerateInputRequestedRegion() + * in order to inform the pipeline execution model. + * \sa ProcessObject::GenerateInputRequestedRegion() */ + virtual void GenerateInputRequestedRegion() ITK_OVERRIDE; + + /** Set up state of filter before multi-threading. + * InterpolatorType::SetInputImage is not thread-safe and hence + * has to be set up before ThreadedGenerateData */ + virtual void BeforeThreadedGenerateData() ITK_OVERRIDE; + + /** Set the state of the filter after multi-threading. */ + virtual void AfterThreadedGenerateData() ITK_OVERRIDE; + + /** Compute the Modified Time based on the changed components. */ + ModifiedTimeType GetMTime(void) const ITK_OVERRIDE; + +#ifdef ITK_USE_CONCEPT_CHECKING + // Begin concept checking + itkConceptMacro( OutputHasNumericTraitsCheck, + ( Concept::HasNumericTraits< PixelComponentType > ) ); + // End concept checking +#endif + +protected: + StitchImageFilter(); + ~StitchImageFilter() ITK_OVERRIDE {} + void PrintSelf(std::ostream & os, Indent indent) const ITK_OVERRIDE; + + /** Override VeriyInputInformation() since this filter's inputs do + * not need to occoupy the same physical space. + * + * \sa ProcessObject::VerifyInputInformation + */ + virtual void VerifyInputInformation() ITK_OVERRIDE { } + + /** StitchImageFilter can be implemented as a multithreaded filter. + * Therefore, this implementation provides a ThreadedGenerateData() + * routine which is called for each processing thread. The output + * image data is allocated automatically by the superclass prior + * to calling ThreadedGenerateData(). + * ThreadedGenerateData can only write to the portion of the output image + * specified by the parameter "outputRegionForThread" + * \sa ImageToImageFilter::ThreadedGenerateData(), + * ImageToImageFilter::GenerateData() */ + virtual void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, + ThreadIdType threadId) ITK_OVERRIDE; + + /** Cast pixel from interpolator output to PixelType. */ + virtual PixelType CastPixelWithBoundsChecking( const InterpolatorOutputType value, + const ComponentType minComponent, + const ComponentType maxComponent) const; + + void SetTransform(unsigned int index, const TransformType* transform); + + /** Helper that ensures that a transform is specified for every input image. + If a input image has no specified transforms, an identity transform will + be created and set as default.*/ + void EnsureTransforms(); + + /** Helper that ensures that an interpolator is specified for every input image. + If a input image has no specified interpolator, a linear interpolator will + be created and set as default.*/ + void EnsureInterpolators(); + + static std::string GetTransformInputName(unsigned int index); + +private: + ITK_DISALLOW_COPY_AND_ASSIGN(StitchImageFilter); + + typedef std::vector InputImageVectorType; + typedef std::map TransformMapType; + typedef std::map InterpolatorMapType; + + InputImageVectorType GetInputs(); + TransformMapType GetTransforms(); + + InterpolatorMapType m_Interpolators; // Image function for + // interpolation + PixelType m_DefaultPixelValue; // default pixel value + // if the point is + // outside the image + SizeType m_Size; // Size of the output image + SpacingType m_OutputSpacing; // output image spacing + OriginPointType m_OutputOrigin; // output image origin + DirectionType m_OutputDirection; // output image direction cosines + IndexType m_OutputStartIndex; // output image start index + bool m_UseReferenceImage; + StitchStrategy m_StitchStrategy; +}; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkStitchImageFilter.tpp" +#endif + +#endif diff --git a/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp b/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp new file mode 100644 index 0000000000..1ec8334ada --- /dev/null +++ b/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp @@ -0,0 +1,639 @@ +/*============================================================================ + +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 itkStitchImageFilter_hxx +#define itkStitchImageFilter_hxx + +#include "itkStitchImageFilter.h" +#include "itkObjectFactory.h" +#include "itkIdentityTransform.h" +#include "itkProgressReporter.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkImageScanlineIterator.h" +#include "itkSpecialCoordinatesImage.h" +#include "itkDefaultConvertPixelTraits.h" +#include "itkSimpleDataObjectDecorator.h" + +#include + +namespace itk +{ + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::StitchImageFilter() : + m_OutputSpacing( 1.0 ), + m_OutputOrigin( 0.0 ), + m_UseReferenceImage( false ), + m_StitchStrategy(StitchStrategy::Mean) +{ + + m_Size.Fill( 0 ); + m_OutputStartIndex.Fill( 0 ); + + m_OutputDirection.SetIdentity(); + + // Pipeline input configuration + + // implicit input index set: + // #1 "ReferenceImage" optional + Self::AddOptionalInputName("ReferenceImage"); + + m_DefaultPixelValue + = NumericTraits::ZeroValue( m_DefaultPixelValue ); +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > + void + StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > + ::SetInput(const InputImageType* image) +{ + this->SetInput(0, image, itk::IdentityTransform< TTransformPrecisionType, ImageDimension>::New().GetPointer(), LinearInterpolatorType::New().GetPointer()); +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +void +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::SetInput(unsigned int index, const InputImageType* image) +{ + this->SetInput(index, image, itk::IdentityTransform< TTransformPrecisionType, ImageDimension>::New().GetPointer(), LinearInterpolatorType::New().GetPointer()); +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +void +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::SetInput(unsigned int index, const InputImageType* image, const TransformType* transform) +{ + this->SetInput(index, image, transform, LinearInterpolatorType::New().GetPointer()); +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +void +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::SetInput(unsigned int index, const InputImageType* image, const TransformType* transform, InterpolatorType* interpolator) +{ + Superclass::SetInput(index, image); + m_Interpolators[image] = interpolator; + + this->SetTransform(index, transform); +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > + void + StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > + ::SetTransform(unsigned int index, const TransformType* transform) +{ + const auto transformName = this->GetTransformInputName(index); + typedef SimpleDataObjectDecorator< TransformPointerType > DecoratorType; + const DecoratorType* oldInput = itkDynamicCastInDebugMode< const DecoratorType* >(this->ProcessObject::GetInput(transformName)); + if (!oldInput || oldInput->Get() != transform) + { + typename DecoratorType::Pointer newInput = DecoratorType::New(); + // Process object is not const-correct so the const_cast is required here + newInput->Set(const_cast(transform)); + this->ProcessObject::SetInput(transformName, newInput); + } +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > + const typename StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType >::TransformType* + StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > + ::GetTransform(unsigned int index) const +{ + typedef SimpleDataObjectDecorator< TransformPointerType > DecoratorType; + const DecoratorType* input = itkDynamicCastInDebugMode< const DecoratorType* >(this->ProcessObject::GetInput(this->GetTransformInputName(index))); + + if (nullptr != input) + { + return input->Get(); + } + + return nullptr; +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > + const typename StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType >::InterpolatorType* + StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > + ::GetInterpolator(unsigned int index) const +{ + auto input = this->GetInput(index); + if (m_Interpolators.find(input) != std::end(m_Interpolators)) + { + return m_Interpolators[input]; + } + + return nullptr; +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +void +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::SetOutputSpacing(const double *spacing) +{ + SpacingType s; + for(unsigned int i = 0; i < TOutputImage::ImageDimension; ++i) + { + s[i] = static_cast< typename SpacingType::ValueType >(spacing[i]); + } + this->SetOutputSpacing(s); +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +void +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::SetOutputOrigin(const double *origin) +{ + OriginPointType p(origin); + + this->SetOutputOrigin(p); +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +void +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::SetOutputParametersFromImage(const ImageBaseType *image) +{ + this->SetOutputOrigin ( image->GetOrigin() ); + this->SetOutputSpacing ( image->GetSpacing() ); + this->SetOutputDirection ( image->GetDirection() ); + this->SetOutputStartIndex ( image->GetLargestPossibleRegion().GetIndex() ); + this->SetSize ( image->GetLargestPossibleRegion().GetSize() ); +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +void +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::BeforeThreadedGenerateData() +{ + this->EnsureInterpolators(); + this->EnsureTransforms(); + + for (const auto& interpolator : m_Interpolators) + { + interpolator.second->SetInputImage(interpolator.first); + } + + unsigned int nComponents + = DefaultConvertPixelTraits::GetNumberOfComponents( + m_DefaultPixelValue ); + + if (nComponents == 0) + { + PixelComponentType zeroComponent + = NumericTraits::ZeroValue( zeroComponent ); + nComponents = this->GetInput()->GetNumberOfComponentsPerPixel(); + NumericTraits::SetLength(m_DefaultPixelValue, nComponents ); + for (unsigned int n=0; n +void +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::AfterThreadedGenerateData() +{ + // Disconnect input image from the interpolator + for (auto& interpolator : m_Interpolators) + { + interpolator.second->SetInputImage(ITK_NULLPTR); + } +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +void +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, + ThreadIdType threadId) +{ + + if( outputRegionForThread.GetNumberOfPixels() == 0 ) + { + return; + } + + // Get the output pointers + OutputImageType* outputPtr = this->GetOutput(); + // Get this input pointers + InputImageVectorType inputs = this->GetInputs(); + TransformMapType transforms = this->GetTransforms(); + + std::map lowerIndices; + std::map upperIndices; + for (const auto& input : inputs) + { + const auto largestRegion = input->GetLargestPossibleRegion(); + lowerIndices[input] = largestRegion.GetIndex(); + upperIndices[input] = largestRegion.GetUpperIndex(); + } + + // Create an iterator that will walk the output region for this thread. + typedef ImageRegionIteratorWithIndex< OutputImageType > OutputIterator; + OutputIterator outIt(outputPtr, outputRegionForThread); + + // Define a few indices that will be used to translate from an input pixel + // to an output pixel + PointType outputPoint; // Coordinates of current output pixel + PointType inputPoint; // Coordinates of current input pixel + + ContinuousInputIndexType inputIndex; + + // Support for progress methods/callbacks + ProgressReporter progress(this, + threadId, + outputRegionForThread.GetNumberOfPixels()); + + // Min/max values of the output pixel type AND these values + // represented as the output type of the interpolator + const PixelComponentType minValue = NumericTraits< PixelComponentType >::NonpositiveMin(); + const PixelComponentType maxValue = NumericTraits< PixelComponentType >::max(); + + typedef typename InterpolatorType::OutputType OutputType; + const ComponentType minOutputValue = static_cast(minValue); + const ComponentType maxOutputValue = static_cast(maxValue); + + // Walk the output region + outIt.GoToBegin(); + + while (!outIt.IsAtEnd()) + { + // Determine the index of the current output pixel + outputPtr->TransformIndexToPhysicalPoint(outIt.GetIndex(), outputPoint); + + std::vector pixvals; + std::vector pixDistance; + + for (const auto& input : inputs) + { + // Compute corresponding input pixel position + inputPoint = transforms[input]->TransformPoint(outputPoint); + const bool isInsideInput = input->TransformPhysicalPointToContinuousIndex(inputPoint, inputIndex); + + // Evaluate input at right position and copy to the output + if (m_Interpolators[input]->IsInsideBuffer(inputIndex) && isInsideInput) + { + OutputType value = m_Interpolators[input]->EvaluateAtContinuousIndex(inputIndex); + pixvals.emplace_back(this->CastPixelWithBoundsChecking(value, minOutputValue, maxOutputValue)); + + ContinuousInputIndexType indexDistance; + const auto spacing = input->GetSpacing(); + + double minBorderDistance = std::numeric_limits::max(); + for (unsigned int i = 0; i < ImageDimension; ++i) + { + minBorderDistance = std::min(minBorderDistance, std::min(std::abs(lowerIndices[input][i] - inputIndex[i]) * spacing[i], std::abs(upperIndices[input][i] - inputIndex[i]) * spacing[i])); + } + pixDistance.emplace_back(minBorderDistance); + } + } + + if (!pixvals.empty()) + { //at least one input provided a value + if (StitchStrategy::Mean == m_StitchStrategy) + { + double sum = std::accumulate(pixvals.begin(), pixvals.end(), 0.0); + outIt.Set(sum / pixvals.size()); + } + else + { + auto finding = std::max_element(pixDistance.begin(), pixDistance.end()); + outIt.Set(pixvals[std::distance(pixDistance.begin(), finding)]); + } + } + else + { + outIt.Set(m_DefaultPixelValue); // default background value + } + + progress.CompletedPixel(); + ++outIt; + } +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +typename StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::PixelType +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::CastPixelWithBoundsChecking(const InterpolatorOutputType value, + const ComponentType minComponent, + const ComponentType maxComponent ) const +{ + const unsigned int nComponents = InterpolatorConvertType::GetNumberOfComponents(value); + PixelType outputValue; + + NumericTraits::SetLength( outputValue, nComponents ); + + for (unsigned int n = 0; n < nComponents; n++) + { + ComponentType component = InterpolatorConvertType::GetNthComponent( n, value ); + + if ( component < minComponent ) + { + PixelConvertType::SetNthComponent( n, outputValue, static_cast( minComponent ) ); + } + else if ( component > maxComponent ) + { + PixelConvertType::SetNthComponent( n, outputValue, static_cast( maxComponent ) ); + } + else + { + PixelConvertType::SetNthComponent(n, outputValue, + static_cast( component ) ); + } + } + + return outputValue; +} + +template +typename StitchImageFilter::InputImageVectorType +StitchImageFilter +::GetInputs() +{ + InputImageVectorType inputs; + for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); ++i) + { + auto input = this->GetInput(i); + if (nullptr != input) + { + inputs.push_back(input); + } + } + return inputs; +} + +template +typename StitchImageFilter::TransformMapType +StitchImageFilter +::GetTransforms() +{ + TransformMapType transforms; + for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); ++i) + { + auto input = this->GetInput(i); + auto transform = this->GetTransform(i); + transforms[input] = transform; + } + return transforms; +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +void +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::GenerateInputRequestedRegion() +{ + // Call the superclass' implementation of this method + Superclass::GenerateInputRequestedRegion(); + + if ( !this->GetInput() ) + { + return; + } + + // Get pointers to the input + auto inputs = this->GetInputs(); + + for (auto& input : inputs) + { + InputImagePointer inputPtr = + const_cast(input); + // Determining the actual input region is non-trivial, especially + // when we cannot assume anything about the transform being used. + // So we do the easy thing and request the entire input image. + // + inputPtr->SetRequestedRegionToLargestPossibleRegion(); + } +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +void +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::GenerateOutputInformation() +{ + // Call the superclass' implementation of this method + Superclass::GenerateOutputInformation(); + + // Get pointers to the input and output + OutputImageType *outputPtr = this->GetOutput(); + if ( !outputPtr ) + { + return; + } + + const ReferenceImageBaseType *referenceImage = this->GetReferenceImage(); + + // Set the size of the output region + if ( m_UseReferenceImage && referenceImage ) + { + outputPtr->SetLargestPossibleRegion( + referenceImage->GetLargestPossibleRegion() ); + } + else + { + typename TOutputImage::RegionType outputLargestPossibleRegion; + outputLargestPossibleRegion.SetSize(m_Size); + outputLargestPossibleRegion.SetIndex(m_OutputStartIndex); + outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion); + } + + // Set spacing and origin + if ( m_UseReferenceImage && referenceImage ) + { + outputPtr->SetSpacing( referenceImage->GetSpacing() ); + outputPtr->SetOrigin( referenceImage->GetOrigin() ); + outputPtr->SetDirection( referenceImage->GetDirection() ); + } + else + { + outputPtr->SetSpacing(m_OutputSpacing); + outputPtr->SetOrigin(m_OutputOrigin); + outputPtr->SetDirection(m_OutputDirection); + } +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +ModifiedTimeType +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::GetMTime(void) const +{ + ModifiedTimeType latestTime = Object::GetMTime(); + + for (const auto& interpolator : m_Interpolators) + { + if (interpolator.second.GetPointer()) + { + if (latestTime < interpolator.second->GetMTime()) + { + latestTime = interpolator.second->GetMTime(); + } + } + } + + return latestTime; +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +void +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + + os << indent << "DefaultPixelValue: " + << static_cast< typename NumericTraits< PixelType >::PrintType > + ( m_DefaultPixelValue ) + << std::endl; + os << indent << "Size: " << m_Size << std::endl; + os << indent << "OutputStartIndex: " << m_OutputStartIndex << std::endl; + os << indent << "OutputSpacing: " << m_OutputSpacing << std::endl; + os << indent << "OutputOrigin: " << m_OutputOrigin << std::endl; + os << indent << "OutputDirection: " << m_OutputDirection << std::endl; + for (const auto& interpolator : m_Interpolators) + { + os << indent << "Interpolator: " << interpolator.second.GetPointer() << std::endl; + } + os << indent << "UseReferenceImage: " << ( m_UseReferenceImage ? "On" : "Off" ) + << std::endl; +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +void +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::EnsureTransforms() +{ + const auto inputCount = this->GetNumberOfIndexedInputs(); + for (unsigned int i = 0; i < inputCount; ++i) + { + auto input = this->GetInput(i); + + if (nullptr == input) + { + itkExceptionMacro(<< "Nth input image is not set (n: " << i << ")."); + } + + auto transform = this->GetTransform(i); + if (nullptr == transform) + { + this->SetTransform(i, itk::IdentityTransform< TTransformPrecisionType, ImageDimension>::New().GetPointer()); + } + } +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +void +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::EnsureInterpolators() +{ + const auto inputCount = this->GetNumberOfIndexedInputs(); + InterpolatorMapType newInterpolatorMap; + + for (unsigned int i = 0; i < inputCount; ++i) + { + auto input = this->GetInput(i); + + if (nullptr == input) + { + itkExceptionMacro(<< "Nth input image is not set (n: " << i << ")."); + } + + if (m_Interpolators[input].IsNull()) + { + newInterpolatorMap[input] = LinearInterpolatorType::New().GetPointer(); + } + else + { + newInterpolatorMap[input] = m_Interpolators[input]; + } + } + m_Interpolators = newInterpolatorMap; +} + +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > +std::string +StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > +::GetTransformInputName(unsigned int index) +{ + return "transform_" + std::to_string(index); +} + +} // end namespace itk + +#endif diff --git a/Modules/MatchPointRegistration/Helper/mitkImageMappingHelper.h b/Modules/MatchPointRegistration/include/mitkImageMappingHelper.h similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkImageMappingHelper.h rename to Modules/MatchPointRegistration/include/mitkImageMappingHelper.h diff --git a/Modules/MatchPointRegistration/include/mitkImageStitchingHelper.h b/Modules/MatchPointRegistration/include/mitkImageStitchingHelper.h new file mode 100644 index 0000000000..8d4b12f3ec --- /dev/null +++ b/Modules/MatchPointRegistration/include/mitkImageStitchingHelper.h @@ -0,0 +1,67 @@ +/*============================================================================ + +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 MITK_IMAGE_STITCHING_HELPER_H +#define MITK_IMAGE_STITCHING_HELPER_H + +#include "mapRegistrationBase.h" +#include "mitkImage.h" +#include "mitkGeometry3D.h" + +#include "mitkMAPRegistrationWrapper.h" +#include "mitkImageMappingHelper.h" + +#include + +#include "MitkMatchPointRegistrationExports.h" + +namespace mitk +{ + /**Helper that stitches a given vector of input images + * @param inputs vector of input images that should be stitched. + * @param registrations vector of registrations that should be used for mapping of the inputs before stitching. + * the method assumes that order of registrations is the same as the order of inputs, thus for the n-th input + * the n-th registration will be used. + * @param resultGeometry Pointer to the Geometry object that specifies the grid of the result image. + * @param paddingValue Indicates the value that should be used if an out of input error occurs (and throwOnOutOfInputAreaError is false). + * @param interpolatorType Indicates the type of interpolation strategy that should be used. + * @param stitchStrategy Strategy used if more than one input can contribute. for more details see the documentation of itk::StitchStrategy. + * @pre inputs must not be empty and contain valid instances + * @pre registration must have same size as inputs and contain valid instances. + * @pre Dimensionality of the registrations must match with the inputs + * @pre resultGeometry must be valid. + * @remark The helper currently only supports 3D images. + * @result Pointer to the resulting mapped image.h*/ + MITKMATCHPOINTREGISTRATION_EXPORT Image::Pointer StitchImages(std::vector inputs, + std::vector<::map::core::RegistrationBase::ConstPointer> registrations, + const BaseGeometry* resultGeometry, + const double& paddingValue = 0, itk::StitchStrategy stitchStrategy = itk::StitchStrategy::Mean, + mitk::ImageMappingInterpolator::Type interpolatorType = mitk::ImageMappingInterpolator::Linear); + + MITKMATCHPOINTREGISTRATION_EXPORT Image::Pointer StitchImages(std::vector inputs, + std::vector registrations, + const BaseGeometry* resultGeometry, + const double& paddingValue = 0, itk::StitchStrategy stitchStrategy = itk::StitchStrategy::Mean, + mitk::ImageMappingInterpolator::Type interpolatorType = mitk::ImageMappingInterpolator::Linear); + + /**@overload + * Convinience version that uses identity transforms form the registrations. + */ + MITKMATCHPOINTREGISTRATION_EXPORT Image::Pointer StitchImages(std::vector inputs, + const BaseGeometry* resultGeometry, + const double& paddingValue = 0, itk::StitchStrategy stitchStrategy = itk::StitchStrategy::Mean, + mitk::ImageMappingInterpolator::Type interpolatorType = mitk::ImageMappingInterpolator::Linear); + +} + +#endif diff --git a/Modules/MatchPointRegistration/Helper/mitkMAPAlgorithmHelper.h b/Modules/MatchPointRegistration/include/mitkMAPAlgorithmHelper.h similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkMAPAlgorithmHelper.h rename to Modules/MatchPointRegistration/include/mitkMAPAlgorithmHelper.h diff --git a/Modules/MatchPointRegistration/mitkMAPRegistrationWrapper.h b/Modules/MatchPointRegistration/include/mitkMAPRegistrationWrapper.h similarity index 100% rename from Modules/MatchPointRegistration/mitkMAPRegistrationWrapper.h rename to Modules/MatchPointRegistration/include/mitkMAPRegistrationWrapper.h diff --git a/Modules/MatchPointRegistration/mitkMAPRegistrationWrapperObjectFactory.h b/Modules/MatchPointRegistration/include/mitkMAPRegistrationWrapperObjectFactory.h similarity index 100% rename from Modules/MatchPointRegistration/mitkMAPRegistrationWrapperObjectFactory.h rename to Modules/MatchPointRegistration/include/mitkMAPRegistrationWrapperObjectFactory.h diff --git a/Modules/MatchPointRegistration/Helper/mitkMaskedAlgorithmHelper.h b/Modules/MatchPointRegistration/include/mitkMaskedAlgorithmHelper.h similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkMaskedAlgorithmHelper.h rename to Modules/MatchPointRegistration/include/mitkMaskedAlgorithmHelper.h diff --git a/Modules/MatchPointRegistration/mitkMatchPointPropertyTags.h b/Modules/MatchPointRegistration/include/mitkMatchPointPropertyTags.h similarity index 100% rename from Modules/MatchPointRegistration/mitkMatchPointPropertyTags.h rename to Modules/MatchPointRegistration/include/mitkMatchPointPropertyTags.h diff --git a/Modules/MatchPointRegistration/Helper/mitkPointSetMappingHelper.h b/Modules/MatchPointRegistration/include/mitkPointSetMappingHelper.h similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkPointSetMappingHelper.h rename to Modules/MatchPointRegistration/include/mitkPointSetMappingHelper.h diff --git a/Modules/MatchPointRegistration/Helper/mitkQMAPAlgorithmModel.h b/Modules/MatchPointRegistration/include/mitkQMAPAlgorithmModel.h similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkQMAPAlgorithmModel.h rename to Modules/MatchPointRegistration/include/mitkQMAPAlgorithmModel.h diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegEvalStyleProperty.h b/Modules/MatchPointRegistration/include/mitkRegEvalStyleProperty.h similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegEvalStyleProperty.h rename to Modules/MatchPointRegistration/include/mitkRegEvalStyleProperty.h diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegEvalWipeStyleProperty.h b/Modules/MatchPointRegistration/include/mitkRegEvalWipeStyleProperty.h similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegEvalWipeStyleProperty.h rename to Modules/MatchPointRegistration/include/mitkRegEvalWipeStyleProperty.h diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegEvaluationMapper2D.h b/Modules/MatchPointRegistration/include/mitkRegEvaluationMapper2D.h similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegEvaluationMapper2D.h rename to Modules/MatchPointRegistration/include/mitkRegEvaluationMapper2D.h diff --git a/Modules/MatchPointRegistration/mitkRegEvaluationObject.h b/Modules/MatchPointRegistration/include/mitkRegEvaluationObject.h similarity index 100% rename from Modules/MatchPointRegistration/mitkRegEvaluationObject.h rename to Modules/MatchPointRegistration/include/mitkRegEvaluationObject.h diff --git a/Modules/MatchPointRegistration/mitkRegEvaluationObjectFactory.h b/Modules/MatchPointRegistration/include/mitkRegEvaluationObjectFactory.h similarity index 100% rename from Modules/MatchPointRegistration/mitkRegEvaluationObjectFactory.h rename to Modules/MatchPointRegistration/include/mitkRegEvaluationObjectFactory.h diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegVisColorStyleProperty.h b/Modules/MatchPointRegistration/include/mitkRegVisColorStyleProperty.h similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegVisColorStyleProperty.h rename to Modules/MatchPointRegistration/include/mitkRegVisColorStyleProperty.h diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegVisDirectionProperty.h b/Modules/MatchPointRegistration/include/mitkRegVisDirectionProperty.h similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegVisDirectionProperty.h rename to Modules/MatchPointRegistration/include/mitkRegVisDirectionProperty.h diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegVisHelper.h b/Modules/MatchPointRegistration/include/mitkRegVisHelper.h similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegVisHelper.h rename to Modules/MatchPointRegistration/include/mitkRegVisHelper.h diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegVisPropertyTags.h b/Modules/MatchPointRegistration/include/mitkRegVisPropertyTags.h similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegVisPropertyTags.h rename to Modules/MatchPointRegistration/include/mitkRegVisPropertyTags.h diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegVisStyleProperty.h b/Modules/MatchPointRegistration/include/mitkRegVisStyleProperty.h similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegVisStyleProperty.h rename to Modules/MatchPointRegistration/include/mitkRegVisStyleProperty.h diff --git a/Modules/MatchPointRegistration/Helper/mitkRegistrationHelper.h b/Modules/MatchPointRegistration/include/mitkRegistrationHelper.h similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkRegistrationHelper.h rename to Modules/MatchPointRegistration/include/mitkRegistrationHelper.h diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegistrationWrapperMapper2D.h b/Modules/MatchPointRegistration/include/mitkRegistrationWrapperMapper2D.h similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegistrationWrapperMapper2D.h rename to Modules/MatchPointRegistration/include/mitkRegistrationWrapperMapper2D.h diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegistrationWrapperMapper3D.h b/Modules/MatchPointRegistration/include/mitkRegistrationWrapperMapper3D.h similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegistrationWrapperMapper3D.h rename to Modules/MatchPointRegistration/include/mitkRegistrationWrapperMapper3D.h diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegistrationWrapperMapperBase.h b/Modules/MatchPointRegistration/include/mitkRegistrationWrapperMapperBase.h similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegistrationWrapperMapperBase.h rename to Modules/MatchPointRegistration/include/mitkRegistrationWrapperMapperBase.h diff --git a/Modules/MatchPointRegistration/Helper/mitkResultNodeGenerationHelper.h b/Modules/MatchPointRegistration/include/mitkResultNodeGenerationHelper.h similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkResultNodeGenerationHelper.h rename to Modules/MatchPointRegistration/include/mitkResultNodeGenerationHelper.h diff --git a/Modules/MatchPointRegistration/Helper/mitkTimeFramesRegistrationHelper.h b/Modules/MatchPointRegistration/include/mitkTimeFramesRegistrationHelper.h similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkTimeFramesRegistrationHelper.h rename to Modules/MatchPointRegistration/include/mitkTimeFramesRegistrationHelper.h diff --git a/Modules/MatchPointRegistration/Helper/mitkUIDHelper.h b/Modules/MatchPointRegistration/include/mitkUIDHelper.h similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkUIDHelper.h rename to Modules/MatchPointRegistration/include/mitkUIDHelper.h diff --git a/Modules/MatchPointRegistration/Helper/QmitkAlgorithmListModel.cpp b/Modules/MatchPointRegistration/src/Helper/QmitkAlgorithmListModel.cpp similarity index 100% rename from Modules/MatchPointRegistration/Helper/QmitkAlgorithmListModel.cpp rename to Modules/MatchPointRegistration/src/Helper/QmitkAlgorithmListModel.cpp diff --git a/Modules/MatchPointRegistration/Helper/QmitkMapPropertyDelegate.cpp b/Modules/MatchPointRegistration/src/Helper/QmitkMapPropertyDelegate.cpp similarity index 100% rename from Modules/MatchPointRegistration/Helper/QmitkMapPropertyDelegate.cpp rename to Modules/MatchPointRegistration/src/Helper/QmitkMapPropertyDelegate.cpp diff --git a/Modules/MatchPointRegistration/Helper/mitkImageMappingHelper.cpp b/Modules/MatchPointRegistration/src/Helper/mitkImageMappingHelper.cpp similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkImageMappingHelper.cpp rename to Modules/MatchPointRegistration/src/Helper/mitkImageMappingHelper.cpp diff --git a/Modules/MatchPointRegistration/src/Helper/mitkImageStitchingHelper.cpp b/Modules/MatchPointRegistration/src/Helper/mitkImageStitchingHelper.cpp new file mode 100644 index 0000000000..e0fcb5f46a --- /dev/null +++ b/Modules/MatchPointRegistration/src/Helper/mitkImageStitchingHelper.cpp @@ -0,0 +1,229 @@ +/*============================================================================ + +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. + +============================================================================*/ + +#include "mitkImageStitchingHelper.h" + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "mapRegistration.h" + +#include "mitkRegistrationHelper.h" + +template +typename ::itk::InterpolateImageFunction< TImage >::Pointer generateInterpolator(mitk::ImageMappingInterpolator::Type interpolatorType) +{ + typedef ::itk::InterpolateImageFunction< TImage > BaseInterpolatorType; + typename BaseInterpolatorType::Pointer result; + + switch (interpolatorType) + { + case mitk::ImageMappingInterpolator::NearestNeighbor: + { + result = ::itk::NearestNeighborInterpolateImageFunction::New(); + break; + } + case mitk::ImageMappingInterpolator::BSpline_3: + { + typename ::itk::BSplineInterpolateImageFunction::Pointer spInterpolator = ::itk::BSplineInterpolateImageFunction::New(); + spInterpolator->SetSplineOrder(3); + result = spInterpolator; + break; + } + case mitk::ImageMappingInterpolator::WSinc_Hamming: + { + result = ::itk::WindowedSincInterpolateImageFunction::New(); + break; + } + case mitk::ImageMappingInterpolator::WSinc_Welch: + { + result = ::itk::WindowedSincInterpolateImageFunction >::New(); + break; + } + default: + { + result = ::itk::LinearInterpolateImageFunction::New(); + break; + } + + } + + return result; +}; + +template +void doMITKStitching(const ::itk::Image* /*input1*/, + mitk::Image::Pointer& result, + std::vector inputs, + std::vector<::map::core::RegistrationBase::ConstPointer> registrations, + const mitk::BaseGeometry* resultGeometry, + const double& paddingValue, itk::StitchStrategy stitchStrategy, mitk::ImageMappingInterpolator::Type interpolatorType) +{ + using ConcreteRegistrationType = ::map::core::Registration; + using ItkImageType = itk::Image; + + using StitchingFilterType = ::itk::StitchImageFilter; + auto stitcher = StitchingFilterType::New(); + + stitcher->SetDefaultPixelValue(paddingValue); + + stitcher->SetOutputOrigin(resultGeometry->GetOrigin()); + + const auto spacing = resultGeometry->GetSpacing(); + stitcher->SetOutputSpacing(spacing); + + typename StitchingFilterType::DirectionType itkDirection; + const auto mitkDirection = resultGeometry->GetIndexToWorldTransform()->GetMatrix(); + for (unsigned int i = 0; i < VImageDimension; ++i) + { + for (unsigned int j = 0; j < VImageDimension; ++j) + { + itkDirection[i][j] = mitkDirection[i][j] / spacing[j]; + } + } + stitcher->SetOutputDirection(itkDirection); + + typename ItkImageType::SizeType size; + size[0] = resultGeometry->GetExtent(0); + size[1] = resultGeometry->GetExtent(1); + size[2] = resultGeometry->GetExtent(2); + stitcher->SetSize(size); + stitcher->SetNumberOfThreads(1); + stitcher->SetStitchStrategy(stitchStrategy); + + auto inputIter = inputs.begin(); + auto regIter = registrations.begin(); + unsigned int index = 0; + + while (inputIter != inputs.end()) + { + auto itkInput = mitk::ImageToItkImage(*inputIter); + + auto castedReg = dynamic_cast(regIter->GetPointer()); + + auto kernel = dynamic_cast* >(&(castedReg->getInverseMapping())); + if (nullptr == kernel) + { + mitkThrow() << "Cannot stitch images. At least passed registration object #"<SetInput(index, itkInput, kernel->getTransformModel(), generateInterpolator< ::itk::Image >(interpolatorType)); + ++inputIter; + ++regIter; + ++index; + } + + stitcher->Update(); + mitk::CastToMitkImage<>(stitcher->GetOutput(),result); +} + +mitk::Image::Pointer +mitk::StitchImages(std::vector inputs, + std::vector<::map::core::RegistrationBase::ConstPointer> registrations, + const BaseGeometry* resultGeometry, + const double& paddingValue, itk::StitchStrategy stitchStrategy, + mitk::ImageMappingInterpolator::Type interpolatorType) +{ + if (inputs.size() != registrations.size()) + { + mitkThrow() << "Cannot stitch images. Passed inputs vector and registrations vector have different sizes."; + } + + if (inputs.empty()) + { + mitkThrow() << "Cannot stitch images. No input images are defined."; + } + + auto inputDim = inputs.front()->GetDimension(); + auto inputPixelType = inputs.front()->GetPixelType(); + + for (const auto& input : inputs) + { + if (input->GetDimension() != inputDim) + { + mitkThrow() << "Cannot stitch images. Images have different dimensions. Dimeonsion of first input: " << inputDim << "; wrong dimension: " << input->GetDimension(); + } + if (input->GetPixelType() != inputPixelType) + { + mitkThrow() << "Cannot stitch images. Input images have different pixeltypes. The current implementation does only support stitching of images with same pixel type. Dimeonsion of first input: " << inputPixelType.GetTypeAsString() << "; wrong dimension: " << input->GetPixelType().GetTypeAsString(); + } + if (input->GetTimeSteps() > 1) + { + mitkThrow() << "Cannot stitch dynamic images. At least one input image has multiple time steps."; + } + } + + for (const auto& reg : registrations) + { + if (reg->getMovingDimensions() != inputDim) + { + mitkThrow() << "Cannot stitch images. At least one registration has a different moving dimension then the inputs. Dimeonsion of inputs: " << inputDim << "; wrong dimension: " << reg->getMovingDimensions(); + } + if (reg->getTargetDimensions() != inputDim) + { + mitkThrow() << "Cannot stitch images. At least one registration has a different target dimension then the inputs. Dimeonsion of inputs: " << inputDim << "; wrong dimension: " << reg->getTargetDimensions(); + } + } + + Image::Pointer result; + + AccessFixedDimensionByItk_n(inputs.front(), doMITKStitching, 3, (result, inputs, registrations, resultGeometry, paddingValue, stitchStrategy, interpolatorType)); + + return result; +} + +mitk::Image::Pointer +mitk::StitchImages(std::vector inputs, + std::vector registrations, + const BaseGeometry* resultGeometry, + const double& paddingValue, itk::StitchStrategy stitchStrategy, + mitk::ImageMappingInterpolator::Type interpolatorType) +{ + + std::vector<::map::core::RegistrationBase::ConstPointer> unwrappedRegs; + for (const auto& reg : registrations) + { + if (!reg) + { + mitkThrow() << "Cannot stitch images. At least one passed registration wrapper pointer is nullptr."; + } + unwrappedRegs.push_back(reg->GetRegistration()); + } + + Image::Pointer result = StitchImages(inputs, unwrappedRegs, resultGeometry, paddingValue, stitchStrategy, interpolatorType); + return result; +} + +mitk::Image::Pointer +mitk::StitchImages(std::vector inputs, + const BaseGeometry* resultGeometry, + const double& paddingValue, itk::StitchStrategy stitchStrategy, + mitk::ImageMappingInterpolator::Type interpolatorType) +{ + auto defaultReg = GenerateIdentityRegistration3D(); + std::vector<::map::core::RegistrationBase::ConstPointer> defaultRegs; + defaultRegs.resize(inputs.size()); + std::fill(defaultRegs.begin(), defaultRegs.end(), defaultReg->GetRegistration()); + + Image::Pointer result = StitchImages(inputs, defaultRegs, resultGeometry, paddingValue, stitchStrategy, interpolatorType); + return result; +} + diff --git a/Modules/MatchPointRegistration/Helper/mitkMAPAlgorithmHelper.cpp b/Modules/MatchPointRegistration/src/Helper/mitkMAPAlgorithmHelper.cpp similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkMAPAlgorithmHelper.cpp rename to Modules/MatchPointRegistration/src/Helper/mitkMAPAlgorithmHelper.cpp diff --git a/Modules/MatchPointRegistration/Helper/mitkMaskedAlgorithmHelper.cpp b/Modules/MatchPointRegistration/src/Helper/mitkMaskedAlgorithmHelper.cpp similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkMaskedAlgorithmHelper.cpp rename to Modules/MatchPointRegistration/src/Helper/mitkMaskedAlgorithmHelper.cpp diff --git a/Modules/MatchPointRegistration/Helper/mitkPointSetMappingHelper.cpp b/Modules/MatchPointRegistration/src/Helper/mitkPointSetMappingHelper.cpp similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkPointSetMappingHelper.cpp rename to Modules/MatchPointRegistration/src/Helper/mitkPointSetMappingHelper.cpp diff --git a/Modules/MatchPointRegistration/Helper/mitkQMAPAlgorithmModel.cpp b/Modules/MatchPointRegistration/src/Helper/mitkQMAPAlgorithmModel.cpp similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkQMAPAlgorithmModel.cpp rename to Modules/MatchPointRegistration/src/Helper/mitkQMAPAlgorithmModel.cpp diff --git a/Modules/MatchPointRegistration/Helper/mitkRegistrationHelper.cpp b/Modules/MatchPointRegistration/src/Helper/mitkRegistrationHelper.cpp similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkRegistrationHelper.cpp rename to Modules/MatchPointRegistration/src/Helper/mitkRegistrationHelper.cpp diff --git a/Modules/MatchPointRegistration/Helper/mitkResultNodeGenerationHelper.cpp b/Modules/MatchPointRegistration/src/Helper/mitkResultNodeGenerationHelper.cpp similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkResultNodeGenerationHelper.cpp rename to Modules/MatchPointRegistration/src/Helper/mitkResultNodeGenerationHelper.cpp diff --git a/Modules/MatchPointRegistration/Helper/mitkTimeFramesRegistrationHelper.cpp b/Modules/MatchPointRegistration/src/Helper/mitkTimeFramesRegistrationHelper.cpp similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkTimeFramesRegistrationHelper.cpp rename to Modules/MatchPointRegistration/src/Helper/mitkTimeFramesRegistrationHelper.cpp diff --git a/Modules/MatchPointRegistration/Helper/mitkUIDHelper.cpp b/Modules/MatchPointRegistration/src/Helper/mitkUIDHelper.cpp similarity index 100% rename from Modules/MatchPointRegistration/Helper/mitkUIDHelper.cpp rename to Modules/MatchPointRegistration/src/Helper/mitkUIDHelper.cpp diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegEvalStyleProperty.cpp b/Modules/MatchPointRegistration/src/Rendering/mitkRegEvalStyleProperty.cpp similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegEvalStyleProperty.cpp rename to Modules/MatchPointRegistration/src/Rendering/mitkRegEvalStyleProperty.cpp diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegEvalWipeStyleProperty.cpp b/Modules/MatchPointRegistration/src/Rendering/mitkRegEvalWipeStyleProperty.cpp similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegEvalWipeStyleProperty.cpp rename to Modules/MatchPointRegistration/src/Rendering/mitkRegEvalWipeStyleProperty.cpp diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegEvaluationMapper2D.cpp b/Modules/MatchPointRegistration/src/Rendering/mitkRegEvaluationMapper2D.cpp similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegEvaluationMapper2D.cpp rename to Modules/MatchPointRegistration/src/Rendering/mitkRegEvaluationMapper2D.cpp diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegVisColorStyleProperty.cpp b/Modules/MatchPointRegistration/src/Rendering/mitkRegVisColorStyleProperty.cpp similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegVisColorStyleProperty.cpp rename to Modules/MatchPointRegistration/src/Rendering/mitkRegVisColorStyleProperty.cpp diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegVisDirectionProperty.cpp b/Modules/MatchPointRegistration/src/Rendering/mitkRegVisDirectionProperty.cpp similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegVisDirectionProperty.cpp rename to Modules/MatchPointRegistration/src/Rendering/mitkRegVisDirectionProperty.cpp diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegVisHelper.cpp b/Modules/MatchPointRegistration/src/Rendering/mitkRegVisHelper.cpp similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegVisHelper.cpp rename to Modules/MatchPointRegistration/src/Rendering/mitkRegVisHelper.cpp diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegVisPropertyTags.cpp b/Modules/MatchPointRegistration/src/Rendering/mitkRegVisPropertyTags.cpp similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegVisPropertyTags.cpp rename to Modules/MatchPointRegistration/src/Rendering/mitkRegVisPropertyTags.cpp diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegVisStyleProperty.cpp b/Modules/MatchPointRegistration/src/Rendering/mitkRegVisStyleProperty.cpp similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegVisStyleProperty.cpp rename to Modules/MatchPointRegistration/src/Rendering/mitkRegVisStyleProperty.cpp diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegistrationWrapperMapper2D.cpp b/Modules/MatchPointRegistration/src/Rendering/mitkRegistrationWrapperMapper2D.cpp similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegistrationWrapperMapper2D.cpp rename to Modules/MatchPointRegistration/src/Rendering/mitkRegistrationWrapperMapper2D.cpp diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegistrationWrapperMapper3D.cpp b/Modules/MatchPointRegistration/src/Rendering/mitkRegistrationWrapperMapper3D.cpp similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegistrationWrapperMapper3D.cpp rename to Modules/MatchPointRegistration/src/Rendering/mitkRegistrationWrapperMapper3D.cpp diff --git a/Modules/MatchPointRegistration/Rendering/mitkRegistrationWrapperMapperBase.cpp b/Modules/MatchPointRegistration/src/Rendering/mitkRegistrationWrapperMapperBase.cpp similarity index 100% rename from Modules/MatchPointRegistration/Rendering/mitkRegistrationWrapperMapperBase.cpp rename to Modules/MatchPointRegistration/src/Rendering/mitkRegistrationWrapperMapperBase.cpp diff --git a/Modules/MatchPointRegistration/mitkMAPRegistrationWrapper.cpp b/Modules/MatchPointRegistration/src/mitkMAPRegistrationWrapper.cpp similarity index 100% rename from Modules/MatchPointRegistration/mitkMAPRegistrationWrapper.cpp rename to Modules/MatchPointRegistration/src/mitkMAPRegistrationWrapper.cpp diff --git a/Modules/MatchPointRegistration/mitkMAPRegistrationWrapperObjectFactory.cpp b/Modules/MatchPointRegistration/src/mitkMAPRegistrationWrapperObjectFactory.cpp similarity index 100% rename from Modules/MatchPointRegistration/mitkMAPRegistrationWrapperObjectFactory.cpp rename to Modules/MatchPointRegistration/src/mitkMAPRegistrationWrapperObjectFactory.cpp diff --git a/Modules/MatchPointRegistration/mitkRegEvaluationObject.cpp b/Modules/MatchPointRegistration/src/mitkRegEvaluationObject.cpp similarity index 100% rename from Modules/MatchPointRegistration/mitkRegEvaluationObject.cpp rename to Modules/MatchPointRegistration/src/mitkRegEvaluationObject.cpp diff --git a/Modules/MatchPointRegistration/mitkRegEvaluationObjectFactory.cpp b/Modules/MatchPointRegistration/src/mitkRegEvaluationObjectFactory.cpp similarity index 100% rename from Modules/MatchPointRegistration/mitkRegEvaluationObjectFactory.cpp rename to Modules/MatchPointRegistration/src/mitkRegEvaluationObjectFactory.cpp diff --git a/Modules/ModelFit/cmdapps/Fuse3Dto4DImageMiniApp.cpp b/Modules/ModelFit/cmdapps/Fuse3Dto4DImageMiniApp.cpp index a2cd40129f..de03cf4055 100644 --- a/Modules/ModelFit/cmdapps/Fuse3Dto4DImageMiniApp.cpp +++ b/Modules/ModelFit/cmdapps/Fuse3Dto4DImageMiniApp.cpp @@ -1,167 +1,167 @@ /*============================================================================ 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. ============================================================================*/ // std includes #include #include // itk includes #include "itksys/SystemTools.hxx" // CTK includes #include "mitkCommandLineParser.h" // MITK includes #include #include #include mitkCommandLineParser::StringContainerType inFilenames; std::string outFileName; std::vector images; std::vector timebounds; void setupParser(mitkCommandLineParser& parser) { // set general information about your MiniApp parser.setCategory("Dynamic Data Analysis Tools"); parser.setTitle("Fuse 3D to 4D Image"); parser.setDescription("MiniApp that allows to fuse several 3D images (with same geometry) into a 3D+t (4D) image that can be processed as dynamic data."); parser.setContributor("DKFZ MIC"); //! [create parser] //! [add arguments] // how should arguments be prefixed parser.setArgumentPrefix("--", "-"); // add each argument, unless specified otherwise each argument is optional // see mitkCommandLineParser::addArgument for more information parser.beginGroup("Required I/O parameters"); parser.addArgument( "inputs", "i", mitkCommandLineParser::StringList, "Input files", "Pathes to the input images that should be fused", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("output", "o", mitkCommandLineParser::File, "Output file path", "Path to the fused 3D+t image.", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.endGroup(); parser.beginGroup("Optional parameters"); parser.addArgument( "time", "t", mitkCommandLineParser::StringList, "Time bounds", "Defines the time geometry of the resulting dynamic image in [ms]. The first number is the start time point of the first time step. All other numbers are the max bound of a time step. So the structure is [minBound0 maxBound1 [maxBound2 [... maxBoundN]]]; e.g. \"2 3.5 10\" encodes a time geometry with two time steps and that starts at 2 ms and the second time step starts at 3.5 ms and ends at 10 ms. If not set e propertional time geometry with 1 ms duration will be generated!", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument("help", "h", mitkCommandLineParser::Bool, "Help:", "Show this help text"); parser.endGroup(); //! [add arguments] } bool configureApplicationSettings(std::map parsedArgs) { if (parsedArgs.size() == 0) return false; inFilenames = us::any_cast(parsedArgs["inputs"]); outFileName = us::any_cast(parsedArgs["output"]); if (parsedArgs.count("time")) { auto timeBoundsStr = us::any_cast(parsedArgs["time"]); for (const auto& timeBoundStr : timeBoundsStr) { std::istringstream stream; stream.imbue(std::locale("C")); stream.str(timeBoundStr); mitk::TimePointType time = 0 ; stream >> time; timebounds.emplace_back(time); } } return true; } int main(int argc, char* argv[]) { mitkCommandLineParser parser; setupParser(parser); - mitk::PreferenceListReaderOptionsFunctor readerFilterFunctor = mitk::PreferenceListReaderOptionsFunctor({ "MITK DICOM Reader v2 (classic config)" }, { "MITK DICOM Reader" }); + mitk::PreferenceListReaderOptionsFunctor readerFilterFunctor = mitk::PreferenceListReaderOptionsFunctor({ "MITK DICOM Reader v2 (autoselect)" }, { "" }); const std::map& parsedArgs = parser.parseArguments(argc, argv); if (!configureApplicationSettings(parsedArgs)) { return EXIT_FAILURE; }; // Show a help message if (parsedArgs.count("help") || parsedArgs.count("h")) { std::cout << parser.helpText(); return EXIT_SUCCESS; } if (timebounds.empty()) { timebounds.resize(inFilenames.size()+1); std::iota(timebounds.begin(), timebounds.end(), 0.0); } else if (inFilenames.size() + 1 != timebounds.size()) { std::cerr << "Cannot fuse images. Explicitly specified time bounds do not match. Use --help for more information on how to specify time bounds correctly."; return EXIT_FAILURE; }; //! [do processing] try { std::cout << "Load images:" << std::endl; auto filter = mitk::TemporalJoinImagesFilter::New(); unsigned int step = 0; for (auto path : inFilenames) { std::cout << "Time step #"<(path, &readerFilterFunctor); images.push_back(image); filter->SetInput(step, image); ++step; } filter->SetFirstMinTimeBound(timebounds[0]); filter->SetMaxTimeBounds({ timebounds.begin() + 1, timebounds.end() }); std::cout << "Fuse the images ..." << std::endl; filter->Update(); auto output = filter->GetOutput(); std::cout << "Save output image: " << outFileName << std::endl; mitk::IOUtil::Save(output, outFileName); std::cout << "Processing finished." << std::endl; return EXIT_SUCCESS; } catch (const std::exception& e) { MITK_ERROR << e.what(); return EXIT_FAILURE; } catch (...) { MITK_ERROR << "Unexpected error encountered."; return EXIT_FAILURE; } } diff --git a/Modules/ModelFit/cmdapps/GenericFittingMiniApp.cpp b/Modules/ModelFit/cmdapps/GenericFittingMiniApp.cpp index 68984e081d..0e45f362ed 100644 --- a/Modules/ModelFit/cmdapps/GenericFittingMiniApp.cpp +++ b/Modules/ModelFit/cmdapps/GenericFittingMiniApp.cpp @@ -1,360 +1,360 @@ /*============================================================================ 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. ============================================================================*/ // std includes #include // itk includes #include "itksys/SystemTools.hxx" // CTK includes #include "mitkCommandLineParser.h" // MITK includes #include #include #include #include #include #include #include #include #include #include #include #include std::string inFilename; std::string outFileName; std::string maskFileName; bool verbose(false); bool roibased(false); std::string functionName; std::string formular; mitk::Image::Pointer image; mitk::Image::Pointer mask; void onFitEvent(::itk::Object* caller, const itk::EventObject & event, void* /*data*/) { itk::ProgressEvent progressEvent; if (progressEvent.CheckEvent(&event)) { mitk::ParameterFitImageGeneratorBase* castedReporter = dynamic_cast(caller); std::cout <GetProgress()*100 << "% "; } } void setupParser(mitkCommandLineParser& parser) { // set general information about your MiniApp parser.setCategory("Dynamic Data Analysis Tools"); parser.setTitle("Generic Fitting"); parser.setDescription("MiniApp that allows to make a pixel based fitting on the intensity signal over time for a given model function."); parser.setContributor("DKFZ MIC"); //! [create parser] //! [add arguments] // how should arguments be prefixed parser.setArgumentPrefix("--", "-"); // add each argument, unless specified otherwise each argument is optional // see mitkCommandLineParser::addArgument for more information parser.beginGroup("Model parameters"); parser.addArgument( "function", "f", mitkCommandLineParser::String, "Model function", "Function that should be used to fit the intensity signals. Options are: \"Linear\" or \"\" (for generic formulas).", us::Any(std::string("Linear"))); parser.addArgument( "formular", "y", mitkCommandLineParser::String, "Generic model function formular", "Formular of a generic model (if selected) that will be parsed and fitted.", us::Any()); parser.endGroup(); parser.beginGroup("Required I/O parameters"); parser.addArgument( "input", "i", mitkCommandLineParser::File, "Input file", "input 3D+t image file", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("output", "o", mitkCommandLineParser::File, "Output file template", "where to save the output parameter images. The specified path will be used as template to determine the format (via extension) and the name \"root\". For each parameter a suffix will be added to the name.", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.endGroup(); parser.beginGroup("Optional parameters"); parser.addArgument( "mask", "m", mitkCommandLineParser::File, "Mask file", "Mask that defines the spatial image region that should be fitted. Must have the same geometry as the input image!", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument( "verbose", "v", mitkCommandLineParser::Bool, "Verbose Output", "Whether to produce verbose output"); parser.addArgument( "roibased", "r", mitkCommandLineParser::Bool, "Roi based fitting", "Will compute a mean intesity signal over the ROI before fitting it. If this mode is used a mask must be specified."); parser.addArgument("help", "h", mitkCommandLineParser::Bool, "Help:", "Show this help text"); parser.endGroup(); //! [add arguments] } bool configureApplicationSettings(std::map parsedArgs) { if (parsedArgs.size() == 0) return false; // parse, cast and set required arguments functionName = "Linear"; if (parsedArgs.count("function")) { functionName = us::any_cast(parsedArgs["function"]); } if (parsedArgs.count("formular")) { formular = us::any_cast(parsedArgs["formular"]); } inFilename = us::any_cast(parsedArgs["input"]); outFileName = us::any_cast(parsedArgs["output"]); verbose = false; if (parsedArgs.count("verbose")) { verbose = us::any_cast(parsedArgs["verbose"]); } roibased = false; if (parsedArgs.count("roibased")) { roibased = us::any_cast(parsedArgs["roibased"]); } if (parsedArgs.count("mask")) { maskFileName = us::any_cast(parsedArgs["mask"]); } return true; } void configureInitialParametersOfParameterizer(mitk::ModelParameterizerBase* parameterizer) { mitk::GenericParamModelParameterizer* genericParameterizer = dynamic_cast(parameterizer); if (genericParameterizer) { genericParameterizer->SetFunctionString(formular); } } mitk::ModelFitFunctorBase::Pointer createDefaultFitFunctor( const mitk::ModelParameterizerBase* parameterizer) { mitk::LevenbergMarquardtModelFitFunctor::Pointer fitFunctor = mitk::LevenbergMarquardtModelFitFunctor::New(); mitk::NormalizedSumOfSquaredDifferencesFitCostFunction::Pointer chi2 = mitk::NormalizedSumOfSquaredDifferencesFitCostFunction::New(); fitFunctor->RegisterEvaluationParameter("Chi^2", chi2); mitk::ModelBase::Pointer refModel = parameterizer->GenerateParameterizedModel(); ::itk::LevenbergMarquardtOptimizer::ScalesType scales; scales.SetSize(refModel->GetNumberOfParameters()); scales.Fill(1.0); fitFunctor->SetScales(scales); fitFunctor->SetDebugParameterMaps(true); return fitFunctor.GetPointer(); } template void generateModelFit_PixelBased(mitk::modelFit::ModelFitInfo::Pointer& /*modelFitInfo*/, mitk::ParameterFitImageGeneratorBase::Pointer& generator) { mitk::PixelBasedParameterFitImageGenerator::Pointer fitGenerator = mitk::PixelBasedParameterFitImageGenerator::New(); typename TParameterizer::Pointer modelParameterizer = TParameterizer::New(); configureInitialParametersOfParameterizer(modelParameterizer); //Specify fitting strategy and criterion parameters mitk::ModelFitFunctorBase::Pointer fitFunctor = createDefaultFitFunctor(modelParameterizer); //Parametrize fit generator fitGenerator->SetModelParameterizer(modelParameterizer); fitGenerator->SetMask(mask); fitGenerator->SetDynamicImage(image); fitGenerator->SetFitFunctor(fitFunctor); generator = fitGenerator.GetPointer(); } template void generateModelFit_ROIBased( mitk::modelFit::ModelFitInfo::Pointer& /*modelFitInfo*/, mitk::ParameterFitImageGeneratorBase::Pointer& generator) { mitk::ROIBasedParameterFitImageGenerator::Pointer fitGenerator = mitk::ROIBasedParameterFitImageGenerator::New(); typename TParameterizer::Pointer modelParameterizer = TParameterizer::New(); configureInitialParametersOfParameterizer(modelParameterizer); //Compute ROI signal mitk::MaskedDynamicImageStatisticsGenerator::Pointer signalGenerator = mitk::MaskedDynamicImageStatisticsGenerator::New(); signalGenerator->SetMask(mask); signalGenerator->SetDynamicImage(image); signalGenerator->Generate(); mitk::MaskedDynamicImageStatisticsGenerator::ResultType roiSignal = signalGenerator->GetMean(); //Specify fitting strategy and criterion parameters mitk::ModelFitFunctorBase::Pointer fitFunctor = createDefaultFitFunctor(modelParameterizer); //Parametrize fit generator fitGenerator->SetModelParameterizer(modelParameterizer); fitGenerator->SetMask(mask); fitGenerator->SetFitFunctor(fitFunctor); fitGenerator->SetSignal(roiSignal); fitGenerator->SetTimeGrid(mitk::ExtractTimeGrid(image)); generator = fitGenerator.GetPointer(); } void doFitting() { mitk::ParameterFitImageGeneratorBase::Pointer generator = nullptr; mitk::modelFit::ModelFitInfo::Pointer fitSession = nullptr; ::itk::CStyleCommand::Pointer command = ::itk::CStyleCommand::New(); command->SetCallback(onFitEvent); bool isLinearFactory = functionName == "Linear"; if (isLinearFactory) { std::cout << "Model: linear" << std::endl; if (!roibased) { generateModelFit_PixelBased(fitSession, generator); } else { generateModelFit_ROIBased(fitSession, generator); } } else { std::cout << "Model: generic (2 parameter)" << std::endl; if (!roibased) { generateModelFit_PixelBased(fitSession, generator); } else { generateModelFit_ROIBased(fitSession, generator); } } if (generator.IsNotNull() ) { std::cout << "Started fitting process..." << std::endl; generator->AddObserver(::itk::AnyEvent(), command); generator->Generate(); std::cout << std::endl << "Finished fitting process" << std::endl; mitk::storeModelFitGeneratorResults(outFileName, generator, fitSession); } else { mitkThrow() << "Fitting error! Could not initalize fitting job."; } } int main(int argc, char* argv[]) { mitkCommandLineParser parser; setupParser(parser); - mitk::PreferenceListReaderOptionsFunctor readerFilterFunctor = mitk::PreferenceListReaderOptionsFunctor({ "MITK DICOM Reader v2 (classic config)" }, { "MITK DICOM Reader" }); + mitk::PreferenceListReaderOptionsFunctor readerFilterFunctor = mitk::PreferenceListReaderOptionsFunctor({ "MITK DICOM Reader v2 (autoselect)" }, { "" }); const std::map& parsedArgs = parser.parseArguments(argc, argv); if (!configureApplicationSettings(parsedArgs)) { return EXIT_FAILURE; }; // Show a help message if (parsedArgs.count("help") || parsedArgs.count("h")) { std::cout << parser.helpText(); return EXIT_SUCCESS; } //! [do processing] try { image = mitk::IOUtil::Load(inFilename, &readerFilterFunctor); std::cout << "Input: " << inFilename << std::endl; if (!maskFileName.empty()) { mask = mitk::IOUtil::Load(maskFileName, &readerFilterFunctor); std::cout << "Mask: " << maskFileName << std::endl; } else { std::cout << "Mask: none" << std::endl; } if (roibased && mask.IsNull()) { mitkThrow() << "Error. Cannot fit. Please specify mask if you select roi based fitting."; } std::cout << "Style: "; if (roibased) { std::cout << "ROI based"; } else { std::cout << "pixel based"; } std::cout << std::endl; doFitting(); std::cout << "Processing finished." << std::endl; return EXIT_SUCCESS; } catch (const itk::ExceptionObject& e) { MITK_ERROR << e.what(); return EXIT_FAILURE; } catch (const std::exception& e) { MITK_ERROR << e.what(); return EXIT_FAILURE; } catch (...) { MITK_ERROR << "Unexpected error encountered."; return EXIT_FAILURE; } } diff --git a/Modules/ModelFit/cmdapps/PixelDumpMiniApp.cpp b/Modules/ModelFit/cmdapps/PixelDumpMiniApp.cpp index 2c50a02d4b..d3323ed319 100644 --- a/Modules/ModelFit/cmdapps/PixelDumpMiniApp.cpp +++ b/Modules/ModelFit/cmdapps/PixelDumpMiniApp.cpp @@ -1,447 +1,447 @@ /*============================================================================ 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. ============================================================================*/ // std includes #include // itk includes #include "itksys/SystemTools.hxx" #include "itkImageRegionConstIteratorWithIndex.h" #include "itkCastImageFilter.h" #include "itkExtractImageFilter.h" // CTK includes #include "mitkCommandLineParser.h" // MITK includes #include #include #include #include "mitkImageAccessByItk.h" #include "mitkImageCast.h" mitkCommandLineParser::StringContainerType inFilenames; std::string outFileName; std::string maskFileName; mitkCommandLineParser::StringContainerType captions; using ImageVectorType = std::vector; ImageVectorType images; mitk::Image::Pointer mask; bool verbose(false); typedef itk::Image InternalImageType; typedef std::map InternalImageMapType; InternalImageMapType internalImages; itk::ImageRegion<3> relevantRegion; InternalImageType::PointType relevantOrigin; InternalImageType::SpacingType relevantSpacing; InternalImageType::DirectionType relevantDirection; typedef itk::Index<3> DumpIndexType; typedef std::vector DumpedValuesType; struct DumpIndexCompare { bool operator() (const DumpIndexType& lhs, const DumpIndexType& rhs) const { if (lhs[0] < rhs[0]) { return true; } else if (lhs[0] > rhs[0]) { return false; } if (lhs[1] < rhs[1]) { return true; } else if (lhs[1] > rhs[1]) { return false; } return lhs[2] < rhs[2]; } }; typedef std::map DumpPixelMapType; DumpPixelMapType dumpedPixels; void setupParser(mitkCommandLineParser& parser) { // set general information about your MiniApp parser.setCategory("Generic Analysis Tools"); parser.setTitle("Pixel Dumper"); parser.setDescription("MiniApp that allows to dump the pixel values of all passed files into a csv. The region of dumping can defined by a mask. All images (and mask) must have the same geometrie."); parser.setContributor("DKFZ MIC"); //! [create parser] //! [add arguments] // how should arguments be prefixed parser.setArgumentPrefix("--", "-"); // add each argument, unless specified otherwise each argument is optional // see mitkCommandLineParser::addArgument for more information parser.beginGroup("Required I/O parameters"); parser.addArgument( "inputs", "i", mitkCommandLineParser::StringList, "Input files", "list of the images that should be dumped.", us::Any(), false); parser.addArgument("output", "o", mitkCommandLineParser::File, "Output file", "where to save the csv.", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.endGroup(); parser.beginGroup("Optional parameters"); parser.addArgument( "mask", "m", mitkCommandLineParser::File, "Mask file", "Mask that defines the spatial image region that should be dumped. Must have the same geometry as the input images!", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument( "captions", "c", mitkCommandLineParser::StringList, "Captions of image columns", "If provided the pixel columns of the csv will be named according to the passed values instead of using the image pathes. Number of images and names must be equal.", us::Any(), false); parser.addArgument("help", "h", mitkCommandLineParser::Bool, "Help:", "Show this help text"); parser.endGroup(); //! [add arguments] } bool configureApplicationSettings(std::map parsedArgs) { if (parsedArgs.size() == 0) return false; // parse, cast and set required arguments inFilenames = us::any_cast(parsedArgs["inputs"]); outFileName = us::any_cast(parsedArgs["output"]); if (parsedArgs.count("mask")) { maskFileName = us::any_cast(parsedArgs["mask"]); } captions = inFilenames; if (parsedArgs.count("captions")) { captions = us::any_cast(parsedArgs["captions"]); } return true; } template < typename TPixel, unsigned int VImageDimension > void ExtractRelevantInformation( const itk::Image< TPixel, VImageDimension > *image) { relevantRegion = image->GetLargestPossibleRegion(); relevantOrigin = image->GetOrigin(); relevantSpacing = image->GetSpacing(); relevantDirection = image->GetDirection(); } template < typename TPixel, unsigned int VImageDimension > void DoInternalImageConversion( const itk::Image< TPixel, VImageDimension > *image, InternalImageType::Pointer& internalImage) { typedef itk::Image< TPixel, VImageDimension > ImageType; //check if image fit to geometry // Make sure that spacing are the same typename ImageType::SpacingType imageSpacing = image->GetSpacing(); typename ImageType::PointType zeroPoint; zeroPoint.Fill(0.0); if ((zeroPoint + imageSpacing).SquaredEuclideanDistanceTo((zeroPoint + relevantSpacing)) > 1e-6) // for the dumper we are not as strict as mitk normally would be (mitk::eps) { mitkThrow() << "Images need to have same spacing! (Image spacing: " << imageSpacing << "; relevant spacing: " << relevantSpacing << ")"; } // Make sure that orientation of mask and image are the same typename ImageType::DirectionType imageDirection = image->GetDirection(); for (unsigned int i = 0; i < imageDirection.RowDimensions; ++i) { for (unsigned int j = 0; j < imageDirection.ColumnDimensions; ++j) { double differenceDirection = imageDirection[i][j] - relevantDirection[i][j]; if (fabs(differenceDirection) > 1e-6) // SD: 1e6 wird hier zum zweiten mal als Magic Number benutzt -> Konstante { // for the dumper we are not as strict as mitk normally would be (mitk::eps) mitkThrow() << "Images need to have same direction! (Image direction: " << imageDirection << "; relevant direction: " << relevantDirection << ")"; } } } // Make sure that origin of mask and image are the same typename ImageType::PointType imageOrigin = image->GetOrigin(); if (imageOrigin.SquaredEuclideanDistanceTo(relevantOrigin) > 1e-6) { // for the dumper we are not as strict as mitk normally would be (mitk::eps) mitkThrow() << "Image need to have same spacing! (Image spacing: " << imageSpacing << "; relevant spacing: " << relevantOrigin << ")"; } typename ImageType::RegionType imageRegion = image->GetLargestPossibleRegion(); if (!imageRegion.IsInside(relevantRegion) && imageRegion != relevantRegion) { mitkThrow() << "Images need to have same region! (Image region: " << imageRegion << "; relevant region: " << relevantRegion << ")"; } //convert to internal image typedef itk::ExtractImageFilter ExtractFilterType; typename ExtractFilterType::Pointer extractFilter = ExtractFilterType::New(); typedef itk::CastImageFilter CastFilterType; typename CastFilterType::Pointer castFilter = CastFilterType::New(); extractFilter->SetInput(image); extractFilter->SetExtractionRegion(relevantRegion); castFilter->SetInput(extractFilter->GetOutput()); castFilter->Update(); internalImage = castFilter->GetOutput(); } template < typename TPixel, unsigned int VImageDimension > void DoMaskedCollecting( const itk::Image< TPixel, VImageDimension > *image) { typedef itk::Image< TPixel, VImageDimension > ImageType; itk::ImageRegionConstIteratorWithIndex it(image, relevantRegion); it.GoToBegin(); while (!it.IsAtEnd()) { if (mask.IsNull() || it.Get() > 0) { DumpedValuesType values; const auto index = it.GetIndex(); for (auto& imagePos : internalImages) { double value = imagePos.second->GetPixel(index); values.push_back(value); } dumpedPixels.insert(std::make_pair(index, values)); } ++it; } } InternalImageMapType ConvertImageTimeSteps(mitk::Image* image) { InternalImageMapType map; InternalImageType::Pointer internalImage; for (unsigned int i = 0; i < image->GetTimeSteps(); ++i) { mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); imageTimeSelector->SetInput(image); imageTimeSelector->SetTimeNr(i); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image::Pointer imageTimePoint = imageTimeSelector->GetOutput(); AccessFixedDimensionByItk_1(imageTimePoint, DoInternalImageConversion, 3, internalImage); std::stringstream stream; stream << "[" << i << "]"; map.insert(std::make_pair(stream.str(), internalImage)); } return map; } void doDumping() { if (mask.IsNotNull() && mask->GetTimeSteps() > 1) { std::cout << "Pixel Dumper: Selected mask has multiple timesteps. Only use first timestep to mask the pixel dumping." << std::endl; mitk::ImageTimeSelector::Pointer maskTimeSelector = mitk::ImageTimeSelector::New(); maskTimeSelector->SetInput(mask); maskTimeSelector->SetTimeNr(0); maskTimeSelector->UpdateLargestPossibleRegion(); mask = maskTimeSelector->GetOutput(); } try { if (mask.IsNotNull()) { // if mask exist, we use the mask because it could be a sub region. AccessFixedDimensionByItk(mask, ExtractRelevantInformation, 3); } else { AccessFixedDimensionByItk(images.front(), ExtractRelevantInformation, 3); } } catch (const std::exception& e) { std::cerr << "Error extracting image geometry. Error text: " << e.what(); throw; } for (unsigned int index = 0; index < images.size(); ++index) { try { InternalImageMapType conversionMap = ConvertImageTimeSteps(images[index]); if (conversionMap.size() == 1) { internalImages.insert(std::make_pair(captions[index], conversionMap.begin()->second)); } else if (conversionMap.size() > 1) { for (auto& pos : conversionMap) { std::stringstream namestream; namestream << captions[index] << " " << pos.first; internalImages.insert(std::make_pair(namestream.str(), pos.second)); } } } catch (const std::exception& e) { std::stringstream errorStr; errorStr << "Inconsistent image \"" << captions[index] << "\" will be excluded from the collection. Error: " << e.what(); std::cerr << errorStr.str() << std::endl; } } if (mask.IsNotNull()) { // if mask exist, we use the mask because it could be a sub region. AccessFixedDimensionByItk(mask, DoMaskedCollecting, 3); } else { AccessFixedDimensionByItk(images.front(), DoMaskedCollecting, 3); } } void storeCSV() { std::ofstream resultfile; resultfile.open(outFileName.c_str()); resultfile << "x,y,z"; for (auto aImage : internalImages) { resultfile << "," << aImage.first; } resultfile << std::endl; for (auto dumpPos : dumpedPixels) { resultfile << dumpPos.first[0] << "," << dumpPos.first[1] << "," << dumpPos.first[2]; for(auto d : dumpPos.second) { resultfile << "," << d; } resultfile << std::endl; } } int main(int argc, char* argv[]) { mitkCommandLineParser parser; setupParser(parser); const std::map& parsedArgs = parser.parseArguments(argc, argv); - mitk::PreferenceListReaderOptionsFunctor readerFilterFunctor = mitk::PreferenceListReaderOptionsFunctor({ "MITK DICOM Reader v2 (classic config)" }, { "MITK DICOM Reader" }); + mitk::PreferenceListReaderOptionsFunctor readerFilterFunctor = mitk::PreferenceListReaderOptionsFunctor({ "MITK DICOM Reader v2 (autoselect)" }, { "" }); if (!configureApplicationSettings(parsedArgs)) { return EXIT_FAILURE; }; // Show a help message if (parsedArgs.count("help") || parsedArgs.count("h")) { std::cout << parser.helpText(); return EXIT_SUCCESS; } if (!captions.empty() && inFilenames.size() != captions.size()) { std::cerr << "Cannot dump images. Number of given captions does not equal number of given images."; return EXIT_FAILURE; }; //! [do processing] try { std::cout << "Load images:" << std::endl; for (auto path : inFilenames) { std::cout << "Input: " << path << std::endl; auto image = mitk::IOUtil::Load(path, &readerFilterFunctor); images.push_back(image); } if (!maskFileName.empty()) { mask = mitk::IOUtil::Load(maskFileName, &readerFilterFunctor); std::cout << "Mask: " << maskFileName << std::endl; } else { std::cout << "Mask: none" << std::endl; } doDumping(); storeCSV(); std::cout << "Processing finished." << std::endl; return EXIT_SUCCESS; } catch (const itk::ExceptionObject& e) { MITK_ERROR << e.what(); return EXIT_FAILURE; } catch (const std::exception& e) { MITK_ERROR << e.what(); return EXIT_FAILURE; } catch (...) { MITK_ERROR << "Unexpected error encountered."; return EXIT_FAILURE; } } diff --git a/Modules/Pharmacokinetics/cmdapps/CurveDescriptorMiniApp.cpp b/Modules/Pharmacokinetics/cmdapps/CurveDescriptorMiniApp.cpp index dd52eec7f8..df8257534f 100644 --- a/Modules/Pharmacokinetics/cmdapps/CurveDescriptorMiniApp.cpp +++ b/Modules/Pharmacokinetics/cmdapps/CurveDescriptorMiniApp.cpp @@ -1,242 +1,242 @@ /*============================================================================ 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. ============================================================================*/ // std includes #include // itk includes #include "itksys/SystemTools.hxx" // CTK includes #include "mitkCommandLineParser.h" // MITK includes #include #include #include #include #include #include #include #include #include std::string inFilename; std::string outFileName; std::string maskFileName; bool verbose(false); bool preview(false); mitk::Image::Pointer image; mitk::Image::Pointer mask; void onFitEvent(::itk::Object* caller, const itk::EventObject & event, void* /*data*/) { itk::ProgressEvent progressEvent; if (progressEvent.CheckEvent(&event)) { mitk::PixelBasedDescriptionParameterImageGenerator* castedReporter = dynamic_cast(caller); std::cout <GetProgress()*100 << "% "; } } void setupParser(mitkCommandLineParser& parser) { // set general information about your MiniApp parser.setCategory("Dynamic Data Analysis Tools"); parser.setTitle("Curve Descriptor"); parser.setDescription("MiniApp that allows to generate curve descriptor maps for dynamic image."); parser.setContributor("DKFZ MIC"); //! [create parser] //! [add arguments] // how should arguments be prefixed parser.setArgumentPrefix("--", "-"); // add each argument, unless specified otherwise each argument is optional // see mitkCommandLineParser::addArgument for more information parser.beginGroup("Required I/O parameters"); parser.addArgument( "input", "i", mitkCommandLineParser::File, "Input file", "input 3D+t image file", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("output", "o", mitkCommandLineParser::File, "Output file template", "where to save the output parameter images. The specified path will be used as template to determine the format (via extension) and the name \"root\". For each parameter a suffix will be added to the name.", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.endGroup(); parser.beginGroup("Optional parameters"); parser.addArgument( "mask", "m", mitkCommandLineParser::File, "Mask file", "Mask that defines the spatial image region that should be fitted. Must have the same geometry as the input image!", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument( "verbose", "v", mitkCommandLineParser::Bool, "Verbose Output", "Whether to produce verbose output"); parser.addArgument( "preview", "p", mitkCommandLineParser::Bool, "Preview outputs", "The application previews the outputs (filename, type) it would produce with the current settings."); parser.addArgument("help", "h", mitkCommandLineParser::Bool, "Help:", "Show this help text"); parser.endGroup(); //! [add arguments] } bool configureApplicationSettings(std::map parsedArgs) { if (parsedArgs.size() == 0) return false; verbose = false; if (parsedArgs.count("verbose")) { verbose = us::any_cast(parsedArgs["verbose"]); } if (parsedArgs.count("mask")) { maskFileName = us::any_cast(parsedArgs["mask"]); } preview = false; if (parsedArgs.count("preview")) { preview = us::any_cast(parsedArgs["preview"]); } inFilename = us::any_cast(parsedArgs["input"]); outFileName = us::any_cast(parsedArgs["output"]); return true; } void ConfigureFunctor(mitk::CurveParameterFunctor* functor) { mitk::CurveDescriptionParameterBase::Pointer parameterFunction = mitk::AreaUnderTheCurveDescriptionParameter::New().GetPointer(); functor->RegisterDescriptionParameter("AUC", parameterFunction); parameterFunction = mitk::AreaUnderFirstMomentDescriptionParameter::New().GetPointer(); functor->RegisterDescriptionParameter("AUMC", parameterFunction); parameterFunction = mitk::MeanResidenceTimeDescriptionParameter::New().GetPointer(); functor->RegisterDescriptionParameter("MRT", parameterFunction); parameterFunction = mitk::TimeToPeakCurveDescriptionParameter::New().GetPointer(); functor->RegisterDescriptionParameter("TimeToPeak", parameterFunction); }; void doDescription() { mitk::PixelBasedDescriptionParameterImageGenerator::Pointer generator = mitk::PixelBasedDescriptionParameterImageGenerator::New(); mitk::CurveParameterFunctor::Pointer functor = mitk::CurveParameterFunctor::New(); ConfigureFunctor(functor); functor->SetGrid(mitk::ExtractTimeGrid(image)); generator->SetFunctor(functor); generator->SetDynamicImage(image); generator->SetMask(mask); ::itk::CStyleCommand::Pointer command = ::itk::CStyleCommand::New(); command->SetCallback(onFitEvent); std::cout << "Started curve descriptor computation process..." << std::endl; generator->AddObserver(::itk::AnyEvent(), command); generator->Generate(); std::cout << std::endl << "Finished computation process" << std::endl; for (auto imageIterator : generator->GetParameterImages()) { mitk::storeParameterResultImage(outFileName, imageIterator.first, imageIterator.second); } } void doPreview() { mitk::CurveParameterFunctor::Pointer functor = mitk::CurveParameterFunctor::New(); ConfigureFunctor(functor); auto pNames = functor->GetDescriptionParameterNames(); for (auto aName : pNames) { auto fullPath = mitk::generateModelFitResultImagePath(outFileName, aName); std::cout << "Store result parameter: " << aName << " -> " << fullPath << std::endl; } } int main(int argc, char* argv[]) { mitkCommandLineParser parser; setupParser(parser); - mitk::PreferenceListReaderOptionsFunctor readerFilterFunctor = mitk::PreferenceListReaderOptionsFunctor({ "MITK DICOM Reader v2 (classic config)" }, { "MITK DICOM Reader" }); + mitk::PreferenceListReaderOptionsFunctor readerFilterFunctor = mitk::PreferenceListReaderOptionsFunctor({ "MITK DICOM Reader v2 (autoselect)" }, { "" }); const std::map& parsedArgs = parser.parseArguments(argc, argv); if (!configureApplicationSettings(parsedArgs)) { return EXIT_FAILURE; }; // Show a help message if (parsedArgs.count("help") || parsedArgs.count("h")) { std::cout << parser.helpText(); return EXIT_SUCCESS; } //! [do processing] try { if (preview) { doPreview(); } else { image = mitk::IOUtil::Load(inFilename, &readerFilterFunctor); std::cout << "Input: " << inFilename << std::endl; if (!maskFileName.empty()) { mask = mitk::IOUtil::Load(maskFileName, &readerFilterFunctor); std::cout << "Mask: " << maskFileName << std::endl; } else { std::cout << "Mask: none" << std::endl; } doDescription(); } std::cout << "Processing finished." << std::endl; return EXIT_SUCCESS; } catch (const itk::ExceptionObject& e) { MITK_ERROR << e.what(); return EXIT_FAILURE; } catch (const std::exception& e) { MITK_ERROR << e.what(); return EXIT_FAILURE; } catch (...) { MITK_ERROR << "Unexpected error encountered."; return EXIT_FAILURE; } } diff --git a/Modules/Pharmacokinetics/cmdapps/MRPerfusionMiniApp.cpp b/Modules/Pharmacokinetics/cmdapps/MRPerfusionMiniApp.cpp index 9fc08615b9..d871108746 100644 --- a/Modules/Pharmacokinetics/cmdapps/MRPerfusionMiniApp.cpp +++ b/Modules/Pharmacokinetics/cmdapps/MRPerfusionMiniApp.cpp @@ -1,890 +1,890 @@ /*============================================================================ 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. ============================================================================*/ // std includes #include // itk includes #include "itksys/SystemTools.hxx" // CTK includes #include "mitkCommandLineParser.h" // MITK includes #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include std::string inFilename; std::string outFileName; std::string maskFileName; std::string aifMaskFileName; std::string aifImageFileName; mitk::Image::Pointer image; mitk::Image::Pointer mask; mitk::Image::Pointer aifImage; mitk::Image::Pointer aifMask; bool useConstraints(false); bool verbose(false); bool roibased(false); bool preview(false); std::string modelName; float aifHematocritLevel(0); float brixInjectionTime(0); const std::string MODEL_NAME_2SL = "2SL"; const std::string MODEL_NAME_3SL = "3SL"; const std::string MODEL_NAME_descriptive = "descriptive"; const std::string MODEL_NAME_tofts = "tofts"; const std::string MODEL_NAME_2CX = "2CX"; void onFitEvent(::itk::Object* caller, const itk::EventObject & event, void* /*data*/) { itk::ProgressEvent progressEvent; if (progressEvent.CheckEvent(&event)) { mitk::ParameterFitImageGeneratorBase* castedReporter = dynamic_cast(caller); std::cout <GetProgress()*100 << "% "; } } void setupParser(mitkCommandLineParser& parser) { // set general information about your MiniApp parser.setCategory("Dynamic Data Analysis Tools"); parser.setTitle("MR Perfusion"); parser.setDescription("MiniApp that allows to fit MRI perfusion models and generates the according parameter maps. IMPORTANT!!!: The app assumes that the input images (signal and AIF) are concentration images. If your images do not hold this assumption, convert the image date before using this app (e.g. by using the signal-to-concentration-converter mini app."); parser.setContributor("DKFZ MIC"); //! [create parser] //! [add arguments] // how should arguments be prefixed parser.setArgumentPrefix("--", "-"); // add each argument, unless specified otherwise each argument is optional // see mitkCommandLineParser::addArgument for more information parser.beginGroup("Model parameters"); parser.addArgument( "model", "l", mitkCommandLineParser::String, "Model function", "Model that should be used to fit the concentration signal. Options are: \""+MODEL_NAME_descriptive+"\" (descriptive pharmacokinetic Brix model),\"" + MODEL_NAME_2SL + "\" (two step linear model),\""+MODEL_NAME_3SL+"\" (three step linear model), \""+MODEL_NAME_tofts+"\" (extended tofts model) or \""+MODEL_NAME_2CX+"\" (two compartment exchange model).", us::Any(std::string(MODEL_NAME_tofts))); parser.addArgument( "injectiontime", "j", mitkCommandLineParser::Float, "Injection time [min]", "Injection time of the bolus. This information is needed for the descriptive pharmacokinetic Brix model.", us::Any()); parser.endGroup(); parser.beginGroup("Required I/O parameters"); parser.addArgument( "input", "i", mitkCommandLineParser::File, "Input file", "input 3D+t image file", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("output", "o", mitkCommandLineParser::File, "Output file template", "where to save the output parameter images. The specified path will be used as template to determine the format (via extension) and the name \"root\". For each parameter a suffix will be added to the name.", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.endGroup(); parser.beginGroup("AIF parameters"); parser.addArgument( "aifmask", "n", mitkCommandLineParser::File, "AIF mask file", "Mask that defines the spatial image region that should be used as AIF for models that need one. Must have the same geometry as the AIF input image!", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument( "aifimage", "a", mitkCommandLineParser::File, "AIF image file", "3D+t image that defines the image that containes the AIF signal. If this flag is not set and the model needs a AIF, the CLI will assume that the AIF is encoded in the normal image. Must have the same geometry as the AIF mask!", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument( "hematocrit", "h", mitkCommandLineParser::Float, "Hematocrit Level", "Value needed for correct AIF computation. Only needed if model needs an AIF. Default value is 0.45.", us::Any(0.45)); parser.endGroup(); parser.beginGroup("Optional parameters"); parser.addArgument( "mask", "m", mitkCommandLineParser::File, "Mask file", "Mask that defines the spatial image region that should be fitted. Must have the same geometry as the input image!", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument( "verbose", "v", mitkCommandLineParser::Bool, "Verbose Output", "Whether to produce verbose output"); parser.addArgument( "roibased", "r", mitkCommandLineParser::Bool, "Roi based fitting", "Will compute a mean intesity signal over the ROI before fitting it. If this mode is used a mask must be specified."); parser.addArgument( "constraints", "c", mitkCommandLineParser::Bool, "Constraints", "Indicates if constraints should be used for the fitting (if flag is set the default contraints will be used.).", us::Any(false)); parser.addArgument( "preview", "p", mitkCommandLineParser::Bool, "Preview outputs", "The application previews the outputs (filename, type) it would produce with the current settings."); parser.addArgument("help", "h", mitkCommandLineParser::Bool, "Help:", "Show this help text"); parser.endGroup(); //! [add arguments] } bool configureApplicationSettings(std::map parsedArgs) { if (parsedArgs.size() == 0) return false; // parse, cast and set required arguments modelName = MODEL_NAME_tofts; if (parsedArgs.count("model")) { modelName = us::any_cast(parsedArgs["model"]); } inFilename = us::any_cast(parsedArgs["input"]); outFileName = us::any_cast(parsedArgs["output"]); if (parsedArgs.count("mask")) { maskFileName = us::any_cast(parsedArgs["mask"]); } if (parsedArgs.count("aifimage")) { aifImageFileName = us::any_cast(parsedArgs["aifimage"]); } if (parsedArgs.count("aifmask")) { aifMaskFileName = us::any_cast(parsedArgs["aifmask"]); } verbose = false; if (parsedArgs.count("verbose")) { verbose = us::any_cast(parsedArgs["verbose"]); } preview = false; if (parsedArgs.count("preview")) { preview = us::any_cast(parsedArgs["preview"]); } roibased = false; if (parsedArgs.count("roibased")) { roibased = us::any_cast(parsedArgs["roibased"]); } useConstraints = false; if (parsedArgs.count("constraints")) { useConstraints = us::any_cast(parsedArgs["constraints"]); } aifHematocritLevel = 0.45; if (parsedArgs.count("hematocrit")) { aifHematocritLevel = us::any_cast(parsedArgs["hematocrit"]); } brixInjectionTime = 0.0; if (parsedArgs.count("injectiontime")) { brixInjectionTime = us::any_cast(parsedArgs["injectiontime"]); } return true; } mitk::ModelFitFunctorBase::Pointer createDefaultFitFunctor( const mitk::ModelParameterizerBase* parameterizer, const mitk::ModelFactoryBase* modelFactory) { mitk::LevenbergMarquardtModelFitFunctor::Pointer fitFunctor = mitk::LevenbergMarquardtModelFitFunctor::New(); mitk::NormalizedSumOfSquaredDifferencesFitCostFunction::Pointer chi2 = mitk::NormalizedSumOfSquaredDifferencesFitCostFunction::New(); fitFunctor->RegisterEvaluationParameter("Chi^2", chi2); if (useConstraints) { fitFunctor->SetConstraintChecker(modelFactory->CreateDefaultConstraints().GetPointer()); } mitk::ModelBase::Pointer refModel = parameterizer->GenerateParameterizedModel(); ::itk::LevenbergMarquardtOptimizer::ScalesType scales; scales.SetSize(refModel->GetNumberOfParameters()); scales.Fill(1.0); fitFunctor->SetScales(scales); fitFunctor->SetDebugParameterMaps(true); return fitFunctor.GetPointer(); } /**Helper that ensures that the mask (if it exists) is always 3D image. If the mask is originally an 4D image, the first time step will be used.*/ mitk::Image::Pointer getMask3D() { mitk::Image::Pointer result; if (mask.IsNotNull()) { result = mask; //mask settings if (mask->GetTimeSteps() > 1) { MITK_INFO << "Selected mask has multiple timesteps. Only use first timestep to mask model fit."; mitk::ImageTimeSelector::Pointer maskedImageTimeSelector = mitk::ImageTimeSelector::New(); maskedImageTimeSelector->SetInput(mask); maskedImageTimeSelector->SetTimeNr(0); maskedImageTimeSelector->UpdateLargestPossibleRegion(); result = maskedImageTimeSelector->GetOutput(); } } return result; } void getAIF(mitk::AIFBasedModelBase::AterialInputFunctionType& aif, mitk::AIFBasedModelBase::AterialInputFunctionType& aifTimeGrid) { if (aifMask.IsNotNull()) { aif.clear(); aifTimeGrid.clear(); mitk::AterialInputFunctionGenerator::Pointer aifGenerator = mitk::AterialInputFunctionGenerator::New(); //Hematocrit level aifGenerator->SetHCL(aifHematocritLevel); std::cout << "AIF hematocrit level: " << aifHematocritLevel << std::endl; mitk::Image::Pointer selectedAIFMask = aifMask; //mask settings if (aifMask->GetTimeSteps() > 1) { MITK_INFO << "Selected AIF mask has multiple timesteps. Only use first timestep to mask model fit."; mitk::ImageTimeSelector::Pointer maskedImageTimeSelector = mitk::ImageTimeSelector::New(); maskedImageTimeSelector->SetInput(aifMask); maskedImageTimeSelector->SetTimeNr(0); maskedImageTimeSelector->UpdateLargestPossibleRegion(); aifMask = maskedImageTimeSelector->GetOutput(); } aifGenerator->SetMask(aifMask); mitk::Image::Pointer selectedAIFImage = image; //image settings if (aifImage.IsNotNull()) { selectedAIFImage = aifImage; } aifGenerator->SetDynamicImage(selectedAIFImage); aif = aifGenerator->GetAterialInputFunction(); aifTimeGrid = aifGenerator->GetAterialInputFunctionTimeGrid(); } else { mitkThrow() << "Cannot generate AIF. AIF mask was not specified or correctly loaded."; } } void generateDescriptiveBrixModel_PixelBased(mitk::modelFit::ModelFitInfo::Pointer& modelFitInfo, mitk::ParameterFitImageGeneratorBase::Pointer& generator) { mitk::PixelBasedParameterFitImageGenerator::Pointer fitGenerator = mitk::PixelBasedParameterFitImageGenerator::New(); mitk::DescriptivePharmacokineticBrixModelParameterizer::Pointer modelParameterizer = mitk::DescriptivePharmacokineticBrixModelParameterizer::New(); mitk::Image::Pointer mask3D = getMask3D(); //Model configuration (static parameters) can be done now modelParameterizer->SetTau(brixInjectionTime); std::cout << "Injection time [min]: " << brixInjectionTime << std::endl; mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); imageTimeSelector->SetInput(image); imageTimeSelector->SetTimeNr(0); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::DescriptivePharmacokineticBrixModelParameterizer::BaseImageType::Pointer baseImage; mitk::CastToItkImage(imageTimeSelector->GetOutput(), baseImage); modelParameterizer->SetBaseImage(baseImage); //Specify fitting strategy and criterion parameters mitk::ModelFactoryBase::Pointer factory = mitk::DescriptivePharmacokineticBrixModelFactory::New().GetPointer(); mitk::ModelFitFunctorBase::Pointer fitFunctor = createDefaultFitFunctor(modelParameterizer, factory); //Parametrize fit generator fitGenerator->SetModelParameterizer(modelParameterizer); std::string roiUID = ""; if (mask3D.IsNotNull()) { fitGenerator->SetMask(mask3D); roiUID = mask->GetUID(); } fitGenerator->SetDynamicImage(image); fitGenerator->SetFitFunctor(fitFunctor); generator = fitGenerator.GetPointer(); //Create model info modelFitInfo = mitk::modelFit::CreateFitInfoFromModelParameterizer(modelParameterizer, image, mitk::ModelFitConstants::FIT_TYPE_VALUE_PIXELBASED(), roiUID); } void generateDescriptiveBrixModel_ROIBased(mitk::modelFit::ModelFitInfo::Pointer& modelFitInfo, mitk::ParameterFitImageGeneratorBase::Pointer& generator) { mitk::Image::Pointer mask3D = getMask3D(); if (mask3D.IsNull()) { return; } mitk::ROIBasedParameterFitImageGenerator::Pointer fitGenerator = mitk::ROIBasedParameterFitImageGenerator::New(); mitk::DescriptivePharmacokineticBrixModelValueBasedParameterizer::Pointer modelParameterizer = mitk::DescriptivePharmacokineticBrixModelValueBasedParameterizer::New(); //Compute ROI signal mitk::MaskedDynamicImageStatisticsGenerator::Pointer signalGenerator = mitk::MaskedDynamicImageStatisticsGenerator::New(); signalGenerator->SetMask(mask3D); signalGenerator->SetDynamicImage(image); signalGenerator->Generate(); mitk::MaskedDynamicImageStatisticsGenerator::ResultType roiSignal = signalGenerator->GetMean(); //Model configuration (static parameters) can be done now modelParameterizer->SetTau(brixInjectionTime); std::cout << "Injection time [min]: " << brixInjectionTime << std::endl; modelParameterizer->SetBaseValue(roiSignal[0]); //Specify fitting strategy and criterion parameters mitk::ModelFactoryBase::Pointer factory = mitk::DescriptivePharmacokineticBrixModelFactory::New().GetPointer(); mitk::ModelFitFunctorBase::Pointer fitFunctor = createDefaultFitFunctor(modelParameterizer, factory); //Parametrize fit generator fitGenerator->SetModelParameterizer(modelParameterizer); fitGenerator->SetMask(mask3D); fitGenerator->SetFitFunctor(fitFunctor); fitGenerator->SetSignal(roiSignal); fitGenerator->SetTimeGrid(mitk::ExtractTimeGrid(image)); generator = fitGenerator.GetPointer(); std::string roiUID = mask->GetUID(); //Create model info modelFitInfo = mitk::modelFit::CreateFitInfoFromModelParameterizer(modelParameterizer, image, mitk::ModelFitConstants::FIT_TYPE_VALUE_ROIBASED(), roiUID); mitk::ScalarListLookupTable::ValueType infoSignal; for (mitk::MaskedDynamicImageStatisticsGenerator::ResultType::const_iterator pos = roiSignal.begin(); pos != roiSignal.end(); ++pos) { infoSignal.push_back(*pos); } modelFitInfo->inputData.SetTableValue("ROI", infoSignal); } template void GenerateLinearModelFit_PixelBased(mitk::modelFit::ModelFitInfo::Pointer& modelFitInfo, mitk::ParameterFitImageGeneratorBase::Pointer& generator) { mitk::PixelBasedParameterFitImageGenerator::Pointer fitGenerator = mitk::PixelBasedParameterFitImageGenerator::New(); typename TParameterizer::Pointer modelParameterizer = TParameterizer::New(); mitk::Image::Pointer mask3D = getMask3D(); //Specify fitting strategy and criterion parameters mitk::ModelFactoryBase::Pointer factory = TFactory::New().GetPointer(); mitk::ModelFitFunctorBase::Pointer fitFunctor = createDefaultFitFunctor(modelParameterizer, factory); //Parametrize fit generator fitGenerator->SetModelParameterizer(modelParameterizer); std::string roiUID = ""; if (mask3D.IsNotNull()) { fitGenerator->SetMask(mask3D); roiUID = mask->GetUID(); } fitGenerator->SetDynamicImage(image); fitGenerator->SetFitFunctor(fitFunctor); generator = fitGenerator.GetPointer(); //Create model info modelFitInfo = mitk::modelFit::CreateFitInfoFromModelParameterizer(modelParameterizer, image, mitk::ModelFitConstants::FIT_TYPE_VALUE_PIXELBASED(), roiUID); } template void GenerateLinearModelFit_ROIBased(mitk::modelFit::ModelFitInfo::Pointer& modelFitInfo, mitk::ParameterFitImageGeneratorBase::Pointer& generator) { mitk::Image::Pointer mask3D = getMask3D(); if (mask3D.IsNull()) { return; } mitk::ROIBasedParameterFitImageGenerator::Pointer fitGenerator = mitk::ROIBasedParameterFitImageGenerator::New(); typename TParameterizer::Pointer modelParameterizer = TParameterizer::New(); //Compute ROI signal mitk::MaskedDynamicImageStatisticsGenerator::Pointer signalGenerator = mitk::MaskedDynamicImageStatisticsGenerator::New(); signalGenerator->SetMask(mask3D); signalGenerator->SetDynamicImage(image); signalGenerator->Generate(); mitk::MaskedDynamicImageStatisticsGenerator::ResultType roiSignal = signalGenerator->GetMean(); //Specify fitting strategy and criterion parameters mitk::ModelFactoryBase::Pointer factory = TFactory::New().GetPointer(); mitk::ModelFitFunctorBase::Pointer fitFunctor = createDefaultFitFunctor(modelParameterizer, factory); //Parametrize fit generator fitGenerator->SetModelParameterizer(modelParameterizer); fitGenerator->SetMask(mask3D); fitGenerator->SetFitFunctor(fitFunctor); fitGenerator->SetSignal(roiSignal); fitGenerator->SetTimeGrid(mitk::ExtractTimeGrid(image)); generator = fitGenerator.GetPointer(); std::string roiUID = mask->GetUID(); //Create model info modelFitInfo = mitk::modelFit::CreateFitInfoFromModelParameterizer(modelParameterizer, image, mitk::ModelFitConstants::FIT_TYPE_VALUE_ROIBASED(), roiUID); mitk::ScalarListLookupTable::ValueType infoSignal; for (mitk::MaskedDynamicImageStatisticsGenerator::ResultType::const_iterator pos = roiSignal.begin(); pos != roiSignal.end(); ++pos) { infoSignal.push_back(*pos); } modelFitInfo->inputData.SetTableValue("ROI", infoSignal); } template void generateAIFbasedModelFit_PixelBased(mitk::modelFit::ModelFitInfo::Pointer& modelFitInfo, mitk::ParameterFitImageGeneratorBase::Pointer& generator) { mitk::PixelBasedParameterFitImageGenerator::Pointer fitGenerator = mitk::PixelBasedParameterFitImageGenerator::New(); typename TParameterizer::Pointer modelParameterizer = TParameterizer::New(); mitk::AIFBasedModelBase::AterialInputFunctionType aif; mitk::AIFBasedModelBase::AterialInputFunctionType aifTimeGrid; getAIF(aif, aifTimeGrid); modelParameterizer->SetAIF(aif); modelParameterizer->SetAIFTimeGrid(aifTimeGrid); mitk::Image::Pointer mask3D = getMask3D(); //Specify fitting strategy and criterion parameters mitk::ModelFactoryBase::Pointer factory = TFactory::New().GetPointer(); mitk::ModelFitFunctorBase::Pointer fitFunctor = createDefaultFitFunctor(modelParameterizer, factory); //Parametrize fit generator fitGenerator->SetModelParameterizer(modelParameterizer); std::string roiUID = ""; if (mask3D.IsNotNull()) { fitGenerator->SetMask(mask3D); roiUID = mask->GetUID(); } fitGenerator->SetDynamicImage(image); fitGenerator->SetFitFunctor(fitFunctor); generator = fitGenerator.GetPointer(); //Create model info modelFitInfo = mitk::modelFit::CreateFitInfoFromModelParameterizer(modelParameterizer, image, mitk::ModelFitConstants::FIT_TYPE_VALUE_PIXELBASED(), roiUID); mitk::ScalarListLookupTable::ValueType infoSignal; for (mitk::AIFBasedModelBase::AterialInputFunctionType::const_iterator pos = aif.begin(); pos != aif.end(); ++pos) { infoSignal.push_back(*pos); } modelFitInfo->inputData.SetTableValue("AIF", infoSignal); } template void generateAIFbasedModelFit_ROIBased( mitk::modelFit::ModelFitInfo::Pointer& modelFitInfo, mitk::ParameterFitImageGeneratorBase::Pointer& generator) { mitk::Image::Pointer mask3D = getMask3D(); if (mask3D.IsNull()) { return; } mitk::ROIBasedParameterFitImageGenerator::Pointer fitGenerator = mitk::ROIBasedParameterFitImageGenerator::New(); typename TParameterizer::Pointer modelParameterizer = TParameterizer::New(); mitk::AIFBasedModelBase::AterialInputFunctionType aif; mitk::AIFBasedModelBase::AterialInputFunctionType aifTimeGrid; getAIF(aif, aifTimeGrid); modelParameterizer->SetAIF(aif); modelParameterizer->SetAIFTimeGrid(aifTimeGrid); //Compute ROI signal mitk::MaskedDynamicImageStatisticsGenerator::Pointer signalGenerator = mitk::MaskedDynamicImageStatisticsGenerator::New(); signalGenerator->SetMask(mask3D); signalGenerator->SetDynamicImage(image); signalGenerator->Generate(); mitk::MaskedDynamicImageStatisticsGenerator::ResultType roiSignal = signalGenerator->GetMean(); //Specify fitting strategy and criterion parameters mitk::ModelFactoryBase::Pointer factory = TFactory::New().GetPointer(); mitk::ModelFitFunctorBase::Pointer fitFunctor = createDefaultFitFunctor(modelParameterizer, factory); //Parametrize fit generator fitGenerator->SetModelParameterizer(modelParameterizer); fitGenerator->SetMask(mask3D); fitGenerator->SetFitFunctor(fitFunctor); fitGenerator->SetSignal(roiSignal); fitGenerator->SetTimeGrid(mitk::ExtractTimeGrid(image)); generator = fitGenerator.GetPointer(); std::string roiUID = mask->GetUID(); //Create model info modelFitInfo = mitk::modelFit::CreateFitInfoFromModelParameterizer(modelParameterizer, image, mitk::ModelFitConstants::FIT_TYPE_VALUE_ROIBASED(), roiUID); mitk::ScalarListLookupTable::ValueType infoSignal; for (mitk::MaskedDynamicImageStatisticsGenerator::ResultType::const_iterator pos = roiSignal.begin(); pos != roiSignal.end(); ++pos) { infoSignal.push_back(*pos); } modelFitInfo->inputData.SetTableValue("ROI", infoSignal); infoSignal.clear(); for (mitk::AIFBasedModelBase::AterialInputFunctionType::const_iterator pos = aif.begin(); pos != aif.end(); ++pos) { infoSignal.push_back(*pos); } modelFitInfo->inputData.SetTableValue("AIF", infoSignal); } void storeResultImage(const std::string& name, mitk::Image* image, mitk::modelFit::Parameter::Type nodeType, const mitk::modelFit::ModelFitInfo* modelFitInfo) { mitk::modelFit::SetModelFitDataProperties(image, name, nodeType, modelFitInfo); std::string ext = ::itksys::SystemTools::GetFilenameLastExtension(outFileName); std::string dir = itksys::SystemTools::GetFilenamePath(outFileName); dir = itksys::SystemTools::ConvertToOutputPath(dir); std::string rootName = itksys::SystemTools::GetFilenameWithoutLastExtension(outFileName); std::string fileName = rootName + "_" + name + ext; std::vector pathElements; pathElements.push_back(dir); pathElements.push_back(fileName); std::string fullOutPath = itksys::SystemTools::ConvertToOutputPath(dir + "/" + fileName); mitk::IOUtil::Save(image, fullOutPath); std::cout << "Store result (parameter: "<(fitSession, generator); } else { GenerateLinearModelFit_ROIBased(fitSession, generator); } } else if (is2SLFactory) { std::cout << "Model: two step linear model" << std::endl; if (!roibased) { GenerateLinearModelFit_PixelBased(fitSession, generator); } else { GenerateLinearModelFit_ROIBased(fitSession, generator); } } else if (isToftsFactory) { std::cout << "Model: extended tofts model" << std::endl; if (!roibased) { generateAIFbasedModelFit_PixelBased(fitSession, generator); } else { generateAIFbasedModelFit_ROIBased(fitSession, generator); } } else if (is2CXMFactory) { std::cout << "Model: two compartment exchange model" << std::endl; if (!roibased) { generateAIFbasedModelFit_PixelBased(fitSession, generator); } else { generateAIFbasedModelFit_ROIBased(fitSession, generator); } } else { std::cerr << "ERROR. Model flag is unknown. Given flag: " << modelName << std::endl; } } void doFitting() { mitk::ParameterFitImageGeneratorBase::Pointer generator = nullptr; mitk::modelFit::ModelFitInfo::Pointer fitSession = nullptr; ::itk::CStyleCommand::Pointer command = ::itk::CStyleCommand::New(); command->SetCallback(onFitEvent); createFitGenerator(fitSession, generator); if (generator.IsNotNull() ) { std::cout << "Started fitting process..." << std::endl; generator->AddObserver(::itk::AnyEvent(), command); generator->Generate(); std::cout << std::endl << "Finished fitting process" << std::endl; mitk::storeModelFitGeneratorResults(outFileName, generator, fitSession); } else { mitkThrow() << "Fitting error! Could not initialize fitting job."; } } void doPreview() { mitk::ParameterFitImageGeneratorBase::Pointer generator = nullptr; mitk::modelFit::ModelFitInfo::Pointer fitSession = nullptr; createFitGenerator(fitSession, generator); if (generator.IsNotNull()) { mitk::previewModelFitGeneratorResults(outFileName, generator); } else { mitkThrow() << "Fitting error! Could not initialize fitting job."; } } int main(int argc, char* argv[]) { mitkCommandLineParser parser; setupParser(parser); const std::map& parsedArgs = parser.parseArguments(argc, argv); - mitk::PreferenceListReaderOptionsFunctor readerFilterFunctor = mitk::PreferenceListReaderOptionsFunctor({ "MITK DICOM Reader v2 (classic config)" }, { "MITK DICOM Reader" }); + mitk::PreferenceListReaderOptionsFunctor readerFilterFunctor = mitk::PreferenceListReaderOptionsFunctor({ "MITK DICOM Reader v2 (autoselect)" }, { "" }); if (!configureApplicationSettings(parsedArgs)) { return EXIT_FAILURE; }; // Show a help message if (parsedArgs.count("help") || parsedArgs.count("h")) { std::cout << parser.helpText(); return EXIT_SUCCESS; } //! [do processing] try { image = mitk::IOUtil::Load(inFilename, &readerFilterFunctor); std::cout << "Input: " << inFilename << std::endl; if (!maskFileName.empty()) { mask = mitk::IOUtil::Load(maskFileName, &readerFilterFunctor); std::cout << "Mask: " << maskFileName << std::endl; } else { std::cout << "Mask: none" << std::endl; } if (modelName != MODEL_NAME_descriptive && modelName != MODEL_NAME_3SL && MODEL_NAME_2SL != modelName) { if (!aifMaskFileName.empty()) { aifMask = mitk::IOUtil::Load(aifMaskFileName, &readerFilterFunctor); std::cout << "AIF mask: " << aifMaskFileName << std::endl; } else { mitkThrow() << "Error. Cannot fit. Choosen model needs an AIF. Please specify AIF mask (--aifmask)."; } if (!aifImageFileName.empty()) { aifImage = mitk::IOUtil::Load(aifImageFileName, &readerFilterFunctor); std::cout << "AIF image: " << aifImageFileName << std::endl; } else { std::cout << "AIF image: none (using signal image)" << std::endl; } } if (roibased && mask.IsNull()) { mitkThrow() << "Error. Cannot fit. Please specify mask if you select roi based fitting."; } std::cout << "Style: "; if (roibased) { std::cout << "ROI based"; } else { std::cout << "pixel based"; } std::cout << std::endl; if (preview) { doPreview(); } else { doFitting(); } std::cout << "Processing finished." << std::endl; return EXIT_SUCCESS; } catch (const itk::ExceptionObject& e) { MITK_ERROR << e.what(); return EXIT_FAILURE; } catch (const std::exception& e) { MITK_ERROR << e.what(); return EXIT_FAILURE; } catch (...) { MITK_ERROR << "Unexpected error encountered."; return EXIT_FAILURE; } } diff --git a/Modules/Pharmacokinetics/cmdapps/MRSignal2ConcentrationMiniApp.cpp b/Modules/Pharmacokinetics/cmdapps/MRSignal2ConcentrationMiniApp.cpp index bd517b3abb..c9057a5bc8 100644 --- a/Modules/Pharmacokinetics/cmdapps/MRSignal2ConcentrationMiniApp.cpp +++ b/Modules/Pharmacokinetics/cmdapps/MRSignal2ConcentrationMiniApp.cpp @@ -1,287 +1,287 @@ /*============================================================================ 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. ============================================================================*/ // std includes #include // itk includes #include "itksys/SystemTools.hxx" // CTK includes #include "mitkCommandLineParser.h" // MITK includes #include #include #include #include #include std::string inFilename; std::string outFileName; mitk::Image::Pointer image; bool verbose(false); bool t1_absolute(false); bool t1_relative(false); bool t1_flash(false); bool t2(false); float k(1.0); float te(0); float rec_time(0); float relaxivity(0); float rel_time(0); void setupParser(mitkCommandLineParser& parser) { // set general information about your MiniApp parser.setCategory("Dynamic Data Analysis Tools"); parser.setTitle("MR Signal to Concentration Converter"); parser.setDescription("MiniApp that allows to convert a T1 or T2 signal image into a concentration image for perfusion analysis."); parser.setContributor("DKFZ MIC"); //! [create parser] //! [add arguments] // how should arguments be prefixed parser.setArgumentPrefix("--", "-"); // add each argument, unless specified otherwise each argument is optional // see mitkCommandLineParser::addArgument for more information parser.beginGroup("Required I/O parameters"); parser.addArgument( "input", "i", mitkCommandLineParser::File, "Input file", "input 3D+t image file", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("output", "o", mitkCommandLineParser::File, "Output file", "where to save the output concentration image.", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.endGroup(); parser.beginGroup("Conversion parameters"); parser.addArgument( "t1-absolute", "", mitkCommandLineParser::Bool, "T1 absolute signal enhancement", "Activate conversion for T1 absolute signal enhancement."); parser.addArgument( "t1-relative", "", mitkCommandLineParser::Bool, "T1 relative signal enhancement", "Activate conversion for T1 relative signal enhancement."); parser.addArgument( "t1-flash", "", mitkCommandLineParser::Bool, "T1 turbo flash", "Activate specific conversion for T1 turbo flash sequences."); parser.addArgument( "t2", "", mitkCommandLineParser::Bool, "T2 signal conversion", "Activate conversion for T2 signal enhancement to concentration."); parser.addArgument( "k", "k", mitkCommandLineParser::Float, "Conversion factor k", "Needed for the following conversion modes: T1-absolute, T1-relative, T2. Default value is 1.", us::Any(1)); parser.addArgument( "recovery-time", "", mitkCommandLineParser::Float, "Recovery time", "Needed for the following conversion modes: T1-flash."); parser.addArgument( "relaxivity", "", mitkCommandLineParser::Float, "Relaxivity", "Needed for the following conversion modes: T1-flash."); parser.addArgument( "relaxation-time", "", mitkCommandLineParser::Float, "Relaxation time", "Needed for the following conversion modes: T1-flash."); parser.addArgument( "te", "", mitkCommandLineParser::Float, "Echo time TE", "Needed for the following conversion modes: T2.", us::Any(1)); parser.beginGroup("Optional parameters"); parser.addArgument( "verbose", "v", mitkCommandLineParser::Bool, "Verbose Output", "Whether to produce verbose output"); parser.addArgument("help", "h", mitkCommandLineParser::Bool, "Help:", "Show this help text"); parser.endGroup(); //! [add arguments] } bool configureApplicationSettings(std::map parsedArgs) { if (parsedArgs.size() == 0) return false; inFilename = us::any_cast(parsedArgs["input"]); outFileName = us::any_cast(parsedArgs["output"]); verbose = false; if (parsedArgs.count("verbose")) { verbose = us::any_cast(parsedArgs["verbose"]); } t1_absolute = false; if (parsedArgs.count("t1-absolute")) { t1_absolute = us::any_cast(parsedArgs["t1-absolute"]); } t1_relative = false; if (parsedArgs.count("t1-relative")) { t1_relative = us::any_cast(parsedArgs["t1-relative"]); } t1_flash = false; if (parsedArgs.count("t1-flash")) { t1_flash = us::any_cast(parsedArgs["t1-flash"]); } t2 = false; if (parsedArgs.count("t2")) { t2 = us::any_cast(parsedArgs["t2"]); } k = 0.0; if (parsedArgs.count("k")) { k = us::any_cast(parsedArgs["k"]); } relaxivity = 0.0; if (parsedArgs.count("relaxivity")) { relaxivity = us::any_cast(parsedArgs["relaxivity"]); } rec_time = 0.0; if (parsedArgs.count("recovery-time")) { rec_time = us::any_cast(parsedArgs["recovery-time"]); } rel_time = 0.0; if (parsedArgs.count("relaxation-time")) { rel_time = us::any_cast(parsedArgs["relaxation-time"]); } te = 0.0; if (parsedArgs.count("te")) { te = us::any_cast(parsedArgs["te"]); } //consistency checks int modeCount = 0; if (t1_absolute) ++modeCount; if (t1_flash) ++modeCount; if (t1_relative) ++modeCount; if (t2) ++modeCount; if (modeCount==0) { mitkThrow() << "Invalid program call. Please select the type of conversion."; } if (modeCount >1) { mitkThrow() << "Invalid program call. Please select only ONE type of conversion."; } if (!k && (t2 || t1_absolute || t1_relative)) { mitkThrow() << "Invalid program call. Please set 'k', if you use t1-absolute, t1-relative or t2."; } if (!te && t2) { mitkThrow() << "Invalid program call. Please set 'te', if you use t2 mode."; } if ((!rec_time||!rel_time||!relaxivity) && t1_flash) { mitkThrow() << "Invalid program call. Please set 'recovery-time', 'relaxation-time' and 'relaxivity', if you use t1-flash mode."; } return true; } void doConversion() { mitk::ConcentrationCurveGenerator::Pointer concentrationGen = mitk::ConcentrationCurveGenerator::New(); concentrationGen->SetDynamicImage(image); concentrationGen->SetisTurboFlashSequence(t1_flash); concentrationGen->SetAbsoluteSignalEnhancement(t1_absolute); concentrationGen->SetRelativeSignalEnhancement(t1_relative); concentrationGen->SetisT2weightedImage(t2); if (t1_flash) { concentrationGen->SetRecoveryTime(rec_time); concentrationGen->SetRelaxationTime(rel_time); concentrationGen->SetRelaxivity(relaxivity); } else if (t2) { concentrationGen->SetT2Factor(k); concentrationGen->SetT2EchoTime(te); } else { concentrationGen->SetFactor(k); } mitk::Image::Pointer concentrationImage = concentrationGen->GetConvertedImage(); mitk::IOUtil::Save(concentrationImage, outFileName); std::cout << "Store result: " << outFileName << std::endl; } int main(int argc, char* argv[]) { mitkCommandLineParser parser; setupParser(parser); const std::map& parsedArgs = parser.parseArguments(argc, argv); if (!configureApplicationSettings(parsedArgs)) { return EXIT_FAILURE; }; - mitk::PreferenceListReaderOptionsFunctor readerFilterFunctor = mitk::PreferenceListReaderOptionsFunctor({ "MITK DICOM Reader v2 (classic config)" }, { "MITK DICOM Reader" }); + mitk::PreferenceListReaderOptionsFunctor readerFilterFunctor = mitk::PreferenceListReaderOptionsFunctor({ "MITK DICOM Reader v2 (autoselect)" }, { "" }); // Show a help message if (parsedArgs.count("help") || parsedArgs.count("h")) { std::cout << parser.helpText(); return EXIT_SUCCESS; } //! [do processing] try { image = mitk::IOUtil::Load(inFilename, &readerFilterFunctor); std::cout << "Input: " << inFilename << std::endl; doConversion(); std::cout << "Processing finished." << std::endl; return EXIT_SUCCESS; } catch (const itk::ExceptionObject& e) { MITK_ERROR << e.what(); return EXIT_FAILURE; } catch (const std::exception& e) { MITK_ERROR << e.what(); return EXIT_FAILURE; } catch (...) { MITK_ERROR << "Unexpected error encountered."; return EXIT_FAILURE; } }