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; private: using FilterType = itk::StitchImageFilter; FilterType::Pointer m_Filter; InputImageType::Pointer m_Input1; InputImageType::Pointer m_Input2; InputImageType::Pointer m_Input3; public: void setUp() override { InputImageType::PointType origin; origin.Fill(0.); 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& pixels) + { + bool result = true; + + itk::ImageRegionConstIteratorWithIndex iter(image, image->GetLargestPossibleRegion()); + auto refIter = pixels.begin(); + iter.GoToBegin(); + while (!iter.IsAtEnd()) + { + if (refIter == pixels.end()) + { + std::cerr << "Error image to check has a different pixel count then the reference pixel value vector."<SetInput(0, m_Input1); m_Filter->SetInput(1, m_Input2); m_Filter->SetInput(2, m_Input3); m_Filter->Update(); auto output = m_Filter->GetOutput(); + + CPPUNIT_ASSERT(CheckPixels(output, {1, 2, 3, 4, 5, 6, 8.5, 14, 19.5, 40, 50, 60, 85, 140, 195, 400, 500, 600, 700, 800, 900, 1000, 1000, 1000, 1000, 1000, 1000})); + } + + void StitchWithNoInterp() + { + m_Filter->SetInput(0, m_Input1); + + using TranslationType = itk::TranslationTransform; + TranslationType::OutputVectorType offset; + offset[0] = 0; + offset[1] = -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; + 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::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 DecoratedTransformType; typedef typename DecoratedTransformType::Pointer DecoratedTransformPointer; /** Interpolator typedef. */ typedef InterpolateImageFunction< InputImageType, TInterpolatorPrecisionType > InterpolatorType; typedef typename InterpolatorType::Pointer InterpolatorPointerType; typedef typename InterpolatorType::OutputType InterpolatorOutputType; typedef DefaultConvertPixelTraits< InterpolatorOutputType > InterpolatorConvertType; typedef typename InterpolatorConvertType::ComponentType ComponentType; typedef LinearInterpolateImageFunction< InputImageType, TInterpolatorPrecisionType > LinearInterpolatorType; typedef typename LinearInterpolatorType::Pointer LinearInterpolatorPointerType; /** Image size typedef. */ typedef Size< itkGetStaticConstMacro(ImageDimension) > SizeType; /** Image index typedef. */ typedef typename TOutputImage::IndexType IndexType; /** Image point typedef. */ typedef typename InterpolatorType::PointType PointType; //typedef typename TOutputImage::PointType PointType; /** Image pixel value typedef. */ typedef typename TOutputImage::PixelType PixelType; typedef typename TInputImage::PixelType InputPixelType; typedef DefaultConvertPixelTraits PixelConvertType; typedef typename PixelConvertType::ComponentType PixelComponentType; /** Input pixel continuous index typdef */ typedef ContinuousIndex< TTransformPrecisionType, ImageDimension > ContinuousInputIndexType; /** Typedef to describe the output image region type. */ typedef typename TOutputImage::RegionType OutputImageRegionType; /** Image spacing,origin and direction typedef */ typedef typename TOutputImage::SpacingType SpacingType; typedef typename TOutputImage::PointType OriginPointType; typedef typename TOutputImage::DirectionType DirectionType; using Superclass::GetInput; /** Typedef the reference image type to be the ImageBase of the OutputImageType */ typedef ImageBase ReferenceImageBaseType; 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 InputImageVectorType; typedef std::map TransformMapType; typedef std::map InterpolatorMapType; InputImageVectorType GetInputs(); TransformMapType GetTransforms(); InterpolatorMapType m_Interpolators; // Image function for // interpolation PixelType m_DefaultPixelValue; // default pixel value // if the point is // outside the image SizeType m_Size; // Size of the output image SpacingType m_OutputSpacing; // output image spacing OriginPointType m_OutputOrigin; // output image origin DirectionType m_OutputDirection; // output image direction cosines IndexType m_OutputStartIndex; // output image start index bool m_UseReferenceImage; }; } // 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 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::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(transform)); this->ProcessObject::SetInput(transformName, newInput); } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > const typename StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType >::TransformType* StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GetTransform(unsigned int index) const { typedef SimpleDataObjectDecorator< TransformPointerType > DecoratorType; const DecoratorType* input = itkDynamicCastInDebugMode< const DecoratorType* >(this->ProcessObject::GetInput(this->GetTransformInputName(index))); if (nullptr != input) { return input->Get(); } return nullptr; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > const typename StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType >::InterpolatorType* StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GetInterpolator(unsigned int index) const { auto input = this->GetInput(index); if (m_Interpolators.find(input) != std::end(m_Interpolators)) { return m_Interpolators[input]; } return nullptr; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetOutputSpacing(const double *spacing) { SpacingType s; for(unsigned int i = 0; i < TOutputImage::ImageDimension; ++i) { s[i] = static_cast< typename SpacingType::ValueType >(spacing[i]); } this->SetOutputSpacing(s); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetOutputOrigin(const double *origin) { OriginPointType p(origin); this->SetOutputOrigin(p); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetOutputParametersFromImage(const ImageBaseType *image) { this->SetOutputOrigin ( image->GetOrigin() ); this->SetOutputSpacing ( image->GetSpacing() ); this->SetOutputDirection ( image->GetDirection() ); this->SetOutputStartIndex ( image->GetLargestPossibleRegion().GetIndex() ); this->SetSize ( image->GetLargestPossibleRegion().GetSize() ); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::BeforeThreadedGenerateData() { this->EnsureInterpolators(); this->EnsureTransforms(); + for (const auto& interpolator : m_Interpolators) + { + interpolator.second->SetInputImage(interpolator.first); + } + unsigned int nComponents = DefaultConvertPixelTraits::GetNumberOfComponents( m_DefaultPixelValue ); if (nComponents == 0) { PixelComponentType zeroComponent = NumericTraits::ZeroValue( zeroComponent ); nComponents = this->GetInput()->GetNumberOfComponentsPerPixel(); NumericTraits::SetLength(m_DefaultPixelValue, nComponents ); for (unsigned int n=0; n void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::AfterThreadedGenerateData() { // Disconnect input image from the interpolator for (auto& interpolator : m_Interpolators) { interpolator.second->SetInputImage(ITK_NULLPTR); } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType threadId) { - // 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::SetLength( outputValue, nComponents ); - - for (unsigned int n = 0; n < nComponents; n++) - { - ComponentType component = InterpolatorConvertType::GetNthComponent( n, value ); - - if ( component < minComponent ) - { - PixelConvertType::SetNthComponent( n, outputValue, static_cast( minComponent ) ); - } - else if ( component > maxComponent ) - { - PixelConvertType::SetNthComponent( n, outputValue, static_cast( maxComponent ) ); - } - else - { - PixelConvertType::SetNthComponent(n, outputValue, - static_cast( component ) ); - } - } - - return outputValue; -} - -template< typename 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(minValue); + const ComponentType maxOutputValue = static_cast(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 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 inputIndices; - - typedef typename PointType::VectorType VectorType; - std::map 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::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 pixvals; - - for (const auto& input : inputs) + PixelConvertType::SetNthComponent( n, outputValue, static_cast( 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( 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( component ) ); } - progress.CompletedPixel(); - outIt.NextLine(); } + + return outputValue; } template typename StitchImageFilter::InputImageVectorType StitchImageFilter ::GetInputs() { InputImageVectorType inputs; for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); ++i) { auto input = this->GetInput(i); if (nullptr != input) { inputs.push_back(input); } } return inputs; } template typename StitchImageFilter::TransformMapType StitchImageFilter ::GetTransforms() { TransformMapType transforms; for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); ++i) { auto input = this->GetInput(i); auto transform = this->GetTransform(i); transforms[input] = transform; } return transforms; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GenerateInputRequestedRegion() { // Call the superclass' implementation of this method Superclass::GenerateInputRequestedRegion(); if ( !this->GetInput() ) { return; } // Get pointers to the input auto inputs = this->GetInputs(); for (auto& input : inputs) { InputImagePointer inputPtr = const_cast(input); // Determining the actual input region is non-trivial, especially // when we cannot assume anything about the transform being used. // So we do the easy thing and request the entire input image. // inputPtr->SetRequestedRegionToLargestPossibleRegion(); } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GenerateOutputInformation() { // Call the superclass' implementation of this method Superclass::GenerateOutputInformation(); // Get pointers to the input and output OutputImageType *outputPtr = this->GetOutput(); if ( !outputPtr ) { return; } const ReferenceImageBaseType *referenceImage = this->GetReferenceImage(); // Set the size of the output region if ( m_UseReferenceImage && referenceImage ) { outputPtr->SetLargestPossibleRegion( referenceImage->GetLargestPossibleRegion() ); } else { typename TOutputImage::RegionType outputLargestPossibleRegion; outputLargestPossibleRegion.SetSize(m_Size); outputLargestPossibleRegion.SetIndex(m_OutputStartIndex); outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion); } // Set spacing and origin if ( m_UseReferenceImage && referenceImage ) { outputPtr->SetSpacing( referenceImage->GetSpacing() ); outputPtr->SetOrigin( referenceImage->GetOrigin() ); outputPtr->SetDirection( referenceImage->GetDirection() ); } else { outputPtr->SetSpacing(m_OutputSpacing); outputPtr->SetOrigin(m_OutputOrigin); outputPtr->SetDirection(m_OutputDirection); } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > ModifiedTimeType StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GetMTime(void) const { ModifiedTimeType latestTime = Object::GetMTime(); for (const auto& interpolator : m_Interpolators) { if (interpolator.second.GetPointer()) { if (latestTime < interpolator.second->GetMTime()) { latestTime = interpolator.second->GetMTime(); } } } return latestTime; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << "DefaultPixelValue: " << static_cast< typename NumericTraits< PixelType >::PrintType > ( m_DefaultPixelValue ) << std::endl; os << indent << "Size: " << m_Size << std::endl; os << indent << "OutputStartIndex: " << m_OutputStartIndex << std::endl; os << indent << "OutputSpacing: " << m_OutputSpacing << std::endl; os << indent << "OutputOrigin: " << m_OutputOrigin << std::endl; os << indent << "OutputDirection: " << m_OutputDirection << std::endl; for (const auto& interpolator : m_Interpolators) { os << indent << "Interpolator: " << interpolator.second.GetPointer() << std::endl; } os << indent << "UseReferenceImage: " << ( m_UseReferenceImage ? "On" : "Off" ) << std::endl; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::EnsureTransforms() { const auto inputCount = this->GetNumberOfIndexedInputs(); for (unsigned int i = 0; i < inputCount; ++i) { auto input = this->GetInput(i); if (nullptr == input) { itkExceptionMacro(<< "Nth input image is not set (n: " << i << ")."); } auto transform = this->GetTransform(i); if (nullptr == transform) { this->SetTransform(i, itk::IdentityTransform< TTransformPrecisionType, ImageDimension>::New().GetPointer()); } } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::EnsureInterpolators() { const auto inputCount = this->GetNumberOfIndexedInputs(); + InterpolatorMapType newInterpolatorMap; + for (unsigned int i = 0; i < inputCount; ++i) { auto input = this->GetInput(i); if (nullptr == input) { itkExceptionMacro(<< "Nth input image is not set (n: " << i << ")."); } if (m_Interpolators[input].IsNull()) { - 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