itkPhasedArray3DSpecialCoordinatesImage.h

Go to the documentation of this file.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkPhasedArray3DSpecialCoordinatesImage.h,v $
  Language:  C++
  Date:      $Date: 2009-02-19 19:41:23 $
  Version:   $Revision: 1.22 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __itkPhasedArray3DSpecialCoordinatesImage_h
#define __itkPhasedArray3DSpecialCoordinatesImage_h

#include "itkSpecialCoordinatesImage.h"
#include "itkImageRegion.h"
#include "itkPoint.h"
#include "itkContinuousIndex.h"
#include "vnl/vnl_math.h"


namespace itk
{

template <class TPixel>
class ITK_EXPORT PhasedArray3DSpecialCoordinatesImage :
public SpecialCoordinatesImage<TPixel,3>
{
public:
  typedef PhasedArray3DSpecialCoordinatesImage Self;
  typedef SpecialCoordinatesImage<TPixel,3>    Superclass;
  typedef SmartPointer<Self>                   Pointer;
  typedef SmartPointer<const Self>             ConstPointer;
  typedef WeakPointer<const Self>              ConstWeakPointer;

  itkNewMacro(Self);

  itkTypeMacro(PhasedArray3DSpecialCoordinatesImage, SpecialCoordinatesImage);

  typedef TPixel PixelType;

  typedef TPixel ValueType;

  typedef TPixel InternalPixelType;

  typedef typename Superclass::IOPixelType   IOPixelType;
  
  typedef DefaultPixelAccessor< PixelType > AccessorType;

  typedef DefaultPixelAccessorFunctor< Self > AccessorFunctorType;

  itkStaticConstMacro(ImageDimension, unsigned int, 3);

  typedef ImportImageContainer<unsigned long, PixelType> PixelContainer;

  typedef typename Superclass::IndexType      IndexType;
  typedef typename Superclass::IndexValueType IndexValueType;

  typedef typename Superclass::OffsetType OffsetType;

  typedef typename Superclass::SizeType   SizeType;

  typedef typename Superclass::RegionType RegionType;

  typedef typename Superclass::SpacingType SpacingType;

  typedef typename Superclass::PointType PointType;

  typedef typename PixelContainer::Pointer      PixelContainerPointer;
  typedef typename PixelContainer::ConstPointer PixelContainerConstPointer;

  template<class TCoordRep>
  bool TransformPhysicalPointToContinuousIndex(
              const Point<TCoordRep, 3>& point,
              ContinuousIndex<TCoordRep, 3>& index   ) const
    {
    RegionType region = this->GetLargestPossibleRegion();
    double maxAzimuth =    region.GetSize(0) - 1;
    double maxElevation =  region.GetSize(1) - 1;
    
    // Convert Cartesian coordinates into angular coordinates
    TCoordRep azimuth   = vcl_atan(point[0] / point[2]);
    TCoordRep elevation = vcl_atan(point[1] / point[2]);
    TCoordRep radius    = vcl_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<class TCoordRep>
  bool TransformPhysicalPointToIndex(
            const Point<TCoordRep, 3>& point,
            IndexType & index                                ) const
    {
    RegionType region = this->GetLargestPossibleRegion();
    double maxAzimuth =    region.GetSize(0) - 1;
    double maxElevation =  region.GetSize(1) - 1;

    // Convert Cartesian coordinates into angular coordinates
    TCoordRep azimuth   = vcl_atan(point[0] / point[2]);
    TCoordRep elevation = vcl_atan(point[1] / point[2]);
    TCoordRep radius    = vcl_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<class TCoordRep>
  void TransformContinuousIndexToPhysicalPoint(
            const ContinuousIndex<TCoordRep, 3>& index,
            Point<TCoordRep, 3>& point        ) const
    {
    RegionType region = this->GetLargestPossibleRegion();
    double maxAzimuth =    region.GetSize(0) - 1;
    double maxElevation =  region.GetSize(1) - 1;

    // Convert the index into proper angular coordinates
    TCoordRep azimuth   = ( index[0] - (maxAzimuth/2.0) )
                          * m_AzimuthAngularSeparation;
    TCoordRep elevation = ( index[1] - (maxElevation/2.0) )
                          * m_ElevationAngularSeparation;
    TCoordRep radius    = (index[2]*m_RadiusSampleSize)+m_FirstSampleDistance;
    
    // Convert the angular coordinates into Cartesian coordinates
    TCoordRep tanOfAzimuth    = vcl_tan(azimuth);
    TCoordRep tanOfElevation  = vcl_tan(elevation);
    point[2] = static_cast<TCoordRep>( radius /
           vcl_sqrt(1 +
                    tanOfAzimuth*tanOfAzimuth +
                    tanOfElevation*tanOfElevation));
    point[1] = static_cast<TCoordRep>( point[2] * tanOfElevation );
    point[0] = static_cast<TCoordRep>( point[2] * tanOfAzimuth );
    }

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

    // Convert the index into proper angular coordinates
    TCoordRep azimuth = 
      (static_cast<double>(index[0]) - (maxAzimuth/2.0) )
      * m_AzimuthAngularSeparation;
    TCoordRep elevation =
      (static_cast<double>(index[1]) - (maxElevation/2.0) )
      * m_ElevationAngularSeparation;
    TCoordRep radius =
      (static_cast<double>(index[2]) * m_RadiusSampleSize)
      + m_FirstSampleDistance;
    
    // Convert the angular coordinates into Cartesian coordinates
    TCoordRep tanOfAzimuth    = vcl_tan(azimuth);
    TCoordRep tanOfElevation  = vcl_tan(elevation);
    point[2] = static_cast<TCoordRep>(
      radius / vcl_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);

#ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION
  template<class TCoordRep>
  void TransformLocalVectorToPhysicalVector(
    const FixedArray<TCoordRep, 3> & inputGradient,
          FixedArray<TCoordRep, 3> & outputGradient ) const
    {
    }
#endif
protected:
  PhasedArray3DSpecialCoordinatesImage()
    {
    m_RadiusSampleSize = 1;
    m_AzimuthAngularSeparation   =  1 * (2.0*vnl_math::pi/360.0); // 1 degree
    m_ElevationAngularSeparation =  1 * (2.0*vnl_math::pi/360.0); // 1 degree
    m_FirstSampleDistance = 0;
    }
  virtual ~PhasedArray3DSpecialCoordinatesImage() {};
  void PrintSelf(std::ostream& os, Indent indent) const;
  
private:
  PhasedArray3DSpecialCoordinatesImage(const Self&);//purposely not implemented
  void operator=(const Self&); //purposely not implemented
  
  double  m_AzimuthAngularSeparation;   // in radians
  double  m_ElevationAngularSeparation; // in radians
  double  m_RadiusSampleSize;
  double  m_FirstSampleDistance;
  
};
} // end namespace itk

// Define instantiation macro for this template.
#define ITK_TEMPLATE_PhasedArray3DSpecialCoordinatesImage(_, EXPORT, x, y) namespace itk { \
  _(1(class EXPORT PhasedArray3DSpecialCoordinatesImage< ITK_TEMPLATE_1 x >)) \
  namespace Templates { typedef PhasedArray3DSpecialCoordinatesImage< ITK_TEMPLATE_1 x > \
                                         PhasedArray3DSpecialCoordinatesImage##y; } \
  }

#if ITK_TEMPLATE_EXPLICIT
# include "Templates/itkPhasedArray3DSpecialCoordinatesImage+-.h"
#endif

#if ITK_TEMPLATE_TXX
# include "itkPhasedArray3DSpecialCoordinatesImage.txx"
#endif

#endif