itkFrequencyHalfHermitianFFTLayoutImageRegionConstIteratorWithIndex.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  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 itkFrequencyHalfHermitianFFTLayoutImageRegionConstIteratorWithIndex_h
#define itkFrequencyHalfHermitianFFTLayoutImageRegionConstIteratorWithIndex_h

#include "itkImageRegionConstIteratorWithIndex.h"

namespace itk
{
template< typename TImage >
class FrequencyHalfHermitianFFTLayoutImageRegionConstIteratorWithIndex:
  public ImageRegionConstIteratorWithIndex< TImage >
{
public:
  using Self = FrequencyHalfHermitianFFTLayoutImageRegionConstIteratorWithIndex;
  using Superclass = ImageRegionConstIteratorWithIndex< TImage >;

  using IndexType = typename Superclass::IndexType;
  using SizeType = typename Superclass::SizeType;
  using OffsetType = typename Superclass::OffsetType;
  using RegionType = typename Superclass::RegionType;
  using ImageType = typename Superclass::ImageType;
  using PixelContainer = typename Superclass::PixelContainer;
  using PixelContainerPointer = typename Superclass::PixelContainerPointer;
  using InternalPixelType = typename Superclass::InternalPixelType;
  using PixelType = typename Superclass::PixelType;
  using AccessorType = typename Superclass::AccessorType;

  using FrequencyType = typename ImageType::SpacingType;
  using FrequencyValueType = typename ImageType::SpacingValueType;
  FrequencyHalfHermitianFFTLayoutImageRegionConstIteratorWithIndex() :
    ImageRegionConstIteratorWithIndex< TImage >()

  {
    this->Init();
  }

  FrequencyHalfHermitianFFTLayoutImageRegionConstIteratorWithIndex(
    const TImage *ptr, const RegionType & region) :
    ImageRegionConstIteratorWithIndex< TImage >(ptr, region),
    m_ActualXDimensionIsOdd(false)
  {
    this->Init();
  }

  explicit FrequencyHalfHermitianFFTLayoutImageRegionConstIteratorWithIndex(
    const Superclass & it) :
    ImageRegionConstIteratorWithIndex< TImage >(it),
    m_ActualXDimensionIsOdd(false)
  {
    this->Init();
  }

  /*
   * Image Index [0, N - 1] returns [0 to N/2] (positive) union [-N/2 + 1, -1] (negative).
   * So index N/2 + 1 returns the bin -N/2 + 1.
   * If first index of the image is not zero, it stills returns values in the same range.
   * f = [0, 1, ...,   N/2-1,     -N/2, ..., -1]  if N is even
   * f = [0, 1, ..., (N-1)/2, -(N-1)/2, ..., -1]  if N is odd
   */
  IndexType GetFrequencyBin() const
  {
    IndexType freqInd;
    freqInd.Fill(0);
    for (unsigned int dim = 0; dim < TImage::ImageDimension; dim++)
      {
      if (this->m_PositionIndex[dim] <= m_LargestPositiveFrequencyIndex[dim])
        {
        freqInd[dim] = this->m_PositionIndex[dim] - this->m_MinIndex[dim];
        }
      else //  -. From -N/2 + 1 (Nyquist if even) to -1 (-df in frequency)
        {
        freqInd[dim] = this->m_PositionIndex[dim] - (this->m_MaxIndex[dim] + 1);
        }
      }
    return freqInd;
  }

  FrequencyType GetFrequency() const
  {
    FrequencyType freq;
    IndexType     freqInd = this->GetFrequencyBin();

    for (unsigned int dim = 0; dim < TImage::ImageDimension; dim++)
      {
      freq[dim] = this->m_FrequencyOrigin[dim]
        + this->m_FrequencySpacing[dim] * freqInd[dim];
      }
    return freq;
  }

  FrequencyValueType GetFrequencyModuloSquare() const
  {
    FrequencyValueType w2(0);
    FrequencyType      w( this->GetFrequency() );

    for (unsigned int dim = 0; dim < TImage::ImageDimension; dim++)
      {
      w2 += w[dim] * w[dim];
      }
    return w2;
  }

  itkGetConstReferenceMacro(LargestPositiveFrequencyIndex, IndexType);

  itkGetConstReferenceMacro(MinIndex, IndexType);

  itkGetConstReferenceMacro(MaxIndex, IndexType);

  itkGetConstReferenceMacro(FrequencyOrigin, FrequencyType);

  itkGetConstReferenceMacro(FrequencySpacing, FrequencyType);

  void SetActualXDimensionIsOdd(bool value)
    {
    this->m_ActualXDimensionIsOdd = value;
    SizeType sizeImage =
      this->m_Image->GetLargestPossibleRegion().GetSize();
    auto size_estimated = 2 * ( sizeImage[0] - 1);
    size_estimated += this->GetActualXDimensionIsOdd() ? 1 : 0;
    this->m_FrequencySpacing[0] = 1.0 / (this->m_Image->GetSpacing()[0] * size_estimated);
    }
  itkGetMacro(ActualXDimensionIsOdd, bool);
  itkBooleanMacro(ActualXDimensionIsOdd);

private:
  void Init()
  {
    SizeType sizeImage =
      this->m_Image->GetLargestPossibleRegion().GetSize();
    this->m_MinIndex =
      this->m_Image->GetLargestPossibleRegion().GetIndex();
    this->m_MaxIndex =
      this->m_Image->GetLargestPossibleRegion().GetUpperIndex();
    for (unsigned int dim = 0; dim < ImageType::ImageDimension; dim++)
      {
      this->m_LargestPositiveFrequencyIndex[dim] = static_cast<FrequencyValueType>(
          this->m_MinIndex[dim] + sizeImage[dim] / 2 );
      // Set frequency metadata.
      // Origin of frequencies is zero after a FFT.
      this->m_FrequencyOrigin[dim] = 0.0;
      // SamplingFrequency = 1.0 / SpatialImageSpacing
      // Freq_BinSize = SamplingFrequency / Size
      this->m_FrequencySpacing[dim] = 1.0 / (this->m_Image->GetSpacing()[dim] * sizeImage[dim]);
      }
    // Corrections for Hermitian
    // The fastest index has no negative frequencies.
    this->m_LargestPositiveFrequencyIndex[0] = static_cast<FrequencyValueType>(
          this->m_MaxIndex[0]);
    // In fastest dimension only:
    // Size_estimated_original = 2 * (Size_current - 1 )
    // Where size_current is the size of the output of RealToHalfHermitianFFT
    // Ex: Size Original = 10, Current = 6, Estimated = 10.
    //     Size Original = 11, Current = 6, Estimated = 10.
    auto size_estimated = 2 * ( sizeImage[0] - 1);
    size_estimated += this->GetActualXDimensionIsOdd() ? 1 : 0;
    this->m_FrequencySpacing[0] = 1.0 / (this->m_Image->GetSpacing()[0] * size_estimated);
  }

  IndexType     m_LargestPositiveFrequencyIndex;
  IndexType     m_MinIndex;
  IndexType     m_MaxIndex;
  FrequencyType m_FrequencyOrigin;
  FrequencyType m_FrequencySpacing;
  bool          m_ActualXDimensionIsOdd{false};
};
} // end namespace itk
#endif