00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifndef __itkPhasedArray3DSpecialCoordinatesImage_h
00018 #define __itkPhasedArray3DSpecialCoordinatesImage_h
00019
00020 #include "itkSpecialCoordinatesImage.h"
00021 #include "itkImageRegion.h"
00022 #include "itkPoint.h"
00023 #include "itkContinuousIndex.h"
00024 #include "vnl/vnl_math.h"
00025
00026
00027 namespace itk
00028 {
00029
00090 template <class TPixel>
00091 class ITK_EXPORT PhasedArray3DSpecialCoordinatesImage :
00092 public SpecialCoordinatesImage<TPixel,3>
00093 {
00094 public:
00096 typedef PhasedArray3DSpecialCoordinatesImage Self;
00097 typedef SpecialCoordinatesImage<TPixel,3> Superclass;
00098 typedef SmartPointer<Self> Pointer;
00099 typedef SmartPointer<const Self> ConstPointer;
00100 typedef WeakPointer<const Self> ConstWeakPointer;
00101
00103 itkNewMacro(Self);
00104
00106 itkTypeMacro(PhasedArray3DSpecialCoordinatesImage, SpecialCoordinatesImage);
00107
00110 typedef TPixel PixelType;
00111
00113 typedef TPixel ValueType;
00114
00119 typedef TPixel InternalPixelType;
00120
00121 typedef typename Superclass::IOPixelType IOPixelType;
00122
00125 typedef DefaultPixelAccessor< PixelType > AccessorType;
00126
00130 typedef DefaultPixelAccessorFunctor< Self > AccessorFunctorType;
00131
00136 itkStaticConstMacro(ImageDimension, unsigned int, 3);
00137
00139 typedef ImportImageContainer<unsigned long, PixelType> PixelContainer;
00140
00142 typedef typename Superclass::IndexType IndexType;
00143
00145 typedef typename Superclass::OffsetType OffsetType;
00146
00148 typedef typename Superclass::SizeType SizeType;
00149
00153 typedef typename Superclass::RegionType RegionType;
00154
00161 typedef typename Superclass::SpacingType SpacingType;
00162
00167 typedef typename Superclass::PointType PointType;
00168
00170 typedef typename PixelContainer::Pointer PixelContainerPointer;
00171 typedef typename PixelContainer::ConstPointer PixelContainerConstPointer;
00172
00177 template<class TCoordRep>
00178 bool TransformPhysicalPointToContinuousIndex(
00179 const Point<TCoordRep, 3>& point,
00180 ContinuousIndex<TCoordRep, 3>& index ) const
00181 {
00182 RegionType region = this->GetLargestPossibleRegion();
00183 double maxAzimuth = region.GetSize(0) - 1;
00184 double maxElevation = region.GetSize(1) - 1;
00185
00186
00187 TCoordRep azimuth = vcl_atan(point[0] / point[2]);
00188 TCoordRep elevation = vcl_atan(point[1] / point[2]);
00189 TCoordRep radius = vcl_sqrt(point[0] * point[0]
00190 + point[1] * point[1]
00191 + point[2] * point[2] );
00192
00193
00194 index[0] = static_cast<TCoordRep>( (azimuth/m_AzimuthAngularSeparation)
00195 + (maxAzimuth/2.0) );
00196 index[1] = static_cast<TCoordRep>( (elevation/m_ElevationAngularSeparation)
00197 + (maxElevation/2.0) );
00198 index[2] = static_cast<TCoordRep>( ( (radius-m_FirstSampleDistance)
00199 / m_RadiusSampleSize) );
00200
00201
00202 const bool isInside = region.IsInside( index );
00203
00204 return isInside;
00205 }
00206
00211 template<class TCoordRep>
00212 bool TransformPhysicalPointToIndex(
00213 const Point<TCoordRep, 3>& point,
00214 IndexType & index ) const
00215 {
00216 typedef typename IndexType::IndexValueType IndexValueType;
00217
00218 RegionType region = this->GetLargestPossibleRegion();
00219 double maxAzimuth = region.GetSize(0) - 1;
00220 double maxElevation = region.GetSize(1) - 1;
00221
00222
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
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
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
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
00268 TCoordRep tanOfAzimuth = vcl_tan(azimuth);
00269 TCoordRep tanOfElevation = vcl_tan(elevation);
00270 point[2] = static_cast<TCoordRep>( radius /
00271 vcl_sqrt(1 +
00272 tanOfAzimuth*tanOfAzimuth +
00273 tanOfElevation*tanOfElevation));
00274 point[1] = static_cast<TCoordRep>( point[2] * tanOfElevation );
00275 point[0] = static_cast<TCoordRep>( point[2] * tanOfAzimuth );
00276 }
00277
00283 template<class TCoordRep>
00284 void TransformIndexToPhysicalPoint(
00285 const IndexType & index,
00286 Point<TCoordRep, 3>& point ) const
00287 {
00288 RegionType region = this->GetLargestPossibleRegion();
00289 double maxAzimuth = region.GetSize(0) - 1;
00290 double maxElevation = region.GetSize(1) - 1;
00292
00293
00294 TCoordRep azimuth =
00295 (static_cast<double>(index[0]) - (maxAzimuth/2.0) )
00296 * m_AzimuthAngularSeparation;
00297 TCoordRep elevation =
00298 (static_cast<double>(index[1]) - (maxElevation/2.0) )
00299 * m_ElevationAngularSeparation;
00300 TCoordRep radius =
00301 (static_cast<double>(index[2]) * m_RadiusSampleSize)
00302 + m_FirstSampleDistance;
00303
00304
00305 TCoordRep tanOfAzimuth = vcl_tan(azimuth);
00306 TCoordRep tanOfElevation = vcl_tan(elevation);
00307 point[2] = static_cast<TCoordRep>(
00308 radius / vcl_sqrt(
00309 1.0 + tanOfAzimuth*tanOfAzimuth + tanOfElevation*tanOfElevation) );
00310 point[1] = static_cast<TCoordRep>( point[2] * tanOfElevation );
00311 point[0] = static_cast<TCoordRep>( point[2] * tanOfAzimuth );
00312 }
00313
00314
00316 itkSetMacro(AzimuthAngularSeparation, double);
00317
00319 itkSetMacro(ElevationAngularSeparation, double);
00320
00322 itkSetMacro(RadiusSampleSize, double);
00323
00325 itkSetMacro(FirstSampleDistance, double);
00326
00327 protected:
00328 PhasedArray3DSpecialCoordinatesImage()
00329 {
00330 m_RadiusSampleSize = 1;
00331 m_AzimuthAngularSeparation = 1 * (2.0*vnl_math::pi/360.0);
00332 m_ElevationAngularSeparation = 1 * (2.0*vnl_math::pi/360.0);
00333 m_FirstSampleDistance = 0;
00334 }
00335 virtual ~PhasedArray3DSpecialCoordinatesImage() {};
00336 void PrintSelf(std::ostream& os, Indent indent) const;
00337
00338 private:
00339 PhasedArray3DSpecialCoordinatesImage(const Self&);
00340 void operator=(const Self&);
00341
00342 double m_AzimuthAngularSeparation;
00343 double m_ElevationAngularSeparation;
00344 double m_RadiusSampleSize;
00345 double m_FirstSampleDistance;
00346
00347 };
00348 }
00349
00350
00351 #define ITK_TEMPLATE_PhasedArray3DSpecialCoordinatesImage(_, EXPORT, x, y) namespace itk { \
00352 _(1(class EXPORT PhasedArray3DSpecialCoordinatesImage< ITK_TEMPLATE_1 x >)) \
00353 namespace Templates { typedef PhasedArray3DSpecialCoordinatesImage< ITK_TEMPLATE_1 x > \
00354 PhasedArray3DSpecialCoordinatesImage##y; } \
00355 }
00356
00357 #if ITK_TEMPLATE_EXPLICIT
00358 # include "Templates/itkPhasedArray3DSpecialCoordinatesImage+-.h"
00359 #endif
00360
00361 #if ITK_TEMPLATE_TXX
00362 # include "itkPhasedArray3DSpecialCoordinatesImage.txx"
00363 #endif
00364
00365 #endif
00366