text: /*=========================================================================
   *
   *  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 itkHessianRecursiveGaussianImageFilter_hxx
  #define itkHessianRecursiveGaussianImageFilter_hxx
  
  #include "itkHessianRecursiveGaussianImageFilter.h"
  #include "itkImageRegionIteratorWithIndex.h"
  #include "itkProgressAccumulator.h"
  
  namespace itk
  {
  /**
   * Constructor
   */
  template< typename TInputImage, typename TOutputImage >
  HessianRecursiveGaussianImageFilter< TInputImage, TOutputImage >
  ::HessianRecursiveGaussianImageFilter()
  {
    m_NormalizeAcrossScale = false;
  
    // note: this is not constant to suppress a warning
    unsigned int numberOfSmoothingFilters = NumberOfSmoothingFilters;
  
    for ( unsigned int i = 0; i < numberOfSmoothingFilters; i++ )
      {
      GaussianFilterPointer filter = GaussianFilterType::New();
      filter->SetOrder(GaussianFilterType::ZeroOrder);
      filter->SetNormalizeAcrossScale(m_NormalizeAcrossScale);
      filter->InPlaceOn();
      filter->ReleaseDataFlagOn();
      m_SmoothingFilters.push_back(filter);
      }
  
    m_DerivativeFilterA = DerivativeFilterAType::New();
    m_DerivativeFilterB = DerivativeFilterBType::New();
  
    m_DerivativeFilterA->SetOrder(DerivativeFilterAType::FirstOrder);
    m_DerivativeFilterA->SetNormalizeAcrossScale(m_NormalizeAcrossScale);
  
    m_DerivativeFilterB->SetOrder(DerivativeFilterBType::FirstOrder);
    m_DerivativeFilterB->SetNormalizeAcrossScale(m_NormalizeAcrossScale);
  
    m_DerivativeFilterA->SetInput( this->GetInput() );
    m_DerivativeFilterB->SetInput( m_DerivativeFilterA->GetOutput() );
  
    m_DerivativeFilterA->InPlaceOff();
    m_DerivativeFilterA->ReleaseDataFlagOff(); // output may be used
                                               // more then once
  
    m_DerivativeFilterB->InPlaceOn();
    m_DerivativeFilterB->ReleaseDataFlagOn(); // output is only used once
  
    // Deal with the 2D case.
    if(numberOfSmoothingFilters > 0)
      {
      m_SmoothingFilters[0]->SetInput( m_DerivativeFilterB->GetOutput() );
      }
    // connect up smoothing filter chain if necessary
    for ( unsigned int i = 1; i < numberOfSmoothingFilters; i++ )
      {
      m_SmoothingFilters[i]->SetInput(m_SmoothingFilters[i - 1]->GetOutput() );
      }
  
    m_ImageAdaptor = OutputImageAdaptorType::New();
  
    this->SetSigma(1.0);
  }
  
  /**
   * Set value of Sigma
   */
  template< typename TInputImage, typename TOutputImage >
  void
  HessianRecursiveGaussianImageFilter< TInputImage, TOutputImage >
  ::SetSigma(RealType sigma)
  {
    unsigned int numberOfSmoothingFilters = NumberOfSmoothingFilters;
  
    for ( unsigned int i = 0; i < numberOfSmoothingFilters; i++ )
      {
      m_SmoothingFilters[i]->SetSigma(sigma);
      }
    m_DerivativeFilterA->SetSigma(sigma);
    m_DerivativeFilterB->SetSigma(sigma);
  
    this->Modified();
  }
  
  /**
   * Get value of Sigma
   */
  template< typename TInputImage, typename TOutputImage >
  typename HessianRecursiveGaussianImageFilter< TInputImage, TOutputImage >
  ::RealType
   HessianRecursiveGaussianImageFilter< TInputImage, TOutputImage >
  ::GetSigma() const
  {
    return m_DerivativeFilterA->GetSigma();
  }
  
  /**
   * Set Normalize Across Scale Space
   */
  template< typename TInputImage, typename TOutputImage >
  void
  HessianRecursiveGaussianImageFilter< TInputImage, TOutputImage >
  ::SetNormalizeAcrossScale(bool normalize)
  {
    m_NormalizeAcrossScale = normalize;
  
    // No normalization across scale is needed for the zero-order
    // gaussian filters. Only the derivatives need normalization.
    m_DerivativeFilterA->SetNormalizeAcrossScale(normalize);
    m_DerivativeFilterB->SetNormalizeAcrossScale(normalize);
  
    this->Modified();
  }
  
  //
  //
  //
  template< typename TInputImage, typename TOutputImage >
  void
  HessianRecursiveGaussianImageFilter< TInputImage, TOutputImage >
  ::GenerateInputRequestedRegion()
  {
    // call the superclass' implementation of this method. this should
    // copy the output requested region to the input requested region
    Superclass::GenerateInputRequestedRegion();
  
    // This filter needs all of the input
    typename HessianRecursiveGaussianImageFilter< TInputImage,
                                                  TOutputImage >::InputImagePointer image =
      const_cast< InputImageType * >( this->GetInput() );
    if ( image )
      {
      image->SetRequestedRegion( this->GetInput()->GetLargestPossibleRegion() );
      }
  }
  
  //
  //
  //
  template< typename TInputImage, typename TOutputImage >
  void
  HessianRecursiveGaussianImageFilter< TInputImage, TOutputImage >
  ::EnlargeOutputRequestedRegion(DataObject *output)
  {
    auto * out = dynamic_cast< TOutputImage * >( output );
  
    if ( out )
      {
      out->SetRequestedRegion( out->GetLargestPossibleRegion() );
      }
  }
  
  /**
   * Compute filter for Gaussian kernel
   */
  template< typename TInputImage, typename TOutputImage >
  void
  HessianRecursiveGaussianImageFilter< TInputImage, TOutputImage >
  ::GenerateData()
  {
    itkDebugMacro(<< "HessianRecursiveGaussianImageFilter generating data ");
  
    // Create a process accumulator for tracking the progress of this
    // minipipeline
    ProgressAccumulator::Pointer progress = ProgressAccumulator::New();
    progress->SetMiniPipelineFilter(this);
  
    // Compute the contribution of each filter to the total progress.
    // This computed weight is an upper bound that assumes that all filters are
    // run for each iteration.
    // However, it is not entirely accurate because it does not take into account
    // that some filters will be reused.  In that case, they will not progress
    // again and will not contribute to the overall progress.
    // This means that the total progress may not reach 1.
    // For example, when run in 3D, m_DerivativeA will not run again when
    // dima = 0 and dimb = 2 because it did not change from when
    // dima = 0 and dimb = 1.
    const double weight =
      1.0 / ( ImageDimension * ( ImageDimension * ( ImageDimension + 1 ) / 2 ) );
  
    // note: this is not constant to suppress a warning
    unsigned int numberOfSmoothingFilters = NumberOfSmoothingFilters;
  
    for ( unsigned int i = 0; i < numberOfSmoothingFilters; i++ )
      {
      progress->RegisterInternalFilter(m_SmoothingFilters[i], weight);
      }
    progress->RegisterInternalFilter(m_DerivativeFilterA, weight);
    progress->RegisterInternalFilter(m_DerivativeFilterB, weight);
  
    const typename TInputImage::ConstPointer inputImage( this->GetInput() );
  
    m_ImageAdaptor->SetImage( this->GetOutput() );
  
    m_ImageAdaptor->SetLargestPossibleRegion(
      inputImage->GetLargestPossibleRegion() );
  
    m_ImageAdaptor->SetBufferedRegion(
      inputImage->GetBufferedRegion() );
  
    m_ImageAdaptor->SetRequestedRegion(
      inputImage->GetRequestedRegion() );
  
    m_ImageAdaptor->Allocate();
  
    m_DerivativeFilterA->SetInput(inputImage);
    m_DerivativeFilterB->SetInput( m_DerivativeFilterA->GetOutput() );
  
    unsigned int element = 0;
  
    for ( unsigned int dima = 0; dima < ImageDimension; dima++ )
      {
      for ( unsigned int dimb = dima; dimb < ImageDimension; dimb++ )
        {
        // Manage the diagonal in a different way in order to avoid
        // applying a double smoothing to this direction, and missing
        // to smooth one of the other directions.
        if ( dimb == dima )
          {
          m_DerivativeFilterA->SetOrder(DerivativeFilterAType::SecondOrder);
          m_DerivativeFilterB->SetOrder(DerivativeFilterBType::ZeroOrder);
  
          // only need the output of m_DerivativeFilterA once
          m_DerivativeFilterB->InPlaceOn();
  
          unsigned int i = 0;
          unsigned int j = 0;
          // find the direction for the first filter.
          while ( j < ImageDimension )
            {
            if ( j != dima )
              {
              m_DerivativeFilterB->SetDirection(j);
              j++;
              break;
              }
            j++;
            }
          // find the direction for all the other filters
          while ( i < numberOfSmoothingFilters )
            {
            while ( j < ImageDimension )
              {
              if ( j != dima )
                {
                m_SmoothingFilters[i]->SetDirection(j);
                j++;
                break;
                }
              j++;
              }
            i++;
            }
  
          m_DerivativeFilterA->SetDirection(dima);
          }
        else
          {
          m_DerivativeFilterA->SetOrder(DerivativeFilterAType::FirstOrder);
          m_DerivativeFilterB->SetOrder(DerivativeFilterBType::FirstOrder);
  
          if ( dimb < ImageDimension - 1 )
            {
            // can reuse the output of m_DerivativeFilterA
            m_DerivativeFilterB->InPlaceOff();
            }
          else
            {
            m_DerivativeFilterB->InPlaceOn();
            }
  
          unsigned int i = 0;
          unsigned int j = 0;
          while ( i < numberOfSmoothingFilters )
            {
            while ( j < ImageDimension )
              {
              if ( j != dima && j != dimb )
                {
                m_SmoothingFilters[i]->SetDirection(j);
                j++;
                break;
                }
              j++;
              }
            i++;
            }
  
          m_DerivativeFilterA->SetDirection(dima);
          m_DerivativeFilterB->SetDirection(dimb);
          }
  
        typename RealImageType::Pointer derivativeImage;
  
        // Deal with the 2D case.
        if ( numberOfSmoothingFilters > 0 )
          {
          int temp_dim = static_cast< int >( ImageDimension ) - 3;
          GaussianFilterPointer lastFilter = m_SmoothingFilters[temp_dim];
          lastFilter->UpdateLargestPossibleRegion();
          derivativeImage = lastFilter->GetOutput();
          }
        else
          {
          m_DerivativeFilterB->UpdateLargestPossibleRegion();
          derivativeImage = m_DerivativeFilterB->GetOutput();
          }
  
        // Copy the results to the corresponding component
        // on the output image of vectors
        m_ImageAdaptor->SelectNthElement(element++);
  
        ImageRegionIteratorWithIndex< RealImageType > it(
          derivativeImage,
          derivativeImage->GetRequestedRegion() );
  
        ImageRegionIteratorWithIndex< OutputImageAdaptorType > ot(
          m_ImageAdaptor,
          m_ImageAdaptor->GetRequestedRegion() );
  
        const RealType spacingA = inputImage->GetSpacing()[dima];
        const RealType spacingB = inputImage->GetSpacing()[dimb];
  
        const RealType factor = spacingA * spacingB;
  
        it.GoToBegin();
        ot.GoToBegin();
        while ( !it.IsAtEnd() )
          {
          ot.Set(it.Get() / factor);
          ++it;
          ++ot;
          }
  
        derivativeImage->ReleaseData();
        }
      }
  
    // manually release memory in last filter in the pipeline
    if (  numberOfSmoothingFilters > 0 )
      {
      m_SmoothingFilters[numberOfSmoothingFilters - 1]->GetOutput()->ReleaseData();
      }
    else
      {
      m_DerivativeFilterB->GetOutput()->ReleaseData();
      }
    m_DerivativeFilterA->GetOutput()->ReleaseData();
  }
  
  template< typename TInputImage, typename TOutputImage >
  void
  HessianRecursiveGaussianImageFilter< TInputImage, TOutputImage >
  ::PrintSelf(std::ostream & os, Indent indent) const
  {
    Superclass::PrintSelf(os, indent);
    os << "NormalizeAcrossScale: " << m_NormalizeAcrossScale << std::endl;
  }
  
  } // end namespace itk
  
  #endif
  