ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkPhasedArray3DSpecialCoordinatesImage.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 #ifndef __itkPhasedArray3DSpecialCoordinatesImage_h
00019 #define __itkPhasedArray3DSpecialCoordinatesImage_h
00020 
00021 #include "itkSpecialCoordinatesImage.h"
00022 #include "itkPoint.h"
00023 #include "vnl/vnl_math.h"
00024 
00025 namespace itk
00026 {
00089 template< class TPixel >
00090 class ITK_EXPORT PhasedArray3DSpecialCoordinatesImage:
00091   public SpecialCoordinatesImage< TPixel, 3 >
00092 {
00093 public:
00095   typedef PhasedArray3DSpecialCoordinatesImage Self;
00096   typedef SpecialCoordinatesImage< TPixel, 3 > Superclass;
00097   typedef SmartPointer< Self >                 Pointer;
00098   typedef SmartPointer< const Self >           ConstPointer;
00099   typedef WeakPointer< const Self >            ConstWeakPointer;
00100 
00102   itkNewMacro(Self);
00103 
00105   itkTypeMacro(PhasedArray3DSpecialCoordinatesImage, SpecialCoordinatesImage);
00106 
00109   typedef TPixel PixelType;
00110 
00112   typedef TPixel ValueType;
00113 
00118   typedef TPixel InternalPixelType;
00119 
00120   typedef typename Superclass::IOPixelType IOPixelType;
00121 
00124   typedef DefaultPixelAccessor< PixelType > AccessorType;
00125 
00129   typedef DefaultPixelAccessorFunctor< Self > AccessorFunctorType;
00130 
00135   itkStaticConstMacro(ImageDimension, unsigned int, 3);
00136 
00138   typedef typename Superclass::IndexType      IndexType;
00139   typedef typename Superclass::IndexValueType IndexValueType;
00140 
00142   typedef typename Superclass::OffsetType OffsetType;
00143 
00145   typedef typename Superclass::SizeType      SizeType;
00146   typedef typename Superclass::SizeValueType SizeValueType;
00147 
00149   typedef ImportImageContainer< SizeValueType, PixelType > PixelContainer;
00150 
00154   typedef typename Superclass::RegionType RegionType;
00155 
00162   typedef typename Superclass::SpacingType SpacingType;
00163 
00168   typedef typename Superclass::PointType PointType;
00169 
00171   typedef typename PixelContainer::Pointer      PixelContainerPointer;
00172   typedef typename PixelContainer::ConstPointer PixelContainerConstPointer;
00173 
00178   template< class TCoordRep >
00179   bool TransformPhysicalPointToContinuousIndex(
00180     const Point< TCoordRep, 3 > & point,
00181     ContinuousIndex< TCoordRep, 3 > & index) const
00182   {
00183     RegionType region = this->GetLargestPossibleRegion();
00184     double     maxAzimuth =    region.GetSize(0) - 1;
00185     double     maxElevation =  region.GetSize(1) - 1;
00186 
00187     // Convert Cartesian coordinates into angular coordinates
00188     TCoordRep azimuth   = vcl_atan(point[0] / point[2]);
00189     TCoordRep elevation = vcl_atan(point[1] / point[2]);
00190     TCoordRep radius    = vcl_sqrt(point[0] * point[0]
00191                                    + point[1] * point[1]
00192                                    + point[2] * point[2]);
00193 
00194     // Convert the "proper" angular coordinates into index format
00195     index[0] = static_cast< TCoordRep >( ( azimuth / m_AzimuthAngularSeparation )
00196                                          + ( maxAzimuth / 2.0 ) );
00197     index[1] = static_cast< TCoordRep >( ( elevation / m_ElevationAngularSeparation )
00198                                          + ( maxElevation / 2.0 ) );
00199     index[2] = static_cast< TCoordRep >( ( ( radius - m_FirstSampleDistance )
00200                                            / m_RadiusSampleSize ) );
00201 
00202     // Now, check to see if the index is within allowed bounds
00203     const bool isInside = region.IsInside(index);
00204 
00205     return isInside;
00206   }
00207 
00212   template< class TCoordRep >
00213   bool TransformPhysicalPointToIndex(
00214     const Point< TCoordRep, 3 > & point,
00215     IndexType & index) const
00216   {
00217     RegionType region = this->GetLargestPossibleRegion();
00218     double     maxAzimuth =    region.GetSize(0) - 1;
00219     double     maxElevation =  region.GetSize(1) - 1;
00221 
00222     // Convert Cartesian coordinates into angular coordinates
00223     TCoordRep azimuth   = vcl_atan(point[0] / point[2]);
00224     TCoordRep elevation = vcl_atan(point[1] / point[2]);
00225     TCoordRep radius    = vcl_sqrt(point[0] * point[0]
00226                                    + point[1] * point[1]
00227                                    + point[2] * point[2]);
00228 
00229     // Convert the "proper" angular coordinates into index format
00230     index[0] = static_cast< IndexValueType >(
00231       ( azimuth / m_AzimuthAngularSeparation )
00232       + ( maxAzimuth / 2.0 ) );
00233     index[1] = static_cast< IndexValueType >(
00234       ( elevation / m_ElevationAngularSeparation )
00235       + ( maxElevation / 2.0 ) );
00236     index[2] = static_cast< IndexValueType >(
00237       ( ( radius - m_FirstSampleDistance )
00238         / m_RadiusSampleSize ) );
00239 
00240     // Now, check to see if the index is within allowed bounds
00241     const bool isInside = region.IsInside(index);
00242 
00243     return isInside;
00244   }
00245 
00250   template< class TCoordRep >
00251   void TransformContinuousIndexToPhysicalPoint(
00252     const ContinuousIndex< TCoordRep, 3 > & index,
00253     Point< TCoordRep, 3 > & point) const
00254   {
00255     RegionType region = this->GetLargestPossibleRegion();
00256     double     maxAzimuth =    region.GetSize(0) - 1;
00257     double     maxElevation =  region.GetSize(1) - 1;
00259 
00260     // Convert the index into proper angular coordinates
00261     TCoordRep azimuth   = ( index[0] - ( maxAzimuth / 2.0 ) )
00262                           * m_AzimuthAngularSeparation;
00263     TCoordRep elevation = ( index[1] - ( maxElevation / 2.0 ) )
00264                           * m_ElevationAngularSeparation;
00265     TCoordRep radius    = ( index[2] * m_RadiusSampleSize ) + m_FirstSampleDistance;
00266 
00267     // Convert the angular coordinates into Cartesian coordinates
00268     TCoordRep tanOfAzimuth    = vcl_tan(azimuth);
00269     TCoordRep tanOfElevation  = vcl_tan(elevation);
00270 
00271     point[2] = static_cast< TCoordRep >( radius
00272                                          / vcl_sqrt(1
00273                                                     + tanOfAzimuth * tanOfAzimuth
00274                                                     + tanOfElevation * tanOfElevation) );
00275     point[1] = static_cast< TCoordRep >( point[2] * tanOfElevation );
00276     point[0] = static_cast< TCoordRep >( point[2] * tanOfAzimuth );
00277   }
00278 
00284   template< class TCoordRep >
00285   void TransformIndexToPhysicalPoint(
00286     const IndexType & index,
00287     Point< TCoordRep, 3 > & point) const
00288   {
00289     RegionType region = this->GetLargestPossibleRegion();
00290     double     maxAzimuth =    region.GetSize(0) - 1;
00291     double     maxElevation =  region.GetSize(1) - 1;
00293 
00294     // Convert the index into proper angular coordinates
00295     TCoordRep azimuth =
00296       ( static_cast< double >( index[0] ) - ( maxAzimuth / 2.0 ) )
00297       * m_AzimuthAngularSeparation;
00298     TCoordRep elevation =
00299       ( static_cast< double >( index[1] ) - ( maxElevation / 2.0 ) )
00300       * m_ElevationAngularSeparation;
00301     TCoordRep radius =
00302       ( static_cast< double >( index[2] ) * m_RadiusSampleSize )
00303       + m_FirstSampleDistance;
00304 
00305     // Convert the angular coordinates into Cartesian coordinates
00306     TCoordRep tanOfAzimuth    = vcl_tan(azimuth);
00307     TCoordRep tanOfElevation  = vcl_tan(elevation);
00308 
00309     point[2] = static_cast< TCoordRep >(
00310       radius / vcl_sqrt(
00311         1.0 + tanOfAzimuth * tanOfAzimuth + tanOfElevation * tanOfElevation) );
00312     point[1] = static_cast< TCoordRep >( point[2] * tanOfElevation );
00313     point[0] = static_cast< TCoordRep >( point[2] * tanOfAzimuth );
00314   }
00315 
00317   itkSetMacro(AzimuthAngularSeparation, double);
00318 
00320   itkSetMacro(ElevationAngularSeparation, double);
00321 
00323   itkSetMacro(RadiusSampleSize, double);
00324 
00326   itkSetMacro(FirstSampleDistance, double);
00327 
00328   template< class TCoordRep >
00329   void TransformLocalVectorToPhysicalVector(
00330     FixedArray< TCoordRep, 3 > & ) const
00331     {}
00332 
00333   template< class TCoordRep >
00334   void TransformPhysicalVectorToLocalVector(
00335     const FixedArray< TCoordRep, 3 > & ,
00336     FixedArray< TCoordRep, 3 > & ) const
00337     {}
00338 
00339 protected:
00340   PhasedArray3DSpecialCoordinatesImage()
00341   {
00342     m_RadiusSampleSize = 1;
00343     m_AzimuthAngularSeparation   =  1 * ( 2.0 * vnl_math::pi / 360.0 ); // 1
00344                                                                         // degree
00345     m_ElevationAngularSeparation =  1 * ( 2.0 * vnl_math::pi / 360.0 ); // 1
00346                                                                         // degree
00347     m_FirstSampleDistance = 0;
00348   }
00349 
00350   virtual ~PhasedArray3DSpecialCoordinatesImage() {}
00351   void PrintSelf(std::ostream & os, Indent indent) const;
00352 
00353 private:
00354   PhasedArray3DSpecialCoordinatesImage(const Self &); //purposely not
00355                                                       // implemented
00356   void operator=(const Self &);                       //purposely not
00357                                                       // implemented
00358 
00359   double m_AzimuthAngularSeparation;    // in radians
00360   double m_ElevationAngularSeparation;  // in radians
00361   double m_RadiusSampleSize;
00362   double m_FirstSampleDistance;
00363 };
00364 } // end namespace itk
00365 
00366 // Define instantiation macro for this template.
00367 #define ITK_TEMPLATE_PhasedArray3DSpecialCoordinatesImage(_, EXPORT, TypeX, TypeY)     \
00368   namespace itk                                                                        \
00369   {                                                                                    \
00370   _( 1 ( class EXPORT PhasedArray3DSpecialCoordinatesImage< ITK_TEMPLATE_1 TypeX > ) ) \
00371   namespace Templates                                                                  \
00372   {                                                                                    \
00373   typedef PhasedArray3DSpecialCoordinatesImage< ITK_TEMPLATE_1 TypeX >                 \
00374   PhasedArray3DSpecialCoordinatesImage##TypeY;                                       \
00375   }                                                                                    \
00376   }
00377 
00378 #if ITK_TEMPLATE_EXPLICIT
00379 #include "Templates/itkPhasedArray3DSpecialCoordinatesImage+-.h"
00380 #endif
00381 
00382 #if ITK_TEMPLATE_TXX
00383 #include "itkPhasedArray3DSpecialCoordinatesImage.hxx"
00384 #endif
00385 
00386 #endif
00387