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..20efdf8766
--- /dev/null
+++ b/Modules/MatchPointRegistration/Testing/itkStitchImageFilterTest.cpp
@@ -0,0 +1,73 @@
+/*============================================================================
+
+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"
+
+class itkStitchImageFilterTestSuite : public mitk::TestFixture
+{
+ CPPUNIT_TEST_SUITE(itkStitchImageFilterTestSuite);
+ // Test the append method
+ MITK_TEST(StitchWithNoTransformAndNoInterp);
+ CPPUNIT_TEST_SUITE_END();
+
+ using InputImageType = mitk::TestImageType;
+ using OutputImageType = itk::Image<double, 2>;
+private:
+ using FilterType = itk::StitchImageFilter<InputImageType, OutputImageType>;
+ 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.);
+ m_Filter = FilterType::New();
+ m_Input1 = mitk::GenerateTestImage(1);
+
+ m_Input2 = mitk::GenerateTestImage(2);
+ origin[1] = 3.;
+ m_Input2->SetOrigin(origin);
+
+ m_Input3 = mitk::GenerateTestImage(3);
+ origin[1] = 6.;
+ m_Input3->SetOrigin(origin);
+
+ FilterType::SizeType size = { 3, 9 };
+ m_Filter->SetSize(size);
+ }
+
+ void tearDown() override
+ {
+ }
+
+ void StitchWithNoTransformAndNoInterp()
+ {
+ m_Filter->SetInput(0, m_Input1);
+ m_Filter->SetInput(1, m_Input2);
+ m_Filter->SetInput(2, m_Input3);
+
+ m_Filter->Update();
+ auto output = m_Filter->GetOutput();
+ }
+
+
+};
+
+MITK_TEST_SUITE_REGISTRATION(itkStitchImageFilter)
diff --git a/Modules/MatchPointRegistration/include/itkStitchImageFilter.h b/Modules/MatchPointRegistration/include/itkStitchImageFilter.h
new file mode 100644
index 0000000000..e9f4c0bcc1
--- /dev/null
+++ b/Modules/MatchPointRegistration/include/itkStitchImageFilter.h
@@ -0,0 +1,321 @@
+/*============================================================================
+
+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
+{
+/** \class StitchImageFilter
+ * \brief ITK filter that resamples/stitches multiple images into a given reference geometry.
+ *
+ * StitchImageFilter is simelar 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 mutliple input images cover the output region a wighted sum of all
+ * mapped input pixel values will be calculated.
+ * It resamples images image through some coordinate
+ * transform, interpolating via some image function. The class is templated
+ * over the types of the input and output images.
+ *
+ * All other behaviors are simelar 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<TransformType> 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<PixelType> 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<ImageDimension> ReferenceImageBaseType;
+
+ virtual void SetInput(unsigned int index, const InputImageType* image);
+ /** 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);
+
+ /** 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;
+
+ /** Default implementation for resampling that works for any
+ * transformation type. */
+ virtual void NonlinearThreadedGenerateData(const OutputImageRegionType &
+ outputRegionForThread,
+ ThreadIdType threadId);
+
+ /** Implementation for resampling that works for with linear
+ * transformation types.
+ */
+ virtual void LinearThreadedGenerateData(const OutputImageRegionType &
+ outputRegionForThread,
+ ThreadIdType threadId);
+
+ /** 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<const InputImageType*> InputImageVectorType;
+ typedef std::map<const InputImageType*, typename TransformType::ConstPointer> TransformMapType;
+ typedef std::map<const InputImageType*, InterpolatorPointerType> 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;
+
+};
+} // 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..d37687f089
--- /dev/null
+++ b/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp
@@ -0,0 +1,754 @@
+/*============================================================================
+
+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 <numeric>
+
+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_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<PixelType>::ZeroValue( m_DefaultPixelValue );
+}
+
+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<TransformType*>(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();
+
+ unsigned int nComponents
+ = DefaultConvertPixelTraits<PixelType>::GetNumberOfComponents(
+ m_DefaultPixelValue );
+
+ if (nComponents == 0)
+ {
+ PixelComponentType zeroComponent
+ = NumericTraits<PixelComponentType>::ZeroValue( zeroComponent );
+ nComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
+ NumericTraits<PixelType>::SetLength(m_DefaultPixelValue, nComponents );
+ for (unsigned int n=0; n<nComponents; n++)
+ {
+ PixelConvertType::SetNthComponent( n, m_DefaultPixelValue,
+ zeroComponent );
+ }
+ }
+}
+
+template< typename TInputImage,
+ typename TOutputImage,
+ typename TInterpolatorPrecisionType,
+ typename TTransformPrecisionType >
+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)
+{
+ // Check whether the input or the output is a
+ // SpecialCoordinatesImage. If either are, then we cannot use the
+ // fast path since index mapping will definitely not be linear.
+ typedef SpecialCoordinatesImage< PixelType, ImageDimension > OutputSpecialCoordinatesImageType;
+ typedef SpecialCoordinatesImage< InputPixelType, InputImageDimension > InputSpecialCoordinatesImageType;
+
+ if( outputRegionForThread.GetNumberOfPixels() == 0 )
+ {
+ return;
+ }
+
+ this->NonlinearThreadedGenerateData(outputRegionForThread, threadId);
+}
+
+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<PixelType>::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<PixelComponentType>( minComponent ) );
+ }
+ else if ( component > maxComponent )
+ {
+ PixelConvertType::SetNthComponent( n, outputValue, static_cast<PixelComponentType>( maxComponent ) );
+ }
+ else
+ {
+ PixelConvertType::SetNthComponent(n, outputValue,
+ static_cast<PixelComponentType>( component ) );
+ }
+ }
+
+ return outputValue;
+}
+
+template< typename TInputImage,
+ typename TOutputImage,
+ typename TInterpolatorPrecisionType,
+ typename TTransformPrecisionType >
+void
+StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType >
+::NonlinearThreadedGenerateData(const OutputImageRegionType &
+ outputRegionForThread,
+ ThreadIdType threadId)
+{
+ // Get the output pointers
+ OutputImageType *outputPtr = this->GetOutput();
+
+ // Get this input pointers
+ InputImageVectorType inputs = this->GetInputs();
+ TransformMapType transforms = this->GetTransforms();
+
+ // Create an iterator that will walk the output region for this thread.
+ typedef ImageRegionIteratorWithIndex< TOutputImage > 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< ComponentType >( minValue );
+ const ComponentType maxOutputValue = static_cast< ComponentType >( 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<PixelType> pixvals;
+
+ 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));
+ }
+ }
+
+ if (!pixvals.empty())
+ { //at least one input provided a value -> compute weighted sum
+ double sum = std::accumulate(pixvals.begin(), pixvals.end(), 0.0);
+ outIt.Set(sum / pixvals.size());
+ }
+ else
+ {
+ outIt.Set(m_DefaultPixelValue); // default background value
+ }
+
+ progress.CompletedPixel();
+ ++outIt;
+ }
+}
+
+template< typename TInputImage,
+ typename TOutputImage,
+ typename TInterpolatorPrecisionType,
+ typename TTransformPrecisionType >
+void
+StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType >
+::LinearThreadedGenerateData(const OutputImageRegionType &
+ outputRegionForThread,
+ ThreadIdType threadId)
+{
+ // Get the output pointers
+ OutputImageType *outputPtr = this->GetOutput();
+
+ // Get this input pointers
+ InputImageVectorType inputs = this->GetInputs();
+ TransformMapType transforms = this->GetTransforms();
+
+ // Create an iterator that will walk the output region for this thread.
+ typedef ImageScanlineIterator< TOutputImage > 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 tmpOutputPoint;
+
+ std::map<const InputImageType*, ContinuousInputIndexType> inputIndices;
+
+ typedef typename PointType::VectorType VectorType;
+ std::map<const InputImageType*, VectorType> deltas; // delta in input continuous index coordinate frame
+
+ IndexType index;
+
+ const typename OutputImageRegionType::SizeType ®ionSize = outputRegionForThread.GetSize();
+ const SizeValueType numberOfLinesToProcess = outputRegionForThread.GetNumberOfPixels() / regionSize[0];
+
+ // Support for progress methods/callbacks
+ ProgressReporter progress( this,
+ threadId,
+ numberOfLinesToProcess );
+
+ typedef typename InterpolatorType::OutputType OutputType;
+
+ // Cache information from the superclass
+ PixelType defaultValue = this->GetDefaultPixelValue();
+
+ // 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< ComponentType >( minValue );
+ const ComponentType maxOutputValue = static_cast< ComponentType >( maxValue );
+
+ // Determine the position of the first pixel in the scanline
+ index = outIt.GetIndex();
+ outputPtr->TransformIndexToPhysicalPoint(index, outputPoint);
+
+ // Compute corresponding input pixel position
+ for (const auto& input : inputs)
+ {
+ auto inputPoint = transforms[input]->TransformPoint(outputPoint);
+
+ ContinuousInputIndexType inputIndex;
+ input->TransformPhysicalPointToContinuousIndex(inputPoint, inputIndex);
+ inputIndices[input] = inputIndex;
+ }
+
+ // As we walk across a scan line in the output image, we trace
+ // an oriented/scaled/translated line in the input image. Cache
+ // the delta along this line in continuous index space of the input
+ // image. This allows us to use vector addition to model the
+ // transformation.
+ //
+ // To find delta, we take two pixels adjacent in a scanline
+ // and determine the continuous indices of these pixels when
+ // mapped to the input coordinate frame. We use the difference
+ // between these two continuous indices as the delta to apply
+ // to an index to trace line in the input image as we move
+ // across the scanline of the output image.
+ //
+ // We determine delta in this manner so that Images
+ // are both handled properly (taking into account the direction cosines).
+ //
+ ++index[0];
+ outputPtr->TransformIndexToPhysicalPoint(index, tmpOutputPoint);
+ for (const auto& input : inputs)
+ {
+ auto tmpInputPoint = transforms[input]->TransformPoint(tmpOutputPoint);
+
+ ContinuousInputIndexType tmpInputIndex;
+ input->TransformPhysicalPointToContinuousIndex(tmpInputPoint,
+ tmpInputIndex);
+ deltas[input] = tmpInputIndex - inputIndices[input];
+ }
+
+ while ( !outIt.IsAtEnd() )
+ {
+ // Determine the continuous index of the first pixel of output
+ // scanline when mapped to the input coordinate frame.
+ //
+
+ // First get the position of the pixel in the output coordinate frame
+ index = outIt.GetIndex();
+ outputPtr->TransformIndexToPhysicalPoint(index, outputPoint);
+
+ // Compute corresponding input pixel continuous index, this index
+ // will incremented in the scanline loop
+ for (const auto& input : inputs)
+ {
+ auto inputPoint = transforms[input]->TransformPoint(outputPoint);
+ input->TransformPhysicalPointToContinuousIndex(inputPoint, inputIndices[input]);
+ }
+
+ while ( !outIt.IsAtEndOfLine() )
+ {
+ std::vector<PixelType> pixvals;
+
+ for (const auto& input : inputs)
+ {
+ if (m_Interpolators[input]->IsInsideBuffer(inputIndices[input]))
+ {
+ auto value = m_Interpolators[input]->EvaluateAtContinuousIndex(inputIndices[input]);
+ pixvals.emplace_back(this->CastPixelWithBoundsChecking(value, minOutputValue, maxOutputValue));
+ }
+ inputIndices[input] += deltas[input];
+ }
+
+ // Evaluate input at right position and copy to the output
+ if ( !pixvals.empty())
+ {
+ double sum = std::accumulate(pixvals.begin(), pixvals.end(), 0.0);
+ outIt.Set(sum / pixvals.size());
+ }
+ else
+ {
+ outIt.Set(defaultValue); // default background value
+ }
+
+ ++outIt;
+ }
+ progress.CompletedPixel();
+ outIt.NextLine();
+ }
+}
+
+template<typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType>
+typename StitchImageFilter<TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType>::InputImageVectorType
+StitchImageFilter<TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType>
+::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 TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType>
+typename StitchImageFilter<TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType>::TransformMapType
+StitchImageFilter<TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType>
+::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<TInputImage*>(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();
+ 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())
+ {
+ m_Interpolators[input] = LinearInterpolatorType::New().GetPointer();
+ }
+ }
+}
+
+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