diff --git a/Modules/MatchPointRegistration/include/itkStitchImageFilter.h b/Modules/MatchPointRegistration/include/itkStitchImageFilter.h index e61e31b905..3519fbacea 100644 --- a/Modules/MatchPointRegistration/include/itkStitchImageFilter.h +++ b/Modules/MatchPointRegistration/include/itkStitchImageFilter.h @@ -1,328 +1,330 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef itkStitchImageFilter_h #define itkStitchImageFilter_h #include "itkFixedArray.h" #include "itkTransform.h" #include "itkImageRegionIterator.h" #include "itkImageToImageFilter.h" #include "itkLinearInterpolateImageFunction.h" #include "itkSize.h" #include "itkDefaultConvertPixelTraits.h" #include "itkDataObjectDecorator.h" namespace itk { enum class StitchStrategy { Mean = 0, //use the mean value of all inputs that can provide a pixel vaule BorderDistance = 1 //use the value that is largest minimal distance to its image borders (use e.g. if vaules tend to be not relyable at borders) }; std::ostream& operator<< (std::ostream& os, const itk::StitchStrategy& strategy) { if (itk::StitchStrategy::Mean == strategy) os << "Mean"; else if (itk::StitchStrategy::BorderDistance == strategy) os << "BorderDistance"; else os << "unkown"; return os; }; /** \class StitchImageFilter * \brief ITK filter that resamples/stitches multiple images into a given reference geometry. * * StitchImageFilter is 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 the behavior depends on * the StitchStragy: * - Mean: a wighted sum of all voxels mapped input pixel values will be calculated. * - BorderDisntance: the voxel will be choosen that have the largest minimal distance to its own image borders. * * 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); + using Superclass::SetInput; + void SetInput(const InputImageType* image) override; + void SetInput(unsigned int index, const InputImageType* image) override; /** Convinience methods that allows setting of input image and its transform in one call.*/ virtual void SetInput(unsigned int index, const InputImageType* image, const TransformType* transform); virtual void SetInput(unsigned int index, const InputImageType* image, const TransformType* transform, InterpolatorType* interpolator); const TransformType* GetTransform(unsigned int index) const; const InterpolatorType* GetInterpolator(unsigned int index) const; /** Get/Set the size of the output image. */ itkSetMacro(Size, SizeType); itkGetConstReferenceMacro(Size, SizeType); /** Get/Set the pixel value when a transformed pixel is outside of the * image. The default default pixel value is 0. */ itkSetMacro(DefaultPixelValue, PixelType); itkGetConstReferenceMacro(DefaultPixelValue, PixelType); /** Set the output image spacing. */ itkSetMacro(OutputSpacing, SpacingType); virtual void SetOutputSpacing(const double *values); /** Get the output image spacing. */ itkGetConstReferenceMacro(OutputSpacing, SpacingType); /** Set the output image origin. */ itkSetMacro(OutputOrigin, OriginPointType); virtual void SetOutputOrigin(const double *values); /** Get the output image origin. */ itkGetConstReferenceMacro(OutputOrigin, OriginPointType); /** Set the output direciton cosine matrix. */ itkSetMacro(OutputDirection, DirectionType); itkGetConstReferenceMacro(OutputDirection, DirectionType); /** Helper method to set the output parameters based on this image. */ void SetOutputParametersFromImage(const ImageBaseType *image); /** Set the start index of the output largest possible region. * The default is an index of all zeros. */ itkSetMacro(OutputStartIndex, IndexType); /** Get the start index of the output largest possible region. */ itkGetConstReferenceMacro(OutputStartIndex, IndexType); /** Set a reference image to use to define the output information. * By default, output information is specificed through the * SetOutputSpacing, Origin, and Direction methods. Alternatively, * this method can be used to specify an image from which to * copy the information. UseReferenceImageOn must be set to utilize the * reference image. */ itkSetInputMacro(ReferenceImage, ReferenceImageBaseType); /** Get the reference image that is defining the output information. */ itkGetInputMacro(ReferenceImage, ReferenceImageBaseType); /** Turn on/off whether a specified reference image should be used to define * the output information. */ itkSetMacro(UseReferenceImage, bool); itkBooleanMacro(UseReferenceImage); itkGetConstMacro(UseReferenceImage, bool); itkSetMacro(StitchStrategy, StitchStrategy); itkGetConstMacro(StitchStrategy, StitchStrategy); /** StitchImageFilter produces an image which is a different size * than its input. As such, it needs to provide an implementation * for GenerateOutputInformation() in order to inform the pipeline * execution model. The original documentation of this method is * below. \sa ProcessObject::GenerateOutputInformaton() */ virtual void GenerateOutputInformation() ITK_OVERRIDE; /** StitchImageFilter needs a different input requested region than * the output requested region. As such, StitchImageFilter needs * to provide an implementation for GenerateInputRequestedRegion() * in order to inform the pipeline execution model. * \sa ProcessObject::GenerateInputRequestedRegion() */ virtual void GenerateInputRequestedRegion() ITK_OVERRIDE; /** Set up state of filter before multi-threading. * InterpolatorType::SetInputImage is not thread-safe and hence * has to be set up before ThreadedGenerateData */ virtual void BeforeThreadedGenerateData() ITK_OVERRIDE; /** Set the state of the filter after multi-threading. */ virtual void AfterThreadedGenerateData() ITK_OVERRIDE; /** Compute the Modified Time based on the changed components. */ ModifiedTimeType GetMTime(void) const ITK_OVERRIDE; #ifdef ITK_USE_CONCEPT_CHECKING // Begin concept checking itkConceptMacro( OutputHasNumericTraitsCheck, ( Concept::HasNumericTraits< PixelComponentType > ) ); // End concept checking #endif protected: StitchImageFilter(); ~StitchImageFilter() ITK_OVERRIDE {} void PrintSelf(std::ostream & os, Indent indent) const ITK_OVERRIDE; /** Override VeriyInputInformation() since this filter's inputs do * not need to occoupy the same physical space. * * \sa ProcessObject::VerifyInputInformation */ virtual void VerifyInputInformation() ITK_OVERRIDE { } /** StitchImageFilter can be implemented as a multithreaded filter. * Therefore, this implementation provides a ThreadedGenerateData() * routine which is called for each processing thread. The output * image data is allocated automatically by the superclass prior * to calling ThreadedGenerateData(). * ThreadedGenerateData can only write to the portion of the output image * specified by the parameter "outputRegionForThread" * \sa ImageToImageFilter::ThreadedGenerateData(), * ImageToImageFilter::GenerateData() */ virtual void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType threadId) ITK_OVERRIDE; /** Cast pixel from interpolator output to PixelType. */ virtual PixelType CastPixelWithBoundsChecking( const InterpolatorOutputType value, const ComponentType minComponent, const ComponentType maxComponent) const; void SetTransform(unsigned int index, const TransformType* transform); /** Helper that ensures that a transform is specified for every input image. If a input image has no specified transforms, an identity transform will be created and set as default.*/ void EnsureTransforms(); /** Helper that ensures that an interpolator is specified for every input image. If a input image has no specified interpolator, a linear interpolator will be created and set as default.*/ void EnsureInterpolators(); static std::string GetTransformInputName(unsigned int index); private: ITK_DISALLOW_COPY_AND_ASSIGN(StitchImageFilter); typedef std::vector InputImageVectorType; typedef std::map TransformMapType; typedef std::map InterpolatorMapType; InputImageVectorType GetInputs(); TransformMapType GetTransforms(); InterpolatorMapType m_Interpolators; // Image function for // interpolation PixelType m_DefaultPixelValue; // default pixel value // if the point is // outside the image SizeType m_Size; // Size of the output image SpacingType m_OutputSpacing; // output image spacing OriginPointType m_OutputOrigin; // output image origin DirectionType m_OutputDirection; // output image direction cosines IndexType m_OutputStartIndex; // output image start index bool m_UseReferenceImage; StitchStrategy m_StitchStrategy; }; } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "itkStitchImageFilter.tpp" #endif #endif diff --git a/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp b/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp index fbff61d979..1ec8334ada 100644 --- a/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp +++ b/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp @@ -1,628 +1,639 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef itkStitchImageFilter_hxx #define itkStitchImageFilter_hxx #include "itkStitchImageFilter.h" #include "itkObjectFactory.h" #include "itkIdentityTransform.h" #include "itkProgressReporter.h" #include "itkImageRegionIteratorWithIndex.h" #include "itkImageScanlineIterator.h" #include "itkSpecialCoordinatesImage.h" #include "itkDefaultConvertPixelTraits.h" #include "itkSimpleDataObjectDecorator.h" #include namespace itk { template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::StitchImageFilter() : m_OutputSpacing( 1.0 ), m_OutputOrigin( 0.0 ), m_UseReferenceImage( false ), m_StitchStrategy(StitchStrategy::Mean) { m_Size.Fill( 0 ); m_OutputStartIndex.Fill( 0 ); m_OutputDirection.SetIdentity(); // Pipeline input configuration // implicit input index set: // #1 "ReferenceImage" optional Self::AddOptionalInputName("ReferenceImage"); m_DefaultPixelValue = NumericTraits::ZeroValue( m_DefaultPixelValue ); } +template< typename TInputImage, + typename TOutputImage, + typename TInterpolatorPrecisionType, + typename TTransformPrecisionType > + void + StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > + ::SetInput(const InputImageType* image) +{ + this->SetInput(0, image, itk::IdentityTransform< TTransformPrecisionType, ImageDimension>::New().GetPointer(), LinearInterpolatorType::New().GetPointer()); +} + template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetInput(unsigned int index, const InputImageType* image) { this->SetInput(index, image, itk::IdentityTransform< TTransformPrecisionType, ImageDimension>::New().GetPointer(), LinearInterpolatorType::New().GetPointer()); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetInput(unsigned int index, const InputImageType* image, const TransformType* transform) { this->SetInput(index, image, transform, LinearInterpolatorType::New().GetPointer()); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetInput(unsigned int index, const InputImageType* image, const TransformType* transform, InterpolatorType* interpolator) { Superclass::SetInput(index, image); m_Interpolators[image] = interpolator; this->SetTransform(index, transform); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetTransform(unsigned int index, const TransformType* transform) { const auto transformName = this->GetTransformInputName(index); typedef SimpleDataObjectDecorator< TransformPointerType > DecoratorType; const DecoratorType* oldInput = itkDynamicCastInDebugMode< const DecoratorType* >(this->ProcessObject::GetInput(transformName)); if (!oldInput || oldInput->Get() != transform) { typename DecoratorType::Pointer newInput = DecoratorType::New(); // Process object is not const-correct so the const_cast is required here newInput->Set(const_cast(transform)); this->ProcessObject::SetInput(transformName, newInput); } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > const typename StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType >::TransformType* StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GetTransform(unsigned int index) const { typedef SimpleDataObjectDecorator< TransformPointerType > DecoratorType; const DecoratorType* input = itkDynamicCastInDebugMode< const DecoratorType* >(this->ProcessObject::GetInput(this->GetTransformInputName(index))); if (nullptr != input) { return input->Get(); } return nullptr; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > const typename StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType >::InterpolatorType* StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GetInterpolator(unsigned int index) const { auto input = this->GetInput(index); if (m_Interpolators.find(input) != std::end(m_Interpolators)) { return m_Interpolators[input]; } return nullptr; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetOutputSpacing(const double *spacing) { SpacingType s; for(unsigned int i = 0; i < TOutputImage::ImageDimension; ++i) { s[i] = static_cast< typename SpacingType::ValueType >(spacing[i]); } this->SetOutputSpacing(s); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetOutputOrigin(const double *origin) { OriginPointType p(origin); this->SetOutputOrigin(p); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetOutputParametersFromImage(const ImageBaseType *image) { this->SetOutputOrigin ( image->GetOrigin() ); this->SetOutputSpacing ( image->GetSpacing() ); this->SetOutputDirection ( image->GetDirection() ); this->SetOutputStartIndex ( image->GetLargestPossibleRegion().GetIndex() ); this->SetSize ( image->GetLargestPossibleRegion().GetSize() ); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::BeforeThreadedGenerateData() { this->EnsureInterpolators(); this->EnsureTransforms(); for (const auto& interpolator : m_Interpolators) { interpolator.second->SetInputImage(interpolator.first); } unsigned int nComponents = DefaultConvertPixelTraits::GetNumberOfComponents( m_DefaultPixelValue ); if (nComponents == 0) { PixelComponentType zeroComponent = NumericTraits::ZeroValue( zeroComponent ); nComponents = this->GetInput()->GetNumberOfComponentsPerPixel(); NumericTraits::SetLength(m_DefaultPixelValue, nComponents ); for (unsigned int n=0; n void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::AfterThreadedGenerateData() { // Disconnect input image from the interpolator for (auto& interpolator : m_Interpolators) { interpolator.second->SetInputImage(ITK_NULLPTR); } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType threadId) { if( outputRegionForThread.GetNumberOfPixels() == 0 ) { return; } // Get the output pointers OutputImageType* outputPtr = this->GetOutput(); // Get this input pointers InputImageVectorType inputs = this->GetInputs(); TransformMapType transforms = this->GetTransforms(); - std::map lowerIndices; - std::map upperIndices; + std::map lowerIndices; + std::map upperIndices; for (const auto& input : inputs) { const auto largestRegion = input->GetLargestPossibleRegion(); lowerIndices[input] = largestRegion.GetIndex(); upperIndices[input] = largestRegion.GetUpperIndex(); } // Create an iterator that will walk the output region for this thread. typedef ImageRegionIteratorWithIndex< OutputImageType > OutputIterator; OutputIterator outIt(outputPtr, outputRegionForThread); // Define a few indices that will be used to translate from an input pixel // to an output pixel PointType outputPoint; // Coordinates of current output pixel PointType inputPoint; // Coordinates of current input pixel ContinuousInputIndexType inputIndex; // Support for progress methods/callbacks ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); // Min/max values of the output pixel type AND these values // represented as the output type of the interpolator const PixelComponentType minValue = NumericTraits< PixelComponentType >::NonpositiveMin(); const PixelComponentType maxValue = NumericTraits< PixelComponentType >::max(); typedef typename InterpolatorType::OutputType OutputType; const ComponentType minOutputValue = static_cast(minValue); const ComponentType maxOutputValue = static_cast(maxValue); // Walk the output region outIt.GoToBegin(); while (!outIt.IsAtEnd()) { // Determine the index of the current output pixel outputPtr->TransformIndexToPhysicalPoint(outIt.GetIndex(), outputPoint); std::vector pixvals; std::vector pixDistance; for (const auto& input : inputs) { // Compute corresponding input pixel position inputPoint = transforms[input]->TransformPoint(outputPoint); const bool isInsideInput = input->TransformPhysicalPointToContinuousIndex(inputPoint, inputIndex); // Evaluate input at right position and copy to the output if (m_Interpolators[input]->IsInsideBuffer(inputIndex) && isInsideInput) { OutputType value = m_Interpolators[input]->EvaluateAtContinuousIndex(inputIndex); pixvals.emplace_back(this->CastPixelWithBoundsChecking(value, minOutputValue, maxOutputValue)); ContinuousInputIndexType indexDistance; const auto spacing = input->GetSpacing(); double minBorderDistance = std::numeric_limits::max(); for (unsigned int i = 0; i < ImageDimension; ++i) { minBorderDistance = std::min(minBorderDistance, std::min(std::abs(lowerIndices[input][i] - inputIndex[i]) * spacing[i], std::abs(upperIndices[input][i] - inputIndex[i]) * spacing[i])); } pixDistance.emplace_back(minBorderDistance); } } if (!pixvals.empty()) { //at least one input provided a value if (StitchStrategy::Mean == m_StitchStrategy) { double sum = std::accumulate(pixvals.begin(), pixvals.end(), 0.0); outIt.Set(sum / pixvals.size()); } else { auto finding = std::max_element(pixDistance.begin(), pixDistance.end()); outIt.Set(pixvals[std::distance(pixDistance.begin(), finding)]); } } else { outIt.Set(m_DefaultPixelValue); // default background value } progress.CompletedPixel(); ++outIt; } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > typename StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::PixelType StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::CastPixelWithBoundsChecking(const InterpolatorOutputType value, const ComponentType minComponent, const ComponentType maxComponent ) const { const unsigned int nComponents = InterpolatorConvertType::GetNumberOfComponents(value); PixelType outputValue; NumericTraits::SetLength( outputValue, nComponents ); for (unsigned int n = 0; n < nComponents; n++) { ComponentType component = InterpolatorConvertType::GetNthComponent( n, value ); if ( component < minComponent ) { PixelConvertType::SetNthComponent( n, outputValue, static_cast( minComponent ) ); } else if ( component > maxComponent ) { PixelConvertType::SetNthComponent( n, outputValue, static_cast( maxComponent ) ); } else { PixelConvertType::SetNthComponent(n, outputValue, static_cast( component ) ); } } return outputValue; } template typename StitchImageFilter::InputImageVectorType StitchImageFilter ::GetInputs() { InputImageVectorType inputs; for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); ++i) { auto input = this->GetInput(i); if (nullptr != input) { inputs.push_back(input); } } return inputs; } template typename StitchImageFilter::TransformMapType StitchImageFilter ::GetTransforms() { TransformMapType transforms; for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); ++i) { auto input = this->GetInput(i); auto transform = this->GetTransform(i); transforms[input] = transform; } return transforms; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GenerateInputRequestedRegion() { // Call the superclass' implementation of this method Superclass::GenerateInputRequestedRegion(); if ( !this->GetInput() ) { return; } // Get pointers to the input auto inputs = this->GetInputs(); for (auto& input : inputs) { InputImagePointer inputPtr = const_cast(input); // Determining the actual input region is non-trivial, especially // when we cannot assume anything about the transform being used. // So we do the easy thing and request the entire input image. // inputPtr->SetRequestedRegionToLargestPossibleRegion(); } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GenerateOutputInformation() { // Call the superclass' implementation of this method Superclass::GenerateOutputInformation(); // Get pointers to the input and output OutputImageType *outputPtr = this->GetOutput(); if ( !outputPtr ) { return; } const ReferenceImageBaseType *referenceImage = this->GetReferenceImage(); // Set the size of the output region if ( m_UseReferenceImage && referenceImage ) { outputPtr->SetLargestPossibleRegion( referenceImage->GetLargestPossibleRegion() ); } else { typename TOutputImage::RegionType outputLargestPossibleRegion; outputLargestPossibleRegion.SetSize(m_Size); outputLargestPossibleRegion.SetIndex(m_OutputStartIndex); outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion); } // Set spacing and origin if ( m_UseReferenceImage && referenceImage ) { outputPtr->SetSpacing( referenceImage->GetSpacing() ); outputPtr->SetOrigin( referenceImage->GetOrigin() ); outputPtr->SetDirection( referenceImage->GetDirection() ); } else { outputPtr->SetSpacing(m_OutputSpacing); outputPtr->SetOrigin(m_OutputOrigin); outputPtr->SetDirection(m_OutputDirection); } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > ModifiedTimeType StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GetMTime(void) const { ModifiedTimeType latestTime = Object::GetMTime(); for (const auto& interpolator : m_Interpolators) { if (interpolator.second.GetPointer()) { if (latestTime < interpolator.second->GetMTime()) { latestTime = interpolator.second->GetMTime(); } } } return latestTime; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << "DefaultPixelValue: " << static_cast< typename NumericTraits< PixelType >::PrintType > ( m_DefaultPixelValue ) << std::endl; os << indent << "Size: " << m_Size << std::endl; os << indent << "OutputStartIndex: " << m_OutputStartIndex << std::endl; os << indent << "OutputSpacing: " << m_OutputSpacing << std::endl; os << indent << "OutputOrigin: " << m_OutputOrigin << std::endl; os << indent << "OutputDirection: " << m_OutputDirection << std::endl; for (const auto& interpolator : m_Interpolators) { os << indent << "Interpolator: " << interpolator.second.GetPointer() << std::endl; } os << indent << "UseReferenceImage: " << ( m_UseReferenceImage ? "On" : "Off" ) << std::endl; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::EnsureTransforms() { const auto inputCount = this->GetNumberOfIndexedInputs(); for (unsigned int i = 0; i < inputCount; ++i) { auto input = this->GetInput(i); if (nullptr == input) { itkExceptionMacro(<< "Nth input image is not set (n: " << i << ")."); } auto transform = this->GetTransform(i); if (nullptr == transform) { this->SetTransform(i, itk::IdentityTransform< TTransformPrecisionType, ImageDimension>::New().GetPointer()); } } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::EnsureInterpolators() { const auto inputCount = this->GetNumberOfIndexedInputs(); InterpolatorMapType newInterpolatorMap; for (unsigned int i = 0; i < inputCount; ++i) { auto input = this->GetInput(i); if (nullptr == input) { itkExceptionMacro(<< "Nth input image is not set (n: " << i << ")."); } if (m_Interpolators[input].IsNull()) { newInterpolatorMap[input] = LinearInterpolatorType::New().GetPointer(); } else { newInterpolatorMap[input] = m_Interpolators[input]; } } m_Interpolators = newInterpolatorMap; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > std::string StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GetTransformInputName(unsigned int index) { return "transform_" + std::to_string(index); } } // end namespace itk #endif diff --git a/Modules/MatchPointRegistration/src/Helper/mitkImageStitchingHelper.cpp b/Modules/MatchPointRegistration/src/Helper/mitkImageStitchingHelper.cpp index b2b139e8ae..e0fcb5f46a 100644 --- a/Modules/MatchPointRegistration/src/Helper/mitkImageStitchingHelper.cpp +++ b/Modules/MatchPointRegistration/src/Helper/mitkImageStitchingHelper.cpp @@ -1,229 +1,229 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkImageStitchingHelper.h" #include #include #include #include #include #include #include #include #include #include #include "mapRegistration.h" #include "mitkRegistrationHelper.h" template typename ::itk::InterpolateImageFunction< TImage >::Pointer generateInterpolator(mitk::ImageMappingInterpolator::Type interpolatorType) { typedef ::itk::InterpolateImageFunction< TImage > BaseInterpolatorType; typename BaseInterpolatorType::Pointer result; switch (interpolatorType) { case mitk::ImageMappingInterpolator::NearestNeighbor: { result = ::itk::NearestNeighborInterpolateImageFunction::New(); break; } case mitk::ImageMappingInterpolator::BSpline_3: { typename ::itk::BSplineInterpolateImageFunction::Pointer spInterpolator = ::itk::BSplineInterpolateImageFunction::New(); spInterpolator->SetSplineOrder(3); result = spInterpolator; break; } case mitk::ImageMappingInterpolator::WSinc_Hamming: { result = ::itk::WindowedSincInterpolateImageFunction::New(); break; } case mitk::ImageMappingInterpolator::WSinc_Welch: { result = ::itk::WindowedSincInterpolateImageFunction >::New(); break; } default: { result = ::itk::LinearInterpolateImageFunction::New(); break; } } return result; }; template -void doMITKStitching(const ::itk::Image* input1, +void doMITKStitching(const ::itk::Image* /*input1*/, mitk::Image::Pointer& result, std::vector inputs, std::vector<::map::core::RegistrationBase::ConstPointer> registrations, const mitk::BaseGeometry* resultGeometry, const double& paddingValue, itk::StitchStrategy stitchStrategy, mitk::ImageMappingInterpolator::Type interpolatorType) { using ConcreteRegistrationType = ::map::core::Registration; using ItkImageType = itk::Image; using StitchingFilterType = ::itk::StitchImageFilter; auto stitcher = StitchingFilterType::New(); stitcher->SetDefaultPixelValue(paddingValue); stitcher->SetOutputOrigin(resultGeometry->GetOrigin()); const auto spacing = resultGeometry->GetSpacing(); stitcher->SetOutputSpacing(spacing); - StitchingFilterType::DirectionType itkDirection; + typename StitchingFilterType::DirectionType itkDirection; const auto mitkDirection = resultGeometry->GetIndexToWorldTransform()->GetMatrix(); for (unsigned int i = 0; i < VImageDimension; ++i) { for (unsigned int j = 0; j < VImageDimension; ++j) { itkDirection[i][j] = mitkDirection[i][j] / spacing[j]; } } stitcher->SetOutputDirection(itkDirection); - ItkImageType::SizeType size; + typename ItkImageType::SizeType size; size[0] = resultGeometry->GetExtent(0); size[1] = resultGeometry->GetExtent(1); size[2] = resultGeometry->GetExtent(2); stitcher->SetSize(size); stitcher->SetNumberOfThreads(1); stitcher->SetStitchStrategy(stitchStrategy); auto inputIter = inputs.begin(); auto regIter = registrations.begin(); unsigned int index = 0; while (inputIter != inputs.end()) { auto itkInput = mitk::ImageToItkImage(*inputIter); auto castedReg = dynamic_cast(regIter->GetPointer()); auto kernel = dynamic_cast* >(&(castedReg->getInverseMapping())); if (nullptr == kernel) { mitkThrow() << "Cannot stitch images. At least passed registration object #"<SetInput(index, itkInput, kernel->getTransformModel(), generateInterpolator< ::itk::Image >(interpolatorType)); ++inputIter; ++regIter; ++index; } stitcher->Update(); mitk::CastToMitkImage<>(stitcher->GetOutput(),result); } mitk::Image::Pointer mitk::StitchImages(std::vector inputs, std::vector<::map::core::RegistrationBase::ConstPointer> registrations, const BaseGeometry* resultGeometry, const double& paddingValue, itk::StitchStrategy stitchStrategy, mitk::ImageMappingInterpolator::Type interpolatorType) { if (inputs.size() != registrations.size()) { mitkThrow() << "Cannot stitch images. Passed inputs vector and registrations vector have different sizes."; } if (inputs.empty()) { mitkThrow() << "Cannot stitch images. No input images are defined."; } auto inputDim = inputs.front()->GetDimension(); auto inputPixelType = inputs.front()->GetPixelType(); for (const auto& input : inputs) { if (input->GetDimension() != inputDim) { mitkThrow() << "Cannot stitch images. Images have different dimensions. Dimeonsion of first input: " << inputDim << "; wrong dimension: " << input->GetDimension(); } if (input->GetPixelType() != inputPixelType) { mitkThrow() << "Cannot stitch images. Input images have different pixeltypes. The current implementation does only support stitching of images with same pixel type. Dimeonsion of first input: " << inputPixelType.GetTypeAsString() << "; wrong dimension: " << input->GetPixelType().GetTypeAsString(); } if (input->GetTimeSteps() > 1) { mitkThrow() << "Cannot stitch dynamic images. At least one input image has multiple time steps."; } } for (const auto& reg : registrations) { if (reg->getMovingDimensions() != inputDim) { mitkThrow() << "Cannot stitch images. At least one registration has a different moving dimension then the inputs. Dimeonsion of inputs: " << inputDim << "; wrong dimension: " << reg->getMovingDimensions(); } if (reg->getTargetDimensions() != inputDim) { mitkThrow() << "Cannot stitch images. At least one registration has a different target dimension then the inputs. Dimeonsion of inputs: " << inputDim << "; wrong dimension: " << reg->getTargetDimensions(); } } Image::Pointer result; AccessFixedDimensionByItk_n(inputs.front(), doMITKStitching, 3, (result, inputs, registrations, resultGeometry, paddingValue, stitchStrategy, interpolatorType)); return result; } mitk::Image::Pointer mitk::StitchImages(std::vector inputs, std::vector registrations, const BaseGeometry* resultGeometry, const double& paddingValue, itk::StitchStrategy stitchStrategy, mitk::ImageMappingInterpolator::Type interpolatorType) { std::vector<::map::core::RegistrationBase::ConstPointer> unwrappedRegs; for (const auto& reg : registrations) { if (!reg) { mitkThrow() << "Cannot stitch images. At least one passed registration wrapper pointer is nullptr."; } unwrappedRegs.push_back(reg->GetRegistration()); } Image::Pointer result = StitchImages(inputs, unwrappedRegs, resultGeometry, paddingValue, stitchStrategy, interpolatorType); return result; } mitk::Image::Pointer mitk::StitchImages(std::vector inputs, const BaseGeometry* resultGeometry, const double& paddingValue, itk::StitchStrategy stitchStrategy, mitk::ImageMappingInterpolator::Type interpolatorType) { auto defaultReg = GenerateIdentityRegistration3D(); std::vector<::map::core::RegistrationBase::ConstPointer> defaultRegs; defaultRegs.resize(inputs.size()); std::fill(defaultRegs.begin(), defaultRegs.end(), defaultReg->GetRegistration()); Image::Pointer result = StitchImages(inputs, defaultRegs, resultGeometry, paddingValue, stitchStrategy, interpolatorType); return result; }