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 &regionSize = 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