itkPhasedArray3DSpecialCoordinatesImage.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 itkPhasedArray3DSpecialCoordinatesImage_h
#define itkPhasedArray3DSpecialCoordinatesImage_h

#include "itkSpecialCoordinatesImage.h"
#include "itkPoint.h"
#include "itkMath.h"
#include "itkDefaultPixelAccessor.h"
#include "itkNeighborhoodAccessorFunctor.h"

namespace itk
{
template< typename TPixel >
class ITK_TEMPLATE_EXPORT PhasedArray3DSpecialCoordinatesImage:
  public SpecialCoordinatesImage< TPixel, 3 >
{
public:
  ITK_DISALLOW_COPY_AND_ASSIGN(PhasedArray3DSpecialCoordinatesImage);

  using Self = PhasedArray3DSpecialCoordinatesImage;
  using Superclass = SpecialCoordinatesImage< TPixel, 3 >;
  using Pointer = SmartPointer< Self >;
  using ConstPointer = SmartPointer< const Self >;
  using ConstWeakPointer = WeakPointer< const Self >;

  itkNewMacro(Self);

  itkTypeMacro(PhasedArray3DSpecialCoordinatesImage, SpecialCoordinatesImage);

  using PixelType = TPixel;

  using ValueType = TPixel;

  using InternalPixelType = TPixel;

  using IOPixelType = typename Superclass::IOPixelType;

  using AccessorType = DefaultPixelAccessor< PixelType >;

  using AccessorFunctorType = DefaultPixelAccessorFunctor< Self >;

  using NeighborhoodAccessorFunctorType = NeighborhoodAccessorFunctor< Self >;

  static constexpr unsigned int ImageDimension = 3;

  using IndexType = typename Superclass::IndexType;
  using IndexValueType = typename Superclass::IndexValueType;

  using OffsetType = typename Superclass::OffsetType;

  using SizeType = typename Superclass::SizeType;
  using SizeValueType = typename Superclass::SizeValueType;

  using PixelContainer = ImportImageContainer< SizeValueType, PixelType >;

  using RegionType = typename Superclass::RegionType;

  using SpacingType = typename Superclass::SpacingType;

  using PointType = typename Superclass::PointType;

  using PixelContainerPointer = typename PixelContainer::Pointer;
  using PixelContainerConstPointer = typename PixelContainer::ConstPointer;

  template< typename TCoordRep, typename TIndexRep >
  bool TransformPhysicalPointToContinuousIndex(
    const Point< TCoordRep, 3 > & point,
    ContinuousIndex< TIndexRep, 3 > & index) const
  {
    const RegionType region = this->GetLargestPossibleRegion();
    const double maxAzimuth = region.GetSize(0) - 1;
    const double maxElevation = region.GetSize(1) - 1;

    // Convert Cartesian coordinates into angular coordinates
    TCoordRep azimuth   = Math::pi_over_2;
    TCoordRep elevation = Math::pi_over_2;
    if( point[2] != 0.0 )
      {
      azimuth   = std::atan(point[0] / point[2]);
      elevation = std::atan(point[1] / point[2]);
      }
    const TCoordRep radius  = std::sqrt(point[0] * point[0]
                                   + point[1] * point[1]
                                   + point[2] * point[2]);

    // Convert the "proper" angular coordinates into index format
    index[0] = static_cast< TCoordRep >( ( azimuth / m_AzimuthAngularSeparation )
                                         + ( maxAzimuth / 2.0 ) );
    index[1] = static_cast< TCoordRep >( ( elevation / m_ElevationAngularSeparation )
                                         + ( maxElevation / 2.0 ) );
    index[2] = static_cast< TCoordRep >( ( ( radius - m_FirstSampleDistance )
                                           / m_RadiusSampleSize ) );

    // Now, check to see if the index is within allowed bounds
    const bool isInside = region.IsInside(index);

    return isInside;
  }

  template< typename TCoordRep >
  bool TransformPhysicalPointToIndex(
    const Point< TCoordRep, 3 > & point,
    IndexType & index) const
  {
    const RegionType region = this->GetLargestPossibleRegion();
    const double maxAzimuth   = region.GetSize(0) - 1;
    const double maxElevation = region.GetSize(1) - 1;

    // Convert Cartesian coordinates into angular coordinates
    TCoordRep azimuth   = Math::pi_over_2;
    TCoordRep elevation = Math::pi_over_2;
    if( point[2] != 0.0 )
      {
      azimuth   = std::atan(point[0] / point[2]);
      elevation = std::atan(point[1] / point[2]);
      }
    const TCoordRep radius = std::sqrt(point[0] * point[0]
                                   + point[1] * point[1]
                                   + point[2] * point[2]);

    // Convert the "proper" angular coordinates into index format
    index[0] = static_cast< IndexValueType >(
      ( azimuth / m_AzimuthAngularSeparation )
      + ( maxAzimuth / 2.0 ) );
    index[1] = static_cast< IndexValueType >(
      ( elevation / m_ElevationAngularSeparation )
      + ( maxElevation / 2.0 ) );
    index[2] = static_cast< IndexValueType >(
      ( ( radius - m_FirstSampleDistance )
        / m_RadiusSampleSize ) );

    // Now, check to see if the index is within allowed bounds
    const bool isInside = region.IsInside(index);

    return isInside;
  }

  template< typename TCoordRep, typename TIndexRep >
  void TransformContinuousIndexToPhysicalPoint(
    const ContinuousIndex< TIndexRep, 3 > & index,
    Point< TCoordRep, 3 > & point) const
  {
    const RegionType region = this->GetLargestPossibleRegion();
    const double maxAzimuth   = region.GetSize(0) - 1;
    const double maxElevation = region.GetSize(1) - 1;

    // Convert the index into proper angular coordinates
    const TCoordRep azimuth   = ( index[0] - ( maxAzimuth / 2.0 ) )
                          * m_AzimuthAngularSeparation;
    const TCoordRep elevation = ( index[1] - ( maxElevation / 2.0 ) )
                          * m_ElevationAngularSeparation;
    const TCoordRep radius    = ( index[2] * m_RadiusSampleSize ) + m_FirstSampleDistance;

    // Convert the angular coordinates into Cartesian coordinates
    const TCoordRep tanOfAzimuth   = std::tan(azimuth);
    const TCoordRep tanOfElevation = std::tan(elevation);

    point[2] = static_cast< TCoordRep >( radius
                                         / std::sqrt(1
                                                    + tanOfAzimuth * tanOfAzimuth
                                                    + tanOfElevation * tanOfElevation) );
    point[1] = static_cast< TCoordRep >( point[2] * tanOfElevation );
    point[0] = static_cast< TCoordRep >( point[2] * tanOfAzimuth );
  }

  template< typename TCoordRep >
  void TransformIndexToPhysicalPoint(
    const IndexType & index,
    Point< TCoordRep, 3 > & point) const
  {
    const RegionType region = this->GetLargestPossibleRegion();
    const double maxAzimuth   = region.GetSize(0) - 1;
    const double maxElevation = region.GetSize(1) - 1;

    // Convert the index into proper angular coordinates
    const TCoordRep azimuth =
      ( static_cast< double >( index[0] ) - ( maxAzimuth / 2.0 ) )
      * m_AzimuthAngularSeparation;
    const TCoordRep elevation =
      ( static_cast< double >( index[1] ) - ( maxElevation / 2.0 ) )
      * m_ElevationAngularSeparation;
    const TCoordRep radius =
      ( static_cast< double >( index[2] ) * m_RadiusSampleSize )
      + m_FirstSampleDistance;

    // Convert the angular coordinates into Cartesian coordinates
    const TCoordRep tanOfAzimuth   = std::tan(azimuth);
    const TCoordRep tanOfElevation = std::tan(elevation);

    point[2] = static_cast< TCoordRep >(
      radius / std::sqrt(
        1.0 + tanOfAzimuth * tanOfAzimuth + tanOfElevation * tanOfElevation) );
    point[1] = static_cast< TCoordRep >( point[2] * tanOfElevation );
    point[0] = static_cast< TCoordRep >( point[2] * tanOfAzimuth );
  }

  itkSetMacro(AzimuthAngularSeparation, double);

  itkSetMacro(ElevationAngularSeparation, double);

  itkSetMacro(RadiusSampleSize, double);

  itkSetMacro(FirstSampleDistance, double);

  template< typename TCoordRep >
  void TransformLocalVectorToPhysicalVector(
    FixedArray< TCoordRep, 3 > & ) const
    {}

  template< typename TCoordRep >
  void TransformPhysicalVectorToLocalVector(
    const FixedArray< TCoordRep, 3 > & ,
    FixedArray< TCoordRep, 3 > & ) const
    {}

  AccessorType GetPixelAccessor()
  { return AccessorType(); }

  const AccessorType GetPixelAccessor() const
  { return AccessorType(); }

  NeighborhoodAccessorFunctorType GetNeighborhoodAccessor()
  { return NeighborhoodAccessorFunctorType(); }

  const NeighborhoodAccessorFunctorType GetNeighborhoodAccessor() const
  { return NeighborhoodAccessorFunctorType(); }

protected:
  PhasedArray3DSpecialCoordinatesImage()
  {
    m_RadiusSampleSize = 1;
    m_AzimuthAngularSeparation   =  1 * ( 2.0 * itk::Math::pi / 360.0 ); // 1
                                                                        // degree
    m_ElevationAngularSeparation =  1 * ( 2.0 * itk::Math::pi / 360.0 ); // 1
                                                                        // degree
    m_FirstSampleDistance = 0;
  }

  ~PhasedArray3DSpecialCoordinatesImage() override = default;
  void PrintSelf(std::ostream & os, Indent indent) const override;

private:
  double m_AzimuthAngularSeparation;    // in radians
  double m_ElevationAngularSeparation;  // in radians
  double m_RadiusSampleSize;
  double m_FirstSampleDistance;
};
} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkPhasedArray3DSpecialCoordinatesImage.hxx"
#endif

#endif