ITK  4.3.0
Insight Segmentation and Registration Toolkit
itkPhasedArray3DSpecialCoordinatesImage.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef __itkPhasedArray3DSpecialCoordinatesImage_h
19 #define __itkPhasedArray3DSpecialCoordinatesImage_h
20 
22 #include "itkPoint.h"
23 #include "vnl/vnl_math.h"
24 
25 namespace itk
26 {
92 template< class TPixel >
94  public SpecialCoordinatesImage< TPixel, 3 >
95 {
96 public:
103 
105  itkNewMacro(Self);
106 
109 
112  typedef TPixel PixelType;
113 
115  typedef TPixel ValueType;
116 
121  typedef TPixel InternalPixelType;
122 
123  typedef typename Superclass::IOPixelType IOPixelType;
124 
128 
133 
138  itkStaticConstMacro(ImageDimension, unsigned int, 3);
139 
141  typedef typename Superclass::IndexType IndexType;
143 
145  typedef typename Superclass::OffsetType OffsetType;
146 
148  typedef typename Superclass::SizeType SizeType;
150 
153 
157  typedef typename Superclass::RegionType RegionType;
158 
165  typedef typename Superclass::SpacingType SpacingType;
166 
171  typedef typename Superclass::PointType PointType;
172 
176 
181  template< class TCoordRep >
182  bool TransformPhysicalPointToContinuousIndex(
183  const Point< TCoordRep, 3 > & point,
184  ContinuousIndex< TCoordRep, 3 > & index) const
185  {
186  RegionType region = this->GetLargestPossibleRegion();
187  double maxAzimuth = region.GetSize(0) - 1;
188  double maxElevation = region.GetSize(1) - 1;
189 
190  // Convert Cartesian coordinates into angular coordinates
191  TCoordRep azimuth = vcl_atan(point[0] / point[2]);
192  TCoordRep elevation = vcl_atan(point[1] / point[2]);
193  TCoordRep radius = vcl_sqrt(point[0] * point[0]
194  + point[1] * point[1]
195  + point[2] * point[2]);
196 
197  // Convert the "proper" angular coordinates into index format
198  index[0] = static_cast< TCoordRep >( ( azimuth / m_AzimuthAngularSeparation )
199  + ( maxAzimuth / 2.0 ) );
200  index[1] = static_cast< TCoordRep >( ( elevation / m_ElevationAngularSeparation )
201  + ( maxElevation / 2.0 ) );
202  index[2] = static_cast< TCoordRep >( ( ( radius - m_FirstSampleDistance )
203  / m_RadiusSampleSize ) );
204 
205  // Now, check to see if the index is within allowed bounds
206  const bool isInside = region.IsInside(index);
207 
208  return isInside;
209  }
210 
215  template< class TCoordRep >
216  bool TransformPhysicalPointToIndex(
217  const Point< TCoordRep, 3 > & point,
218  IndexType & index) const
219  {
220  RegionType region = this->GetLargestPossibleRegion();
221  double maxAzimuth = region.GetSize(0) - 1;
222  double maxElevation = region.GetSize(1) - 1;
224 
225  // Convert Cartesian coordinates into angular coordinates
226  TCoordRep azimuth = vcl_atan(point[0] / point[2]);
227  TCoordRep elevation = vcl_atan(point[1] / point[2]);
228  TCoordRep radius = vcl_sqrt(point[0] * point[0]
229  + point[1] * point[1]
230  + point[2] * point[2]);
231 
232  // Convert the "proper" angular coordinates into index format
233  index[0] = static_cast< IndexValueType >(
234  ( azimuth / m_AzimuthAngularSeparation )
235  + ( maxAzimuth / 2.0 ) );
236  index[1] = static_cast< IndexValueType >(
237  ( elevation / m_ElevationAngularSeparation )
238  + ( maxElevation / 2.0 ) );
239  index[2] = static_cast< IndexValueType >(
240  ( ( radius - m_FirstSampleDistance )
241  / m_RadiusSampleSize ) );
242 
243  // Now, check to see if the index is within allowed bounds
244  const bool isInside = region.IsInside(index);
245 
246  return isInside;
247  }
248 
253  template< class TCoordRep >
254  void TransformContinuousIndexToPhysicalPoint(
255  const ContinuousIndex< TCoordRep, 3 > & index,
256  Point< TCoordRep, 3 > & point) const
257  {
258  RegionType region = this->GetLargestPossibleRegion();
259  double maxAzimuth = region.GetSize(0) - 1;
260  double maxElevation = region.GetSize(1) - 1;
262 
263  // Convert the index into proper angular coordinates
264  TCoordRep azimuth = ( index[0] - ( maxAzimuth / 2.0 ) )
265  * m_AzimuthAngularSeparation;
266  TCoordRep elevation = ( index[1] - ( maxElevation / 2.0 ) )
267  * m_ElevationAngularSeparation;
268  TCoordRep radius = ( index[2] * m_RadiusSampleSize ) + m_FirstSampleDistance;
269 
270  // Convert the angular coordinates into Cartesian coordinates
271  TCoordRep tanOfAzimuth = vcl_tan(azimuth);
272  TCoordRep tanOfElevation = vcl_tan(elevation);
273 
274  point[2] = static_cast< TCoordRep >( radius
275  / vcl_sqrt(1
276  + tanOfAzimuth * tanOfAzimuth
277  + tanOfElevation * tanOfElevation) );
278  point[1] = static_cast< TCoordRep >( point[2] * tanOfElevation );
279  point[0] = static_cast< TCoordRep >( point[2] * tanOfAzimuth );
280  }
281 
287  template< class TCoordRep >
288  void TransformIndexToPhysicalPoint(
289  const IndexType & index,
290  Point< TCoordRep, 3 > & point) const
291  {
292  RegionType region = this->GetLargestPossibleRegion();
293  double maxAzimuth = region.GetSize(0) - 1;
294  double maxElevation = region.GetSize(1) - 1;
296 
297  // Convert the index into proper angular coordinates
298  TCoordRep azimuth =
299  ( static_cast< double >( index[0] ) - ( maxAzimuth / 2.0 ) )
300  * m_AzimuthAngularSeparation;
301  TCoordRep elevation =
302  ( static_cast< double >( index[1] ) - ( maxElevation / 2.0 ) )
303  * m_ElevationAngularSeparation;
304  TCoordRep radius =
305  ( static_cast< double >( index[2] ) * m_RadiusSampleSize )
306  + m_FirstSampleDistance;
307 
308  // Convert the angular coordinates into Cartesian coordinates
309  TCoordRep tanOfAzimuth = vcl_tan(azimuth);
310  TCoordRep tanOfElevation = vcl_tan(elevation);
311 
312  point[2] = static_cast< TCoordRep >(
313  radius / vcl_sqrt(
314  1.0 + tanOfAzimuth * tanOfAzimuth + tanOfElevation * tanOfElevation) );
315  point[1] = static_cast< TCoordRep >( point[2] * tanOfElevation );
316  point[0] = static_cast< TCoordRep >( point[2] * tanOfAzimuth );
317  }
318 
320  itkSetMacro(AzimuthAngularSeparation, double);
321 
323  itkSetMacro(ElevationAngularSeparation, double);
324 
326  itkSetMacro(RadiusSampleSize, double);
327 
329  itkSetMacro(FirstSampleDistance, double);
330 
331  template< class TCoordRep >
332  void TransformLocalVectorToPhysicalVector(
334  {}
335 
336  template< class TCoordRep >
337  void TransformPhysicalVectorToLocalVector(
340  {}
341 
342 protected:
344  {
345  m_RadiusSampleSize = 1;
346  m_AzimuthAngularSeparation = 1 * ( 2.0 * vnl_math::pi / 360.0 ); // 1
347  // degree
348  m_ElevationAngularSeparation = 1 * ( 2.0 * vnl_math::pi / 360.0 ); // 1
349  // degree
350  m_FirstSampleDistance = 0;
351  }
352 
354  void PrintSelf(std::ostream & os, Indent indent) const;
355 
356 private:
357  PhasedArray3DSpecialCoordinatesImage(const Self &); //purposely not
358  // implemented
359  void operator=(const Self &); //purposely not
360  // implemented
361 
362  double m_AzimuthAngularSeparation; // in radians
363  double m_ElevationAngularSeparation; // in radians
366 };
367 } // end namespace itk
368 
369 #ifndef ITK_MANUAL_INSTANTIATION
370 #include "itkPhasedArray3DSpecialCoordinatesImage.hxx"
371 #endif
372 
373 #endif
374