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<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);
+ 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<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;
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 <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_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<PixelType>::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<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)
{
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<const InputImageType*, typedef InputImageType::IndexType> lowerIndices;
- std::map<const InputImageType*, typedef InputImageType::IndexType> upperIndices;
+ std::map<const InputImageType*, typename InputImageType::IndexType> lowerIndices;
+ std::map<const InputImageType*, typename InputImageType::IndexType> 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<ComponentType>(minValue);
const ComponentType maxOutputValue = static_cast<ComponentType>(maxValue);
// Walk the output region
outIt.GoToBegin();
while (!outIt.IsAtEnd())
{
// Determine the index of the current output pixel
outputPtr->TransformIndexToPhysicalPoint(outIt.GetIndex(), outputPoint);
std::vector<PixelType> pixvals;
std::vector<double> 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<double>::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<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>
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())
{
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 <itkInterpolateImageFunction.h>
#include <itkNearestNeighborInterpolateImageFunction.h>
#include <itkLinearInterpolateImageFunction.h>
#include <itkBSplineInterpolateImageFunction.h>
#include <itkWindowedSincInterpolateImageFunction.h>
#include <mitkImageAccessByItk.h>
#include <mitkImageCast.h>
#include <mitkGeometry3D.h>
#include <mitkImageToItk.h>
#include <mitkMAPAlgorithmHelper.h>
#include "mapRegistration.h"
#include "mitkRegistrationHelper.h"
template <typename TImage >
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<TImage>::New();
break;
}
case mitk::ImageMappingInterpolator::BSpline_3:
{
typename ::itk::BSplineInterpolateImageFunction<TImage>::Pointer spInterpolator = ::itk::BSplineInterpolateImageFunction<TImage>::New();
spInterpolator->SetSplineOrder(3);
result = spInterpolator;
break;
}
case mitk::ImageMappingInterpolator::WSinc_Hamming:
{
result = ::itk::WindowedSincInterpolateImageFunction<TImage,4>::New();
break;
}
case mitk::ImageMappingInterpolator::WSinc_Welch:
{
result = ::itk::WindowedSincInterpolateImageFunction<TImage,4,::itk::Function::WelchWindowFunction<4> >::New();
break;
}
default:
{
result = ::itk::LinearInterpolateImageFunction<TImage>::New();
break;
}
}
return result;
};
template <typename TPixelType, unsigned int VImageDimension >
-void doMITKStitching(const ::itk::Image<TPixelType,VImageDimension>* input1,
+void doMITKStitching(const ::itk::Image<TPixelType,VImageDimension>* /*input1*/,
mitk::Image::Pointer& result,
std::vector<mitk::Image::ConstPointer> 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<VImageDimension, VImageDimension>;
using ItkImageType = itk::Image<TPixelType, VImageDimension>;
using StitchingFilterType = ::itk::StitchImageFilter<ItkImageType, ItkImageType>;
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<TPixelType, VImageDimension>(*inputIter);
auto castedReg = dynamic_cast<const ConcreteRegistrationType*>(regIter->GetPointer());
auto kernel = dynamic_cast<const ::map::core::RegistrationKernel<VImageDimension, VImageDimension>* >(&(castedReg->getInverseMapping()));
if (nullptr == kernel)
{
mitkThrow() << "Cannot stitch images. At least passed registration object #"<<index<<" doesn't have a valid inverse mapping registration kernel.";
}
stitcher->SetInput(index, itkInput, kernel->getTransformModel(), generateInterpolator< ::itk::Image<TPixelType, VImageDimension> >(interpolatorType));
++inputIter;
++regIter;
++index;
}
stitcher->Update();
mitk::CastToMitkImage<>(stitcher->GetOutput(),result);
}
mitk::Image::Pointer
mitk::StitchImages(std::vector<Image::ConstPointer> 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<Image::ConstPointer> inputs,
std::vector<MAPRegistrationWrapper::ConstPointer> 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<Image::ConstPointer> 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;
}