diff --git a/Modules/Classification/CLUtilities/include/itkLocalStatisticFilter.hxx b/Modules/Classification/CLUtilities/include/itkLocalStatisticFilter.hxx index a8fabc92fe..d6f18f07e5 100644 --- a/Modules/Classification/CLUtilities/include/itkLocalStatisticFilter.hxx +++ b/Modules/Classification/CLUtilities/include/itkLocalStatisticFilter.hxx @@ -1,113 +1,114 @@ #ifndef itkLocalStatisticFilter_cpp #define itkLocalStatisticFilter_cpp #include #include #include #include #include "itkMinimumMaximumImageCalculator.h" #include template< class TInputImageType, class TOuputImageType> itk::LocalStatisticFilter::LocalStatisticFilter(): m_Size(5), m_Bins(5) { + this->DynamicMultiThreadingOff(); this->SetNumberOfRequiredOutputs(m_Bins); this->SetNumberOfRequiredInputs(0); for (int i = 0; i < m_Bins; ++i) { this->SetNthOutput( i, this->MakeOutput(i) ); } } template< class TInputImageType, class TOuputImageType> void itk::LocalStatisticFilter::BeforeThreadedGenerateData() { InputImagePointer input = this->GetInput(0); for (int i = 0; i < m_Bins; ++i) { CreateOutputImage(input, this->GetOutput(i)); } } template< class TInputImageType, class TOuputImageType> void itk::LocalStatisticFilter::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType /*threadId*/) { typedef itk::ImageRegionIterator IteratorType; typedef itk::ConstNeighborhoodIterator ConstIteratorType; typename TInputImageType::SizeType size; size.Fill(m_Size); InputImagePointer input = this->GetInput(0); if (TInputImageType::ImageDimension == 3) { size[2] = 0; } // MITK_INFO << "Creating output iterator"; std::vector iterVector; for (int i = 0; i < m_Bins; ++i) { IteratorType iter(this->GetOutput(i), outputRegionForThread); iterVector.push_back(iter); } ConstIteratorType inputIter(size, input, outputRegionForThread); while (!inputIter.IsAtEnd()) { for (int i = 0; i < m_Bins; ++i) { iterVector[i].Set(0); } double min = std::numeric_limits::max(); double max = std::numeric_limits::lowest(); double mean = 0; double std = 0; for (unsigned int i = 0; i < inputIter.Size(); ++i) { double value = inputIter.GetPixel(i); min = std::min(min, value); max = std::max(max, value); mean += value / inputIter.Size(); std += (value*value) / inputIter.Size(); } iterVector[0].Value() = min; iterVector[1].Value() = max; iterVector[2].Value() = mean; iterVector[3].Value() = std::sqrt(std - mean*mean); iterVector[4].Value() = max-min; for (int i = 0; i < m_Bins; ++i) { ++(iterVector[i]); } ++inputIter; } } template< class TInputImageType, class TOuputImageType> itk::ProcessObject::DataObjectPointer itk::LocalStatisticFilter::MakeOutput(itk::ProcessObject::DataObjectPointerArraySizeType /*idx*/) { itk::ProcessObject::DataObjectPointer output; output = ( TOuputImageType::New() ).GetPointer(); return output; } template< class TInputImageType, class TOuputImageType> void itk::LocalStatisticFilter::CreateOutputImage(InputImagePointer input, OutputImagePointer output) { output->SetRegions(input->GetLargestPossibleRegion()); output->Allocate(); } #endif //itkLocalStatisticFilter_cpp diff --git a/Modules/Classification/CLUtilities/include/itkMultiHistogramFilter.cpp b/Modules/Classification/CLUtilities/include/itkMultiHistogramFilter.cpp index b18a65904c..2eabf31065 100644 --- a/Modules/Classification/CLUtilities/include/itkMultiHistogramFilter.cpp +++ b/Modules/Classification/CLUtilities/include/itkMultiHistogramFilter.cpp @@ -1,124 +1,125 @@ /*============================================================================ 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 itkMultiHistogramFilter_cpp #define itkMultiHistogramFilter_cpp #include #include #include #include #include "itkMinimumMaximumImageCalculator.h" template< class TInputImageType, class TOuputImageType> itk::MultiHistogramFilter::MultiHistogramFilter(): m_Delta(0.6), m_Offset(-3.0), m_Bins(11), m_Size(5), m_UseImageIntensityRange(false) { + this->DynamicMultiThreadingOff(); this->SetNumberOfRequiredOutputs(m_Bins); this->SetNumberOfRequiredInputs(0); for (int i = 0; i < m_Bins; ++i) { this->SetNthOutput( i, this->MakeOutput(i) ); } } template< class TInputImageType, class TOuputImageType> void itk::MultiHistogramFilter::BeforeThreadedGenerateData() { typedef itk::MinimumMaximumImageCalculator ImageCalculatorFilterType; if (m_UseImageIntensityRange) { typename ImageCalculatorFilterType::Pointer imageCalculatorFilter = ImageCalculatorFilterType::New(); imageCalculatorFilter->SetImage(this->GetInput(0)); imageCalculatorFilter->Compute(); m_Offset = imageCalculatorFilter->GetMinimum(); m_Delta = 1.0*(imageCalculatorFilter->GetMaximum() - imageCalculatorFilter->GetMinimum()) / (1.0*m_Bins); } InputImagePointer input = this->GetInput(0); for (int i = 0; i < m_Bins; ++i) { CreateOutputImage(input, this->GetOutput(i)); } } template< class TInputImageType, class TOuputImageType> void itk::MultiHistogramFilter::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType /*threadId*/) { double offset = m_Offset;// -3.0; double delta = m_Delta;// 0.6; typedef itk::ImageRegionIterator IteratorType; typedef itk::ConstNeighborhoodIterator ConstIteratorType; typename TInputImageType::SizeType size; size.Fill(m_Size); InputImagePointer input = this->GetInput(0); // MITK_INFO << "Creating output iterator"; std::vector iterVector; for (int i = 0; i < m_Bins; ++i) { IteratorType iter(this->GetOutput(i), outputRegionForThread); iterVector.push_back(iter); } ConstIteratorType inputIter(size, input, outputRegionForThread); while (!inputIter.IsAtEnd()) { for (int i = 0; i < m_Bins; ++i) { iterVector[i].Set(0); } for (unsigned int i = 0; i < inputIter.Size(); ++i) { double value = inputIter.GetPixel(i); value -= offset; value /= delta; auto pos = (int)(value); pos = std::max(0, std::min(m_Bins-1, pos)); iterVector[pos].Value() += 1;// (iterVector[pos].GetCenterPixel() + 1); } for (int i = 0; i < m_Bins; ++i) { ++(iterVector[i]); } ++inputIter; } } template< class TInputImageType, class TOuputImageType> itk::ProcessObject::DataObjectPointer itk::MultiHistogramFilter::MakeOutput(itk::ProcessObject::DataObjectPointerArraySizeType /*idx*/) { itk::ProcessObject::DataObjectPointer output; output = ( TOuputImageType::New() ).GetPointer(); return output; } template< class TInputImageType, class TOuputImageType> void itk::MultiHistogramFilter::CreateOutputImage(InputImagePointer input, OutputImagePointer output) { output->SetRegions(input->GetLargestPossibleRegion()); output->Allocate(); } #endif //itkMultiHistogramFilter_cpp diff --git a/Modules/Classification/CLUtilities/include/itkNeighborhoodFunctorImageFilter.h b/Modules/Classification/CLUtilities/include/itkNeighborhoodFunctorImageFilter.h index 6942ce74c4..229a132ebd 100644 --- a/Modules/Classification/CLUtilities/include/itkNeighborhoodFunctorImageFilter.h +++ b/Modules/Classification/CLUtilities/include/itkNeighborhoodFunctorImageFilter.h @@ -1,152 +1,153 @@ /*============================================================================ 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 itkNeighborhoodFunctorImageFilter_h #define itkNeighborhoodFunctorImageFilter_h #include "itkImageToImageFilter.h" #include "itkZeroFluxNeumannBoundaryCondition.h" #include "itkConstNeighborhoodIterator.h" #include "itkImage.h" #include #include #include "itkHistogram.h" namespace itk { template class NeighborhoodFunctorImageFilter : public ImageToImageFilter< TInputImageType, TFeatureImageType> { public: typedef NeighborhoodFunctorImageFilter Self; typedef ImageToImageFilter Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; itkNewMacro(Self); itkTypeMacro(NeighborhoodFunctorImageFilter, ImageToImageFilter); /** Extract some information from the image types. Dimensionality * of the two images is assumed to be the same. */ itkStaticConstMacro(ImageDimension, unsigned int, TFeatureImageType::ImageDimension); itkStaticConstMacro(InputImageDimension, unsigned int, TInputImageType::ImageDimension); typedef TInputImageType InputImageType; typedef typename TInputImageType::PixelType InputImagePixelType; typedef itk::Image MaskImageType; typedef typename MaskImageType::PixelType MaskImagePixelType; typedef TFeatureImageType FeatureImageType; typedef typename FeatureImageType::PixelType FeaturePixelType; typedef itk::Size SizeType; /** Typedef for generic boundary condition pointer. */ typedef ImageBoundaryCondition< InputImageType > * ImageBoundaryConditionPointerType; /** Typedef for the default boundary condition */ typedef ZeroFluxNeumannBoundaryCondition< InputImageType > DefaultBoundaryCondition; /** Superclass typedefs. */ typedef typename Superclass::OutputImageRegionType OutputImageRegionType; typedef Neighborhood< InputImagePixelType, InputImageDimension > NeighborhoodType; /** Allows a user to override the internal boundary condition. Care should be * be taken to ensure that the overriding boundary condition is a persistent * object during the time it is referenced. The overriding condition * can be of a different type than the default type as long as it is * a subclass of ImageBoundaryCondition. */ void OverrideBoundaryCondition(const ImageBoundaryConditionPointerType i) { m_BoundsCondition = i; } /** Get the boundary condition specified */ ImageBoundaryConditionPointerType GetBoundaryCondition() { return m_BoundsCondition; } void SetNeighborhoodSize(SizeType size){m_Size = size;} void SetNeighborhoodSize(unsigned int size){m_Size.Fill(size);} SizeType GetNeighborhoodSize(){return m_Size;} void SetMask(const typename MaskImageType::Pointer & ptr){m_MaskImage = ptr;} const FunctorType & GetFunctorReference() const { return m_Functor; } FunctorType & GetFunctorReference() { return m_Functor; } void SetFunctor(const FunctorType & func) { m_Functor = func; } protected: NeighborhoodFunctorImageFilter() { + this->DynamicMultiThreadingOff(); m_Size.Fill(0); m_MaskImage = nullptr; m_BoundsCondition = static_cast< ImageBoundaryConditionPointerType >( &m_DefaultBoundaryCondition ); this->SetNumberOfIndexedOutputs(FunctorType::OutputCount); } ~NeighborhoodFunctorImageFilter() override{} void BeforeThreadedGenerateData() override; void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType threadId) override; /** NeighborhoodFunctorImageFilter needs a larger input requested * region than the output requested region. As such, * NeighborhoodOperatorImageFilter needs to provide an implementation for * GenerateInputRequestedRegion() in order to inform the pipeline * execution model. * * \sa ProcessObject::GenerateInputRequestedRegion() */ void GenerateInputRequestedRegion() override; private: NeighborhoodFunctorImageFilter(const Self &); // purposely not implemented void operator=(const Self &); // purposely not implemented /** Pointer to a persistent boundary condition object used * for the image iterator. */ ImageBoundaryConditionPointerType m_BoundsCondition; /** Default boundary condition */ DefaultBoundaryCondition m_DefaultBoundaryCondition; /** Internal operator used to filter the image. */ FunctorType m_Functor; itk::Size m_Size; typename MaskImageType::Pointer m_MaskImage; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "../src/Features/itkNeighborhoodFunctorImageFilter.cpp" #endif #endif // itkFeatureImageFilter_h diff --git a/Modules/Classification/CLUtilities/src/Features/itkLineHistogramBasedMassImageFilter.cpp b/Modules/Classification/CLUtilities/src/Features/itkLineHistogramBasedMassImageFilter.cpp index 1d5281bd08..626150db4b 100644 --- a/Modules/Classification/CLUtilities/src/Features/itkLineHistogramBasedMassImageFilter.cpp +++ b/Modules/Classification/CLUtilities/src/Features/itkLineHistogramBasedMassImageFilter.cpp @@ -1,255 +1,256 @@ /*============================================================================ 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 ITKLINEHISTOGRAMBASEDMASSIMAGEFILTER_CPP #define ITKLINEHISTOGRAMBASEDMASSIMAGEFILTER_CPP #include #include #include #include #include #include template< class TInputImageType, class TOutputImageType, class TMaskImageType> void itk::LineHistogramBasedMassImageFilter ::BeforeThreadedGenerateData() { if(m_ImageMask.IsNull()) { itkExceptionMacro("No binary mask image provided.") } if(m_BinaryContour.IsNull()){ typename TMaskImageType::RegionType region = m_ImageMask->GetLargestPossibleRegion(); unsigned int xdim = region.GetSize(0); unsigned int ydim = region.GetSize(1); unsigned int zdim = region.GetSize(2); // ensure border pixels are zero so that an closed contour is created typename TMaskImageType::Pointer mask_copy = TMaskImageType::New(); mask_copy->SetSpacing(m_ImageMask->GetSpacing()); mask_copy->SetDirection(m_ImageMask->GetDirection()); mask_copy->SetOrigin(m_ImageMask->GetOrigin()); mask_copy->SetRegions(region); mask_copy->Allocate(); mask_copy->FillBuffer(0); itk::ImageRegionIteratorWithIndex oit(mask_copy, mask_copy->GetLargestPossibleRegion()); itk::ImageRegionConstIteratorWithIndex mit(m_ImageMask, m_ImageMask->GetLargestPossibleRegion()); while(!mit.IsAtEnd()) { if(mit.Value() != 0 //is under the mask && mit.GetIndex()[0] != 0 //is not at the border && mit.GetIndex()[1] != 0 && mit.GetIndex()[2] != 0 && mit.GetIndex()[0] != xdim-1 && mit.GetIndex()[1] != ydim-1 && mit.GetIndex()[2] != zdim-1) { oit.Set(1); } ++mit; ++oit; } typedef itk::BinaryContourImageFilter BinaryContourImagefilterType; typename BinaryContourImagefilterType::Pointer filter = BinaryContourImagefilterType::New(); filter->SetInput(mask_copy); filter->SetBackgroundValue(0); filter->SetForegroundValue(1); filter->SetFullyConnected(true); filter->Update(); m_BinaryContour = filter->GetOutput(); mitk::Image::Pointer outimg; mitk::CastToMitkImage(m_BinaryContour,outimg); } if(m_BinaryContour.IsNull()) { itkExceptionMacro("No binary contour image provided.") } m_CenterOfMask = GetCenterOfMass(m_ImageMask); if(m_CenterOfMask.is_zero()) { itkExceptionMacro("Center of mass is corrupt.") } } template< class TInputImageType, class TOutputImageType, class TMaskImageType> vnl_vector itk::LineHistogramBasedMassImageFilter ::GetCenterOfMass(const TMaskImageType * maskImage) { mitk::Image::Pointer img; mitk::CastToMitkImage(maskImage, img); typedef itk::ImageRegionConstIterator< TMaskImageType > IteratorType; IteratorType iter( maskImage, maskImage->GetLargestPossibleRegion() ); iter.GoToBegin(); vnl_vector_fixed mean_index; mean_index.fill(0); unsigned int count = 0; while ( !iter.IsAtEnd() ) { if ( iter.Get() != 0 ) { mitk::Point3D current_index_pos; img->GetGeometry()->IndexToWorld(iter.GetIndex(),current_index_pos); mean_index += current_index_pos.GetVnlVector(); count++; } ++iter; } mean_index /= count; return mean_index.as_ref(); } template< class TInputImageType, class TOutputImageType, class TMaskImageType> void itk::LineHistogramBasedMassImageFilter ::ThreadedGenerateData(const typename Superclass::OutputImageRegionType &outputRegionForThread, ThreadIdType /*threadId*/) { TOutputImageType * outimage = this->GetOutput(); itk::ImageRegionConstIteratorWithIndex cit(m_BinaryContour, outputRegionForThread); itk::ImageRegionConstIteratorWithIndex iit(this->GetInput(), outputRegionForThread); itk::ImageRegionIteratorWithIndex oit(outimage,outputRegionForThread); typedef typename itk::ImageRegionIteratorWithIndex::IndexType IndexType; std::vector target_world_indices; while(!cit.IsAtEnd()) { if(cit.Value() != 0 ) target_world_indices.push_back(cit.GetIndex()); ++cit; } mitk::Image::Pointer image; mitk::CastToMitkImage(this->GetInput(), image); mitk::BaseGeometry * transform = image->GetGeometry(); while(!target_world_indices.empty()) { vnl_vector u_v(3,1), x_v(3,1), p_v(3,1); double line_histo_area = 0; // w->i skull point mitk::Point3D skull_point, tmp_point; image->GetGeometry()->IndexToWorld(target_world_indices.back(),skull_point); // project centerpoint to slice // x = p + s*u u_v = skull_point.GetVnlVector() - m_CenterOfMask; u_v.normalize(); p_v = m_CenterOfMask; //step width double s_s = 0.1; // i->w center point mitk::Point3D center_point; center_point[0] = m_CenterOfMask[0]; center_point[1] = m_CenterOfMask[1]; center_point[2] = m_CenterOfMask[2]; IndexType tmp_index; transform->WorldToIndex(center_point,tmp_index); oit.SetIndex(tmp_index); std::vector under_line_indices; while(true) { // store current_index under_line_indices.push_back(tmp_index); // set histo value iit.SetIndex(tmp_index); line_histo_area += iit.Value(); // break in end reached if(tmp_index == target_world_indices.back()) { break; } // get next voxel index under the line while(tmp_index == oit.GetIndex()) // if new voxel index reached brake { x_v = p_v + s_s * u_v; // next tmp_point.GetVnlVector().set(x_v.data_block()); transform->WorldToIndex(tmp_point, tmp_index); s_s+=0.1; } oit.SetIndex(tmp_index); } while (!under_line_indices.empty()) { IndexType current_index = under_line_indices.back(); under_line_indices.pop_back(); oit.SetIndex(current_index); oit.Set(line_histo_area / (s_s * u_v).magnitude() ); } target_world_indices.pop_back(); } } template< class TInputImageType, class TOutputImageType, class TMaskImageType> void itk::LineHistogramBasedMassImageFilter::SetImageMask(TMaskImageType * maskimage) { this->m_ImageMask = maskimage; } template< class TInputImageType, class TOutputImageType, class TMaskImageType> void itk::LineHistogramBasedMassImageFilter::SetBinaryContour(TMaskImageType * binarycontour) { this->m_BinaryContour = binarycontour; } template< class TInputImageType, class TOutputImageType, class TMaskImageType> itk::LineHistogramBasedMassImageFilter::LineHistogramBasedMassImageFilter() { + this->DynamicMultiThreadingOff(); this->SetNumberOfIndexedOutputs(1); this->SetNumberOfIndexedInputs(1); m_CenterOfMask.fill(0); } template< class TInputImageType, class TOutputImageType, class TMaskImageType> itk::LineHistogramBasedMassImageFilter::~LineHistogramBasedMassImageFilter() { } #endif