diff --git a/Modules/MatchPointRegistration/Testing/itkStitchImageFilterTest.cpp b/Modules/MatchPointRegistration/Testing/itkStitchImageFilterTest.cpp
index 20efdf8766..916e0c0cad 100644
--- a/Modules/MatchPointRegistration/Testing/itkStitchImageFilterTest.cpp
+++ b/Modules/MatchPointRegistration/Testing/itkStitchImageFilterTest.cpp
@@ -1,73 +1,152 @@
/*============================================================================
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<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 = mitk::GenerateTestImage(10);
+ origin[1] = 2.;
m_Input2->SetOrigin(origin);
- m_Input3 = mitk::GenerateTestImage(3);
- origin[1] = 6.;
+ m_Input3 = mitk::GenerateTestImage(100);
+ origin[1] = 4.;
m_Input3->SetOrigin(origin);
FilterType::SizeType size = { 3, 9 };
+ m_Filter->SetDefaultPixelValue(1000);
m_Filter->SetSize(size);
}
void tearDown() override
{
}
+ bool CheckPixels(const OutputImageType* image, const std::vector<double>& pixels)
+ {
+ bool result = true;
+
+ itk::ImageRegionConstIteratorWithIndex<OutputImageType> 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."<<std::endl;
+ return false;
+ }
+
+ if (*refIter != iter.Get())
+ {
+ std::cerr << "Checked image differs from reference. Index: " << iter.GetIndex() << "; value: " << iter.Get() << "; ref: " << *refIter <<std::endl;
+ result = false;
+ }
+ ++iter;
+ ++refIter;
+ }
+
+ return result;
+ }
+
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();
+
+ 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<double, 2>;
+ TranslationType::OutputVectorType offset;
+ offset[0] = 0;
+ offset[1] = -1;
+ auto translation1 = TranslationType::New();
+ translation1->SetOffset(offset);
+ m_Filter->SetInput(1, m_Input2, translation1);
+
+ offset[1] = -2;
+ 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<double, 2>;
+ TranslationType::OutputVectorType offset;
+ offset[0] = 0;
+ offset[1] = -1.5;
+ auto translation1 = TranslationType::New();
+ translation1->SetOffset(offset);
+
+ offset[1] = -2.5;
+ auto translation2 = TranslationType::New();
+ translation2->SetOffset(offset);
+
+ m_Filter->SetInput(0, m_Input1);
+ m_Filter->SetInput(1, m_Input2, translation1, itk::NearestNeighborInterpolateImageFunction<InputImageType>::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/include/itkStitchImageFilter.h b/Modules/MatchPointRegistration/include/itkStitchImageFilter.h
index e9f4c0bcc1..b42d255e1c 100644
--- a/Modules/MatchPointRegistration/include/itkStitchImageFilter.h
+++ b/Modules/MatchPointRegistration/include/itkStitchImageFilter.h
@@ -1,321 +1,308 @@
/*============================================================================
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
index d37687f089..fa1ffaa4b4 100644
--- a/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp
+++ b/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp
@@ -1,754 +1,600 @@
/*============================================================================
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();
+ for (const auto& interpolator : m_Interpolators)
+ {
+ interpolator.second->SetInputImage(interpolator.first);
+ }
+
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();
+ 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() );
+ 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();
+ 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 );
+ const ComponentType minOutputValue = static_cast<ComponentType>(minValue);
+ const ComponentType maxOutputValue = static_cast<ComponentType>(maxValue);
// Walk the output region
outIt.GoToBegin();
- while ( !outIt.IsAtEnd() )
+ 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
+typename StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType >
+::PixelType
StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType >
-::LinearThreadedGenerateData(const OutputImageRegionType &
- outputRegionForThread,
- ThreadIdType threadId)
+::CastPixelWithBoundsChecking(const InterpolatorOutputType value,
+ const ComponentType minComponent,
+ const ComponentType maxComponent ) const
{
- // 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.
- //
+ const unsigned int nComponents = InterpolatorConvertType::GetNumberOfComponents(value);
+ PixelType outputValue;
- // First get the position of the pixel in the output coordinate frame
- index = outIt.GetIndex();
- outputPtr->TransformIndexToPhysicalPoint(index, outputPoint);
+ NumericTraits<PixelType>::SetLength( outputValue, nComponents );
- // Compute corresponding input pixel continuous index, this index
- // will incremented in the scanline loop
- for (const auto& input : inputs)
+ for (unsigned int n = 0; n < nComponents; n++)
{
- auto inputPoint = transforms[input]->TransformPoint(outputPoint);
- input->TransformPhysicalPointToContinuousIndex(inputPoint, inputIndices[input]);
- }
+ ComponentType component = InterpolatorConvertType::GetNthComponent( n, value );
- while ( !outIt.IsAtEndOfLine() )
+ if ( component < minComponent )
{
- std::vector<PixelType> pixvals;
-
- for (const auto& input : inputs)
+ PixelConvertType::SetNthComponent( n, outputValue, static_cast<PixelComponentType>( minComponent ) );
+ }
+ else if ( component > maxComponent )
{
- 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];
+ PixelConvertType::SetNthComponent( n, outputValue, static_cast<PixelComponentType>( maxComponent ) );
}
-
- // 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;
+ else
+ {
+ PixelConvertType::SetNthComponent(n, outputValue,
+ static_cast<PixelComponentType>( component ) );
}
- progress.CompletedPixel();
- outIt.NextLine();
}
+
+ return outputValue;
}
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();
+ 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())
{
- m_Interpolators[input] = LinearInterpolatorType::New().GetPointer();
+ 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