diff --git a/Modules/ImageDenoising/itkLocalVariationImageFilter.txx b/Modules/ImageDenoising/itkLocalVariationImageFilter.txx index 2cfa5e98d6..3048232a6c 100644 --- a/Modules/ImageDenoising/itkLocalVariationImageFilter.txx +++ b/Modules/ImageDenoising/itkLocalVariationImageFilter.txx @@ -1,193 +1,193 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef _itkLocalVariationImageFilter_txx #define _itkLocalVariationImageFilter_txx #include "itkLocalVariationImageFilter.h" #include "itkConstShapedNeighborhoodIterator.h" #include "itkImageRegionIterator.h" #include "itkNeighborhoodAlgorithm.h" #include "itkNeighborhoodInnerProduct.h" #include "itkOffset.h" #include "itkProgressReporter.h" #include "itkVectorImage.h" #include "itkZeroFluxNeumannBoundaryCondition.h" #include #include namespace itk { template LocalVariationImageFilter::LocalVariationImageFilter() { } template void LocalVariationImageFilter::GenerateInputRequestedRegion() throw( InvalidRequestedRegionError) { // call the superclass' implementation of this method Superclass::GenerateInputRequestedRegion(); // get pointers to the input and output typename Superclass::InputImagePointer inputPtr = const_cast(this->GetInput()); typename Superclass::OutputImagePointer outputPtr = this->GetOutput(); if (!inputPtr || !outputPtr) { return; } // get a copy of the input requested region (should equal the output // requested region) typename TInputImage::RegionType inputRequestedRegion; inputRequestedRegion = inputPtr->GetRequestedRegion(); // pad the input requested region by 1 inputRequestedRegion.PadByRadius(1); // crop the input requested region at the input's largest possible region if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion())) { inputPtr->SetRequestedRegion(inputRequestedRegion); return; } else { // Couldn't crop the region (requested region is outside the largest // possible region). Throw an exception. // store what we tried to request (prior to trying to crop) inputPtr->SetRequestedRegion(inputRequestedRegion); // build an exception InvalidRequestedRegionError e(__FILE__, __LINE__); e.SetLocation(ITK_LOCATION); e.SetDescription("Requested region outside possible region."); e.SetDataObject(inputPtr); throw e; } } template <> double SquaredEuclideanMetric>::Calc(itk::VariableLengthVector p) { return p.GetSquaredNorm(); } template <> double SquaredEuclideanMetric>::Calc(itk::VariableLengthVector p) { return p.GetSquaredNorm(); } template double SquaredEuclideanMetric::Calc(TPixelType p) { return p * p; } template void LocalVariationImageFilter::ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId) { // Allocate output typename OutputImageType::Pointer output = this->GetOutput(); typename InputImageType::ConstPointer input = this->GetInput(); itk::Size size; - for (int i = 0; i < InputImageDimension; i++) + for (unsigned int i = 0; i < InputImageDimension; i++) size[i] = 1; // Find the data-set boundary "faces" NeighborhoodAlgorithm::ImageBoundaryFacesCalculator bC; typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator::FaceListType faceList = bC(input, outputRegionForThread, size); // support progress methods/callbacks ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); ZeroFluxNeumannBoundaryCondition nbc; std::vector pixels; // Process each of the boundary faces. These are N-d regions which border // the edge of the buffer. for (auto fit = faceList.begin(); fit != faceList.end(); ++fit) { // iterators over output and input ImageRegionIterator output_image_it(output, *fit); ImageRegionConstIterator input_image_it(input.GetPointer(), *fit); // neighborhood iterator for input image ConstShapedNeighborhoodIterator input_image_neighbors_it(size, input, *fit); typename ConstShapedNeighborhoodIterator::OffsetType offset; input_image_neighbors_it.OverrideBoundaryCondition(&nbc); input_image_neighbors_it.ClearActiveList(); - for (int i = 0; i < InputImageDimension; i++) + for (unsigned int i = 0; i < InputImageDimension; i++) { offset.Fill(0); offset[i] = -1; input_image_neighbors_it.ActivateOffset(offset); offset[i] = 1; input_image_neighbors_it.ActivateOffset(offset); } input_image_neighbors_it.GoToBegin(); // const unsigned int neighborhoodSize = InputImageDimension*2; while (!input_image_neighbors_it.IsAtEnd()) { // collect all the pixels in the neighborhood, note that we use // GetPixel on the NeighborhoodIterator to honor the boundary conditions typename OutputImageType::PixelType locVariation = 0; typename ConstShapedNeighborhoodIterator::ConstIterator input_neighbors_it; for (input_neighbors_it = input_image_neighbors_it.Begin(); !input_neighbors_it.IsAtEnd(); input_neighbors_it++) { typename TInputImage::PixelType diffVec = input_neighbors_it.Get() - input_image_it.Get(); locVariation += SquaredEuclideanMetric::Calc(diffVec); } locVariation = sqrt(locVariation + 0.0001); output_image_it.Set(locVariation); // update iterators ++input_image_neighbors_it; ++output_image_it; ++input_image_it; // report progress progress.CompletedPixel(); } } } /** * Standard "PrintSelf" method */ template void LocalVariationImageFilter::PrintSelf(std::ostream &os, Indent indent) const { Superclass::PrintSelf(os, indent); } } // end namespace itk #endif //_itkLocalVariationImageFilter_txx diff --git a/Modules/ImageDenoising/itkTotalVariationSingleIterationImageFilter.txx b/Modules/ImageDenoising/itkTotalVariationSingleIterationImageFilter.txx index 74e988d84e..13bd039e27 100644 --- a/Modules/ImageDenoising/itkTotalVariationSingleIterationImageFilter.txx +++ b/Modules/ImageDenoising/itkTotalVariationSingleIterationImageFilter.txx @@ -1,256 +1,256 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef _itkTotalVariationSingleIterationImageFilter_txx #define _itkTotalVariationSingleIterationImageFilter_txx #include "itkTotalVariationSingleIterationImageFilter.h" // itk includes #include "itkConstShapedNeighborhoodIterator.h" #include "itkImageRegionIterator.h" #include "itkLocalVariationImageFilter.h" #include "itkNeighborhoodAlgorithm.h" #include "itkNeighborhoodInnerProduct.h" #include "itkOffset.h" #include "itkProgressReporter.h" #include "itkZeroFluxNeumannBoundaryCondition.h" // other includes #include #include namespace itk { /** * constructor */ template TotalVariationSingleIterationImageFilter::TotalVariationSingleIterationImageFilter() { m_Lambda = 1.0; m_LocalVariation = LocalVariationImageType::New(); } /** * generate requested region */ template void TotalVariationSingleIterationImageFilter::GenerateInputRequestedRegion() throw( InvalidRequestedRegionError) { // call the superclass' implementation of this method Superclass::GenerateInputRequestedRegion(); // get pointers to the input and output typename Superclass::InputImagePointer inputPtr = const_cast(this->GetInput()); typename Superclass::OutputImagePointer outputPtr = this->GetOutput(); if (!inputPtr || !outputPtr) { return; } // get a copy of the input requested region (should equal the output // requested region) typename TInputImage::RegionType inputRequestedRegion; inputRequestedRegion = inputPtr->GetRequestedRegion(); // pad the input requested region by 1 inputRequestedRegion.PadByRadius(1); // crop the input requested region at the input's largest possible region if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion())) { inputPtr->SetRequestedRegion(inputRequestedRegion); return; } else { // Couldn't crop the region (requested region is outside the largest // possible region). Throw an exception. // store what we tried to request (prior to trying to crop) inputPtr->SetRequestedRegion(inputRequestedRegion); // build an exception InvalidRequestedRegionError e(__FILE__, __LINE__); e.SetLocation(ITK_LOCATION); e.SetDescription("Requested region outside possible region."); e.SetDataObject(inputPtr); throw e; } } /** * generate output */ template void TotalVariationSingleIterationImageFilter::ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId) { typename OutputImageType::Pointer output = this->GetOutput(); typename InputImageType::ConstPointer input = this->GetInput(); // Find the data-set boundary "faces" itk::Size size; - for (int i = 0; i < InputImageDimension; i++) + for (unsigned int i = 0; i < InputImageDimension; i++) size[i] = 1; NeighborhoodAlgorithm::ImageBoundaryFacesCalculator bC; typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator::FaceListType faceList = bC(input, outputRegionForThread, size); NeighborhoodAlgorithm::ImageBoundaryFacesCalculator lv_bC; typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator::FaceListType lv_faceList = lv_bC(m_LocalVariation, outputRegionForThread, size); // support progress methods/callbacks ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); ZeroFluxNeumannBoundaryCondition nbc; ZeroFluxNeumannBoundaryCondition lv_nbc; std::vector ws; std::vector hs; auto lv_fit = lv_faceList.begin(); // Process each of the boundary faces. These are N-d regions which border // the edge of the buffer. for (auto fit = faceList.begin(); fit != faceList.end(); ++fit) { // iterators over output, input, original and local variation image ImageRegionIterator output_image_it = ImageRegionIterator(output, *fit); ImageRegionConstIterator input_image_it = ImageRegionConstIterator(input, *fit); ImageRegionConstIterator orig_image_it = ImageRegionConstIterator(m_OriginalImage, *fit); ImageRegionConstIterator loc_var_image_it = ImageRegionConstIterator(m_LocalVariation, *fit); // neighborhood in input image ConstShapedNeighborhoodIterator input_image_neighbors_it(size, input, *fit); typename ConstShapedNeighborhoodIterator::OffsetType offset; input_image_neighbors_it.OverrideBoundaryCondition(&nbc); input_image_neighbors_it.ClearActiveList(); - for (int i = 0; i < InputImageDimension; i++) + for (unsigned int i = 0; i < InputImageDimension; i++) { offset.Fill(0); offset[i] = -1; input_image_neighbors_it.ActivateOffset(offset); offset[i] = 1; input_image_neighbors_it.ActivateOffset(offset); } input_image_neighbors_it.GoToBegin(); // neighborhood in local variation image ConstShapedNeighborhoodIterator loc_var_image_neighbors_it( size, m_LocalVariation, *lv_fit); loc_var_image_neighbors_it.OverrideBoundaryCondition(&lv_nbc); loc_var_image_neighbors_it.ClearActiveList(); - for (int i = 0; i < InputImageDimension; i++) + for (unsigned int i = 0; i < InputImageDimension; i++) { offset.Fill(0); offset[i] = -1; loc_var_image_neighbors_it.ActivateOffset(offset); offset[i] = 1; loc_var_image_neighbors_it.ActivateOffset(offset); } loc_var_image_neighbors_it.GoToBegin(); const unsigned int neighborhoodSize = InputImageDimension * 2; ws.resize(neighborhoodSize); while (!output_image_it.IsAtEnd()) { // 1 / ||nabla_alpha(u)||_a double locvar_alpha_inv = 1.0 / loc_var_image_it.Get(); // compute w_alphabetas int count = 0; double wsum = 0; typename ConstShapedNeighborhoodIterator::ConstIterator loc_var_neighbors_it; for (loc_var_neighbors_it = loc_var_image_neighbors_it.Begin(); !loc_var_neighbors_it.IsAtEnd(); loc_var_neighbors_it++) { // w_alphabeta(u) = // 1 / ||nabla_alpha(u)||_a + 1 / ||nabla_beta(u)||_a ws[count] = locvar_alpha_inv + (1.0 / (double)loc_var_neighbors_it.Get()); wsum += ws[count++]; } // h_alphaalpha * u_alpha^zero typename OutputImageType::PixelType res = static_cast( ((typename OutputImageType::PixelType)orig_image_it.Get()) * (m_Lambda / (m_Lambda + wsum))); // add the different h_alphabeta * u_beta count = 0; typename ConstShapedNeighborhoodIterator::ConstIterator input_neighbors_it; for (input_neighbors_it = input_image_neighbors_it.Begin(); !input_neighbors_it.IsAtEnd(); input_neighbors_it++) { res += input_neighbors_it.Get() * (ws[count++] / (m_Lambda + wsum)); } // set output result output_image_it.Set(res); // increment iterators ++output_image_it; ++input_image_it; ++orig_image_it; ++loc_var_image_it; ++input_image_neighbors_it; ++loc_var_image_neighbors_it; // report progress progress.CompletedPixel(); } ++lv_fit; } } /** * first calculate local variation in the image */ template void TotalVariationSingleIterationImageFilter::BeforeThreadedGenerateData() { typedef typename itk::LocalVariationImageFilter FilterType; typename FilterType::Pointer filter = FilterType::New(); filter->SetInput(this->GetInput(0)); filter->SetNumberOfThreads(this->GetNumberOfThreads()); filter->Update(); this->m_LocalVariation = filter->GetOutput(); } /** * Standard "PrintSelf" method */ template void TotalVariationSingleIterationImageFilter::PrintSelf(std::ostream &os, Indent indent) const { Superclass::PrintSelf(os, indent); } } // end namespace itk #endif