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
00085 template <class TPixel>
00086 class ITK_EXPORT PhasedArray3DSpecialCoordinatesImage :
00087 public SpecialCoordinatesImage<TPixel,3>
00088 {
00089 public:
00091 typedef PhasedArray3DSpecialCoordinatesImage Self;
00092 typedef SpecialCoordinatesImage<TPixel,3> Superclass;
00093 typedef SmartPointer<Self> Pointer;
00094 typedef SmartPointer<const Self> ConstPointer;
00095 typedef WeakPointer<const Self> ConstWeakPointer;
00096
00098 itkNewMacro(Self);
00099
00101 itkTypeMacro(PhasedArray3DSpecialCoordinatesImage, SpecialCoordinatesImage);
00102
00105 typedef TPixel PixelType;
00106
00108 typedef TPixel ValueType ;
00109
00114 typedef TPixel InternalPixelType;
00115
00116 typedef typename Superclass::IOPixelType IOPixelType;
00117
00120 typedef DefaultPixelAccessor< PixelType > AccessorType;
00121
00125 typedef DefaultPixelAccessorFunctor< Self > AccessorFunctorType;
00126
00131 itkStaticConstMacro(ImageDimension, unsigned int, 3);
00132
00134 typedef ImportImageContainer<unsigned long, PixelType> PixelContainer;
00135
00137 typedef typename Superclass::IndexType IndexType;
00138
00140 typedef typename Superclass::OffsetType OffsetType;
00141
00143 typedef typename Superclass::SizeType SizeType;
00144
00146 typedef typename Superclass::RegionType RegionType;
00147
00152 typedef typename Superclass::SpacingType SpacingType;
00153
00156 typedef typename Superclass::PointType PointType;
00157
00159 typedef typename PixelContainer::Pointer PixelContainerPointer;
00160 typedef typename PixelContainer::ConstPointer PixelContainerConstPointer;
00161
00166 template<class TCoordRep>
00167 bool TransformPhysicalPointToContinuousIndex(
00168 const Point<TCoordRep, 3>& point,
00169 ContinuousIndex<TCoordRep, 3>& index ) const
00170 {
00171 RegionType region = this->GetLargestPossibleRegion();
00172 double maxAzimuth = region.GetSize(0) - 1;
00173 double maxElevation = region.GetSize(1) - 1;
00174
00175
00176 TCoordRep azimuth = vcl_atan(point[0] / point[2]);
00177 TCoordRep elevation = vcl_atan(point[1] / point[2]);
00178 TCoordRep radius = vcl_sqrt(point[0] * point[0]
00179 + point[1] * point[1]
00180 + point[2] * point[2] );
00181
00182
00183 index[0] = static_cast<TCoordRep>( (azimuth/m_AzimuthAngularSeparation)
00184 + (maxAzimuth/2.0) );
00185 index[1] = static_cast<TCoordRep>( (elevation/m_ElevationAngularSeparation)
00186 + (maxElevation/2.0) );
00187 index[2] = static_cast<TCoordRep>( ( (radius-m_FirstSampleDistance)
00188 / m_RadiusSampleSize) );
00189
00190
00191 const bool isInside = region.IsInside( index );
00192
00193 return isInside;
00194 }
00195
00200 template<class TCoordRep>
00201 bool TransformPhysicalPointToIndex(
00202 const Point<TCoordRep, 3>& point,
00203 IndexType & index ) const
00204 {
00205 typedef typename IndexType::IndexValueType IndexValueType;
00206
00207 RegionType region = this->GetLargestPossibleRegion();
00208 double maxAzimuth = region.GetSize(0) - 1;
00209 double maxElevation = region.GetSize(1) - 1;
00210
00211
00212 TCoordRep azimuth = vcl_atan(point[0] / point[2]);
00213 TCoordRep elevation = vcl_atan(point[1] / point[2]);
00214 TCoordRep radius = vcl_sqrt(point[0] * point[0]
00215 + point[1] * point[1]
00216 + point[2] * point[2] );
00217
00218
00219 index[0] = static_cast<IndexValueType>( (azimuth/m_AzimuthAngularSeparation)
00220 + (maxAzimuth/2.0) );
00221 index[1] = static_cast<IndexValueType>( (elevation/m_ElevationAngularSeparation)
00222 + (maxElevation/2.0) );
00223 index[2] = static_cast<IndexValueType>( ( (radius-m_FirstSampleDistance)
00224 / m_RadiusSampleSize ) );
00225
00226
00227 const bool isInside = region.IsInside( index );
00228
00229 return isInside;
00230 }
00231
00236 template<class TCoordRep>
00237 void TransformContinuousIndexToPhysicalPoint(
00238 const ContinuousIndex<TCoordRep, 3>& index,
00239 Point<TCoordRep, 3>& point ) const
00240 {
00241 RegionType region = this->GetLargestPossibleRegion();
00242 double maxAzimuth = region.GetSize(0) - 1;
00243 double maxElevation = region.GetSize(1) - 1;
00245
00246
00247 TCoordRep azimuth = ( index[0] - (maxAzimuth/2.0) )
00248 * m_AzimuthAngularSeparation;
00249 TCoordRep elevation = ( index[1] - (maxElevation/2.0) )
00250 * m_ElevationAngularSeparation;
00251 TCoordRep radius = (index[2]*m_RadiusSampleSize)+m_FirstSampleDistance;
00252
00253
00254 TCoordRep tanOfAzimuth = vcl_tan(azimuth);
00255 TCoordRep tanOfElevation = vcl_tan(elevation);
00256 point[2] = static_cast<TCoordRep>( radius /
00257 vcl_sqrt(1 + tanOfAzimuth*tanOfAzimuth + tanOfElevation*tanOfElevation));
00258 point[1] = static_cast<TCoordRep>( point[2] * tanOfElevation );
00259 point[0] = static_cast<TCoordRep>( point[2] * tanOfAzimuth );
00260 }
00261
00267 template<class TCoordRep>
00268 void TransformIndexToPhysicalPoint(
00269 const IndexType & index,
00270 Point<TCoordRep, 3>& point ) const
00271 {
00272 RegionType region = this->GetLargestPossibleRegion();
00273 double maxAzimuth = region.GetSize(0) - 1;
00274 double maxElevation = region.GetSize(1) - 1;
00276
00277
00278 TCoordRep azimuth = ( static_cast<double>(index[0]) - (maxAzimuth/2.0) )
00279 * m_AzimuthAngularSeparation;
00280 TCoordRep elevation = ( static_cast<double>(index[1]) - (maxElevation/2.0) )
00281 * m_ElevationAngularSeparation;
00282 TCoordRep radius = (static_cast<double>(index[2]) * m_RadiusSampleSize)
00283 + m_FirstSampleDistance;
00284
00285
00286 TCoordRep tanOfAzimuth = vcl_tan(azimuth);
00287 TCoordRep tanOfElevation = vcl_tan(elevation);
00288 point[2] = static_cast<TCoordRep>( radius / vcl_sqrt(
00289 1.0 + tanOfAzimuth*tanOfAzimuth + tanOfElevation*tanOfElevation) );
00290 point[1] = static_cast<TCoordRep>( point[2] * tanOfElevation );
00291 point[0] = static_cast<TCoordRep>( point[2] * tanOfAzimuth );
00292 }
00293
00294
00296 itkSetMacro(AzimuthAngularSeparation, double);
00297
00299 itkSetMacro(ElevationAngularSeparation, double);
00300
00302 itkSetMacro(RadiusSampleSize, double);
00303
00305 itkSetMacro(FirstSampleDistance, double);
00306
00307 protected:
00308 PhasedArray3DSpecialCoordinatesImage()
00309 {
00310 m_RadiusSampleSize = 1;
00311 m_AzimuthAngularSeparation = 1 * (2.0*vnl_math::pi/360.0);
00312 m_ElevationAngularSeparation = 1 * (2.0*vnl_math::pi/360.0);
00313 m_FirstSampleDistance = 0;
00314 }
00315 virtual ~PhasedArray3DSpecialCoordinatesImage() {};
00316 void PrintSelf(std::ostream& os, Indent indent) const;
00317
00318 private:
00319 PhasedArray3DSpecialCoordinatesImage(const Self&);
00320 void operator=(const Self&);
00321
00322 double m_AzimuthAngularSeparation;
00323 double m_ElevationAngularSeparation;
00324 double m_RadiusSampleSize;
00325 double m_FirstSampleDistance;
00326
00327 };
00328 }
00329 #ifndef ITK_MANUAL_INSTANTIATION
00330 #include "itkPhasedArray3DSpecialCoordinatesImage.txx"
00331 #endif
00332
00333 #endif
00334
00335