itkPhasedArray3DSpecialCoordinatesImage.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  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
 *
 *         https://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_MOVE(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);

  itkOverrideGetNameOfClassMacro(PhasedArray3DSpecialCoordinatesImage);

  using PixelType = TPixel;

  using ValueType = TPixel;

  using InternalPixelType = TPixel;
 
  using typename Superclass::IOPixelType;

  using AccessorType = DefaultPixelAccessor<PixelType>;

  using AccessorFunctorType = DefaultPixelAccessorFunctor<Self>;

  using NeighborhoodAccessorFunctorType = NeighborhoodAccessorFunctor<Self>;

  static constexpr unsigned int ImageDimension = 3;

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

  using typename Superclass::OffsetType;

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

  using PixelContainer = ImportImageContainer<SizeValueType, PixelType>;

  using typename Superclass::RegionType;

  using typename Superclass::SpacingType;

  using typename Superclass::PointType;

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

  template <typename TIndexRep, typename TCoordinate>
  [[nodiscard]] ContinuousIndex<TIndexRep, 3>
  TransformPhysicalPointToContinuousIndex(const Point<TCoordinate, 3> & point) 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
    TCoordinate azimuth = Math::pi_over_2;
    TCoordinate 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 TCoordinate radius = std::sqrt(point[0] * point[0] + point[1] * point[1] + point[2] * point[2]);
 
    // Convert the "proper" angular coordinates into index format
    ContinuousIndex<TIndexRep, 3> index;
    index[0] = static_cast<TCoordinate>((azimuth / m_AzimuthAngularSeparation) + (maxAzimuth / 2.0));
    index[1] = static_cast<TCoordinate>((elevation / m_ElevationAngularSeparation) + (maxElevation / 2.0));
    index[2] = static_cast<TCoordinate>(((radius - m_FirstSampleDistance) / m_RadiusSampleSize));
    return index;
  }

  template <typename TCoordinate, typename TIndexRep>
  ITK_NODISCARD("Call the overload which has the point as the only parameter and returns the index")
  bool TransformPhysicalPointToContinuousIndex(const Point<TCoordinate, 3> &   point,
                                               ContinuousIndex<TIndexRep, 3> & index) const
  {
    index = this->TransformPhysicalPointToContinuousIndex<TIndexRep>(point);
 
    // Now, check to see if the index is within allowed bounds
    const bool isInside = this->GetLargestPossibleRegion().IsInside(index);
 
    return isInside;
  }

  template <typename TCoordinate>
  [[nodiscard]] IndexType
  TransformPhysicalPointToIndex(const Point<TCoordinate, 3> & point) 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
    TCoordinate azimuth = Math::pi_over_2;
    TCoordinate 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 TCoordinate radius = std::sqrt(point[0] * point[0] + point[1] * point[1] + point[2] * point[2]);
 
    // Convert the "proper" angular coordinates into index format
    IndexType index;
    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));
    return index;
  }

  template <typename TCoordinate>
  ITK_NODISCARD("Call the overload which has the point as the only parameter and returns the index")
  bool TransformPhysicalPointToIndex(const Point<TCoordinate, 3> & point, IndexType & index) const
  {
    index = this->TransformPhysicalPointToIndex(point);
 
    // Now, check to see if the index is within allowed bounds
    const bool isInside = this->GetLargestPossibleRegion().IsInside(index);
 
    return isInside;
  }

  template <typename TCoordinate, typename TIndexRep>
  void
  TransformContinuousIndexToPhysicalPoint(const ContinuousIndex<TIndexRep, 3> & index,
                                          Point<TCoordinate, 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 TCoordinate azimuth = (index[0] - (maxAzimuth / 2.0)) * m_AzimuthAngularSeparation;
    const TCoordinate elevation = (index[1] - (maxElevation / 2.0)) * m_ElevationAngularSeparation;
    const TCoordinate radius = (index[2] * m_RadiusSampleSize) + m_FirstSampleDistance;
 
    // Convert the angular coordinates into Cartesian coordinates
    const TCoordinate tanOfAzimuth = std::tan(azimuth);
    const TCoordinate tanOfElevation = std::tan(elevation);
 
    point[2] =
      static_cast<TCoordinate>(radius / std::sqrt(1 + tanOfAzimuth * tanOfAzimuth + tanOfElevation * tanOfElevation));
    point[1] = static_cast<TCoordinate>(point[2] * tanOfElevation);
    point[0] = static_cast<TCoordinate>(point[2] * tanOfAzimuth);
  }

  template <typename TCoordinate, typename TIndexRep>
  [[nodiscard]] Point<TCoordinate, 3>
  TransformContinuousIndexToPhysicalPoint(const ContinuousIndex<TIndexRep, 3> & index) const
  {
    Point<TCoordinate, 3> point;
    this->TransformContinuousIndexToPhysicalPoint(index, point);
    return point;
  }


  template <typename TCoordinate>
  void
  TransformIndexToPhysicalPoint(const IndexType & index, Point<TCoordinate, 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 TCoordinate azimuth = (static_cast<double>(index[0]) - (maxAzimuth / 2.0)) * m_AzimuthAngularSeparation;
    const TCoordinate elevation = (static_cast<double>(index[1]) - (maxElevation / 2.0)) * m_ElevationAngularSeparation;
    const TCoordinate radius = (static_cast<double>(index[2]) * m_RadiusSampleSize) + m_FirstSampleDistance;
 
    // Convert the angular coordinates into Cartesian coordinates
    const TCoordinate tanOfAzimuth = std::tan(azimuth);
    const TCoordinate tanOfElevation = std::tan(elevation);
 
    point[2] =
      static_cast<TCoordinate>(radius / std::sqrt(1.0 + tanOfAzimuth * tanOfAzimuth + tanOfElevation * tanOfElevation));
    point[1] = static_cast<TCoordinate>(point[2] * tanOfElevation);
    point[0] = static_cast<TCoordinate>(point[2] * tanOfAzimuth);
  }

  template <typename TCoordinate>
  [[nodiscard]] Point<TCoordinate, 3>
  TransformIndexToPhysicalPoint(const IndexType & index) const
  {
    Point<TCoordinate, 3> point;
    this->TransformIndexToPhysicalPoint(index, point);
    return point;
  }


  itkSetMacro(AzimuthAngularSeparation, double);

  itkSetMacro(ElevationAngularSeparation, double);

  itkSetMacro(RadiusSampleSize, double);

  itkSetMacro(FirstSampleDistance, double);
 
  template <typename TCoordinate>
  void
  TransformLocalVectorToPhysicalVector(FixedArray<TCoordinate, 3> &) const
  {}
 
  template <typename TCoordinate>
  void
  TransformPhysicalVectorToLocalVector(const FixedArray<TCoordinate, 3> &, FixedArray<TCoordinate, 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