itkPhasedArray3DSpecialCoordinatesImage.h

Go to the documentation of this file.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkPhasedArray3DSpecialCoordinatesImage.h,v $
  Language:  C++
  Date:      $Date: 2005/09/25 19:02:54 $
  Version:   $Revision: 1.11 $

  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 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::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   = atan(point[0] / point[2]);
    TCoordRep elevation = atan(point[1] / point[2]);
    TCoordRep radius    = 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
    {
    typedef typename IndexType::IndexValueType IndexValueType;
    
    RegionType region = this->GetLargestPossibleRegion();
    double maxAzimuth =    region.GetSize(0) - 1;
    double maxElevation =  region.GetSize(1) - 1;
    
    // Convert Cartesian coordinates into angular coordinates
    TCoordRep azimuth   = atan(point[0] / point[2]);
    TCoordRep elevation = atan(point[1] / point[2]);
    TCoordRep radius    = 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    = tan(azimuth);
    TCoordRep tanOfElevation  = tan(elevation);
    point[2] = static_cast<TCoordRep>( radius /
           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    = tan(azimuth);
    TCoordRep tanOfElevation  = tan(elevation);
    point[2] = static_cast<TCoordRep>( radius / 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);
  
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;
  
};
#ifdef ITK_EXPLICIT_INSTANTIATION
   extern template class PhasedArray3DSpecialCoordinatesImage<float         >;
   extern template class PhasedArray3DSpecialCoordinatesImage<double        >;
   extern template class PhasedArray3DSpecialCoordinatesImage<unsigned char >;
   extern template class PhasedArray3DSpecialCoordinatesImage<unsigned short>;
   extern template class PhasedArray3DSpecialCoordinatesImage<unsigned int  >;
   extern template class PhasedArray3DSpecialCoordinatesImage<signed char   >;
   extern template class PhasedArray3DSpecialCoordinatesImage<signed short  >;
   extern template class PhasedArray3DSpecialCoordinatesImage<signed int    >;
#endif
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkPhasedArray3DSpecialCoordinatesImage.txx"
#endif

#endif

---