diff --git a/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp b/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp index 98e8f0dc4e..2b535e174f 100644 --- a/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp +++ b/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp @@ -1,638 +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) { + this->DynamicMultiThreadingOff(); 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(); 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; 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)); 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/ModelFit/include/itkMaskedStatisticsImageFilter.hxx b/Modules/ModelFit/include/itkMaskedStatisticsImageFilter.hxx index 3057714651..63a3902636 100644 --- a/Modules/ModelFit/include/itkMaskedStatisticsImageFilter.hxx +++ b/Modules/ModelFit/include/itkMaskedStatisticsImageFilter.hxx @@ -1,393 +1,395 @@ /*========================================================================= * * Copyright Insight Software Consortium * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0.txt * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * *=========================================================================*/ #ifndef __itkMaskedStatisticsImageFilter_hxx #define __itkMaskedStatisticsImageFilter_hxx #include "itkMaskedStatisticsImageFilter.h" #include "itkImageScanlineIterator.h" #include "itkProgressReporter.h" namespace itk { template< typename TInputImage, typename TMaskImage > MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::MaskedStatisticsImageFilter():m_ThreadSum(1), m_SumOfSquares(1), m_Count(1), m_ThreadMin(1), m_ThreadMax(1) { + this->DynamicMultiThreadingOff(); + // first output is a copy of the image, DataObject created by // superclass // // allocate the data objects for the outputs which are // just decorators around pixel types for ( int i = 1; i < 3; ++i ) { typename PixelObjectType::Pointer output = static_cast< PixelObjectType * >( this->MakeOutput(i).GetPointer() ); this->ProcessObject::SetNthOutput( i, output.GetPointer() ); } // allocate the data objects for the outputs which are // just decorators around real types for ( int i = 3; i < 7; ++i ) { typename RealObjectType::Pointer output = static_cast< RealObjectType * >( this->MakeOutput(i).GetPointer() ); this->ProcessObject::SetNthOutput( i, output.GetPointer() ); } this->GetMinimumOutput()->Set( NumericTraits< PixelType >::max() ); this->GetMaximumOutput()->Set( NumericTraits< PixelType >::NonpositiveMin() ); this->GetMeanOutput()->Set( NumericTraits< RealType >::max() ); this->GetSigmaOutput()->Set( NumericTraits< RealType >::max() ); this->GetVarianceOutput()->Set( NumericTraits< RealType >::max() ); this->GetSumOutput()->Set(NumericTraits< RealType >::Zero); } template< typename TInputImage, typename TMaskImage > DataObject::Pointer MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::MakeOutput(DataObjectPointerArraySizeType output) { switch ( output ) { case 0: return TInputImage::New().GetPointer(); break; case 1: return PixelObjectType::New().GetPointer(); break; case 2: return PixelObjectType::New().GetPointer(); break; case 3: case 4: case 5: case 6: return RealObjectType::New().GetPointer(); break; default: // might as well make an image return TInputImage::New().GetPointer(); break; } } template< typename TInputImage, typename TMaskImage > typename MaskedStatisticsImageFilter< TInputImage, TMaskImage >::PixelObjectType * MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::GetMinimumOutput() { return static_cast< PixelObjectType * >( this->ProcessObject::GetOutput(1) ); } template< typename TInputImage, typename TMaskImage > const typename MaskedStatisticsImageFilter< TInputImage, TMaskImage >::PixelObjectType * MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::GetMinimumOutput() const { return static_cast< const PixelObjectType * >( this->ProcessObject::GetOutput(1) ); } template< typename TInputImage, typename TMaskImage > typename MaskedStatisticsImageFilter< TInputImage, TMaskImage >::PixelObjectType * MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::GetMaximumOutput() { return static_cast< PixelObjectType * >( this->ProcessObject::GetOutput(2) ); } template< typename TInputImage, typename TMaskImage > const typename MaskedStatisticsImageFilter< TInputImage, TMaskImage >::PixelObjectType * MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::GetMaximumOutput() const { return static_cast< const PixelObjectType * >( this->ProcessObject::GetOutput(2) ); } template< typename TInputImage, typename TMaskImage > typename MaskedStatisticsImageFilter< TInputImage, TMaskImage >::RealObjectType * MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::GetMeanOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(3) ); } template< typename TInputImage, typename TMaskImage > const typename MaskedStatisticsImageFilter< TInputImage, TMaskImage >::RealObjectType * MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::GetMeanOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(3) ); } template< typename TInputImage, typename TMaskImage > typename MaskedStatisticsImageFilter< TInputImage, TMaskImage >::RealObjectType * MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::GetSigmaOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(4) ); } template< typename TInputImage, typename TMaskImage > const typename MaskedStatisticsImageFilter< TInputImage, TMaskImage >::RealObjectType * MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::GetSigmaOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(4) ); } template< typename TInputImage, typename TMaskImage > typename MaskedStatisticsImageFilter< TInputImage, TMaskImage >::RealObjectType * MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::GetVarianceOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(5) ); } template< typename TInputImage, typename TMaskImage > const typename MaskedStatisticsImageFilter< TInputImage, TMaskImage >::RealObjectType * MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::GetVarianceOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(5) ); } template< typename TInputImage, typename TMaskImage > typename MaskedStatisticsImageFilter< TInputImage, TMaskImage >::RealObjectType * MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::GetSumOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(6) ); } template< typename TInputImage, typename TMaskImage > const typename MaskedStatisticsImageFilter< TInputImage, TMaskImage >::RealObjectType * MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::GetSumOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(6) ); } template< typename TInputImage, typename TMaskImage > void MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); if ( this->GetInput() ) { InputImagePointer image = const_cast< typename Superclass::InputImageType * >( this->GetInput() ); image->SetRequestedRegionToLargestPossibleRegion(); } } template< typename TInputImage, typename TMaskImage > void MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::EnlargeOutputRequestedRegion(DataObject *data) { Superclass::EnlargeOutputRequestedRegion(data); data->SetRequestedRegionToLargestPossibleRegion(); } template< typename TInputImage, typename TMaskImage > void MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::AllocateOutputs() { // Pass the input through as the output InputImagePointer image = const_cast< TInputImage * >( this->GetInput() ); this->GraftOutput(image); // Nothing that needs to be allocated for the remaining outputs } template< typename TInputImage, typename TMaskImage > void MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::BeforeThreadedGenerateData() { ThreadIdType numberOfThreads = this->GetNumberOfWorkUnits(); // Resize the thread temporaries m_Count.SetSize(numberOfThreads); m_SumOfSquares.SetSize(numberOfThreads); m_ThreadSum.SetSize(numberOfThreads); m_ThreadMin.SetSize(numberOfThreads); m_ThreadMax.SetSize(numberOfThreads); // Initialize the temporaries m_Count.Fill(NumericTraits< SizeValueType >::Zero); m_ThreadSum.Fill(NumericTraits< RealType >::Zero); m_SumOfSquares.Fill(NumericTraits< RealType >::Zero); m_ThreadMin.Fill( NumericTraits< PixelType >::max() ); m_ThreadMax.Fill( NumericTraits< PixelType >::NonpositiveMin() ); } template< typename TInputImage, typename TMaskImage > void MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::AfterThreadedGenerateData() { ThreadIdType i; SizeValueType count; RealType sumOfSquares; ThreadIdType numberOfThreads = this->GetNumberOfWorkUnits(); PixelType minimum; PixelType maximum; RealType mean; RealType sigma; RealType variance; RealType sum; sum = sumOfSquares = NumericTraits< RealType >::Zero; count = 0; // Find the min/max over all threads and accumulate count, sum and // sum of squares minimum = NumericTraits< PixelType >::max(); maximum = NumericTraits< PixelType >::NonpositiveMin(); for ( i = 0; i < numberOfThreads; i++ ) { count += m_Count[i]; sum += m_ThreadSum[i]; sumOfSquares += m_SumOfSquares[i]; if ( m_ThreadMin[i] < minimum ) { minimum = m_ThreadMin[i]; } if ( m_ThreadMax[i] > maximum ) { maximum = m_ThreadMax[i]; } } // compute statistics mean = sum / static_cast< RealType >( count ); // unbiased estimate variance = ( sumOfSquares - ( sum * sum / static_cast< RealType >( count ) ) ) / ( static_cast< RealType >( count ) - 1 ); sigma = std::sqrt(variance); // Set the outputs this->GetMinimumOutput()->Set(minimum); this->GetMaximumOutput()->Set(maximum); this->GetMeanOutput()->Set(mean); this->GetSigmaOutput()->Set(sigma); this->GetVarianceOutput()->Set(variance); this->GetSumOutput()->Set(sum); } template< typename TInputImage, typename TMaskImage > void MaskedStatisticsImageFilter< TInputImage, TMaskImage > ::ThreadedGenerateData(const RegionType & outputRegionForThread, ThreadIdType threadId) { const SizeValueType size0 = outputRegionForThread.GetSize(0); if( size0 == 0) { return; } RealType realValue; PixelType value; RealType sum = NumericTraits< RealType >::Zero; RealType sumOfSquares = NumericTraits< RealType >::Zero; SizeValueType count = NumericTraits< SizeValueType >::Zero; PixelType min = NumericTraits< PixelType >::max(); PixelType max = NumericTraits< PixelType >::NonpositiveMin(); ImageScanlineConstIterator< TInputImage > it (this->GetInput(), outputRegionForThread); // support progress methods/callbacks const size_t numberOfLinesToProcess = outputRegionForThread.GetNumberOfPixels() / size0; ProgressReporter progress( this, threadId, numberOfLinesToProcess ); // do the work while ( !it.IsAtEnd() ) { while ( !it.IsAtEndOfLine() ) { bool isValid = true; if(m_Mask.IsNotNull()) { typename InputImageType::IndexType index = it.GetIndex(); typename InputImageType::PointType point; this->GetInput()->TransformIndexToPhysicalPoint(index, point); if (this->m_Mask->TransformPhysicalPointToIndex(point, index)) { isValid = this->m_Mask->GetPixel(index) > 0.0; }; } if (isValid) { value = it.Get(); realValue = static_cast< RealType >( value ); if ( value < min ) { min = value; } if ( value > max ) { max = value; } sum += realValue; sumOfSquares += ( realValue * realValue ); ++count; } ++it; } it.NextLine(); progress.CompletedPixel(); } m_ThreadSum[threadId] = sum; m_SumOfSquares[threadId] = sumOfSquares; m_Count[threadId] = count; m_ThreadMin[threadId] = min; m_ThreadMax[threadId] = max; } template< typename TImage, typename TMaskImage > void MaskedStatisticsImageFilter< TImage, TMaskImage > ::PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << "Minimum: " << static_cast< typename NumericTraits< PixelType >::PrintType >( this->GetMinimum() ) << std::endl; os << indent << "Maximum: " << static_cast< typename NumericTraits< PixelType >::PrintType >( this->GetMaximum() ) << std::endl; os << indent << "Sum: " << this->GetSum() << std::endl; os << indent << "Mean: " << this->GetMean() << std::endl; os << indent << "Sigma: " << this->GetSigma() << std::endl; os << indent << "Variance: " << this->GetVariance() << std::endl; } } // end namespace itk #endif diff --git a/Modules/ModelFit/include/itkMultiOutputNaryFunctorImageFilter.tpp b/Modules/ModelFit/include/itkMultiOutputNaryFunctorImageFilter.tpp index 09d16df573..6f2f9adb75 100644 --- a/Modules/ModelFit/include/itkMultiOutputNaryFunctorImageFilter.tpp +++ b/Modules/ModelFit/include/itkMultiOutputNaryFunctorImageFilter.tpp @@ -1,231 +1,233 @@ /*============================================================================ 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 __itkMultiOutputNaryFunctorImageFilter_hxx #define __itkMultiOutputNaryFunctorImageFilter_hxx #include "itkMultiOutputNaryFunctorImageFilter.h" #include "itkImageRegionIterator.h" #include "itkProgressReporter.h" namespace itk { /** * Constructor */ template< class TInputImage, class TOutputImage, class TFunction, class TMaskImage > MultiOutputNaryFunctorImageFilter< TInputImage, TOutputImage, TFunction, TMaskImage > ::MultiOutputNaryFunctorImageFilter() { + this->DynamicMultiThreadingOff(); + // This number will be incremented each time an image // is added over the two minimum required this->SetNumberOfRequiredInputs(1); this->ActualizeOutputs(); } template< class TInputImage, class TOutputImage, class TFunction, class TMaskImage > void MultiOutputNaryFunctorImageFilter< TInputImage, TOutputImage, TFunction, TMaskImage > ::ActualizeOutputs() { this->SetNumberOfRequiredOutputs(m_Functor.GetNumberOfOutputs()); for (typename Superclass::DataObjectPointerArraySizeType i = this->GetNumberOfIndexedOutputs(); i< m_Functor.GetNumberOfOutputs(); ++i) { this->SetNthOutput( i, this->MakeOutput(i) ); } while(this->GetNumberOfIndexedOutputs() > m_Functor.GetNumberOfOutputs()) { this->RemoveOutput(this->GetNumberOfIndexedOutputs()-1); } }; /** * ThreadedGenerateData Performs the pixel-wise addition */ template< class TInputImage, class TOutputImage, class TFunction, class TMaskImage > void MultiOutputNaryFunctorImageFilter< TInputImage, TOutputImage, TFunction, TMaskImage > ::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType threadId) { ProgressReporter progress( this, threadId, outputRegionForThread.GetNumberOfPixels() ); const unsigned int numberOfInputImages = static_cast< unsigned int >( this->GetNumberOfIndexedInputs() ); const unsigned int numberOfOutputImages = static_cast< unsigned int >( this->GetNumberOfIndexedOutputs() ); typedef ImageRegionConstIterator< TInputImage > ImageRegionConstIteratorType; std::vector< ImageRegionConstIteratorType * > inputItrVector; inputItrVector.reserve(numberOfInputImages); typedef ImageRegionIterator< TOutputImage > OutputImageRegionIteratorType; std::vector< OutputImageRegionIteratorType * > outputItrVector; outputItrVector.reserve(numberOfOutputImages); //check if mask image is set and generate iterator if mask is valid typedef ImageRegionConstIterator< TMaskImage > MaskImageRegionIteratorType; MaskImageRegionIteratorType* pMaskIterator = nullptr; if (m_Mask.IsNotNull()) { if (!m_Mask->GetLargestPossibleRegion().IsInside(outputRegionForThread)) { itkExceptionMacro("Mask of filter is set but does not cover region of thread. Mask region: "<< m_Mask->GetLargestPossibleRegion() <<"Thread region: "<( ProcessObject::GetInput(i) ); if ( inputPtr ) { inputItrVector.push_back( new ImageRegionConstIteratorType(inputPtr, outputRegionForThread) ); } } // go through the outputs and add iterators for non-null outputs for ( unsigned int i = 0; i < numberOfOutputImages; ++i ) { OutputImagePointer outputPtr = dynamic_cast< TOutputImage * >( ProcessObject::GetOutput(i) ); if ( outputPtr ) { outputItrVector.push_back( new OutputImageRegionIteratorType(outputPtr, outputRegionForThread) ); } } typename std::vector< ImageRegionConstIteratorType * >::iterator regionInputIterators; const typename std::vector< ImageRegionConstIteratorType * >::const_iterator regionInputItEnd = inputItrVector.end(); typename std::vector< OutputImageRegionIteratorType * >::iterator regionOutputIterators; const typename std::vector< OutputImageRegionIteratorType * >::const_iterator regionOutputItEnd = outputItrVector.end(); const unsigned int numberOfValidInputImages = inputItrVector.size(); const unsigned int numberOfValidOutputImages = outputItrVector.size(); if ( (numberOfValidInputImages != 0) && ( numberOfValidOutputImages != 0)) { try { while ( !(outputItrVector.front()->IsAtEnd()) ) { typename NaryInputArrayType::iterator arrayInIt; typename NaryOutputArrayType::iterator arrayOutIt; NaryInputArrayType naryInputArray(numberOfValidInputImages); NaryOutputArrayType naryOutputArray(numberOfValidOutputImages); bool isValid = true; if (pMaskIterator) { isValid = pMaskIterator->Get() > 0; ++(*pMaskIterator); } arrayInIt = naryInputArray.begin(); regionInputIterators = inputItrVector.begin(); typename ImageRegionConstIteratorType::IndexType currentIndex; if(regionInputIterators != regionInputItEnd) { currentIndex = ( *regionInputIterators )->GetIndex(); } while ( regionInputIterators != regionInputItEnd ) { *arrayInIt++ = ( *regionInputIterators )->Get(); ++( *( *regionInputIterators ) ); ++regionInputIterators; } if (isValid) { naryOutputArray = m_Functor(naryInputArray, currentIndex); if (numberOfValidOutputImages != naryOutputArray.size()) { itkExceptionMacro("Error. Number of valid output images do not equal number of outputs required by functor. Number of valid outputs: "<< numberOfValidOutputImages << "; needed output number:" << this->m_Functor.GetNumberOfOutputs()); } } else { for (typename NaryOutputArrayType::iterator pos = naryOutputArray.begin(); pos!= naryOutputArray.end(); ++pos) { *pos = 0.0; } } arrayOutIt = naryOutputArray.begin(); regionOutputIterators = outputItrVector.begin(); while ( regionOutputIterators != regionOutputItEnd ) { ( *regionOutputIterators )->Set(*arrayOutIt++); ++( *( *regionOutputIterators ) ); ++regionOutputIterators; } progress.CompletedPixel(); } } catch(...) { // Free memory in case of exceptions regionInputIterators = inputItrVector.begin(); while ( regionInputIterators != regionInputItEnd ) { delete ( *regionInputIterators++ ); } regionOutputIterators = outputItrVector.begin(); while ( regionOutputIterators != regionOutputItEnd ) { delete ( *regionOutputIterators++ ); } delete pMaskIterator; throw; } } // Free memory regulary regionInputIterators = inputItrVector.begin(); while ( regionInputIterators != regionInputItEnd ) { delete ( *regionInputIterators++ ); } regionOutputIterators = outputItrVector.begin(); while ( regionOutputIterators != regionOutputItEnd ) { delete ( *regionOutputIterators++ ); } delete pMaskIterator; } } // end namespace itk #endif