← Back to branch summary

~ubuntu-branches/ubuntu/precise/insighttoolkit/precise

« back to all changes in this revision

Viewing changes to Code/Algorithms/itkRayCastInterpolateImageFunction.txx

  • Committer: Bazaar Package Importer
  • Author(s): Steve M. Robbins
  • Date: 2008-05-31 12:07:29 UTC
  • mfrom: (3.1.3 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080531120729-7g7layn480le43ko
Tags: 3.6.0-3
http://bugs.debian.org/478500
http://bugs.debian.org/478940
http://bugs.debian.org/483312
* debian/patches/gccxml-workaround.patch: New.  Work around gccxml issue
  with #include_next; c.f. http://www.gccxml.org/Bug/view.php?id=7134.  
* debian/patches/gcc43.patch: include <cstring> in itkNeighbourhood.h.
  This only showed up in the tcl wrapping step.

* Above two entries fix FTBFS for GCC 4.3-based systems.
  Closes: #478500.

* debian/patches/sharedforward.patch: New.  Ensure that linux/sparc
  systems are not also configured as a SUN sparc system, which requires
  SUN header sys/isa_defs.h.  Closes: #478940, #483312.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Code/Algorithms/itkRayCastInterpolateImageFunction.txx
3
3
Program:   Insight Segmentation & Registration Toolkit
4
4
Module:    $RCSfile: itkRayCastInterpolateImageFunction.txx,v $
5
5
Language:  C++
6
 
Date:      $Date: 2006/03/19 04:36:55 $
7
 
Version:   $Revision: 1.19 $
 
6
Date:      $Date: 2007-04-26 20:20:40 $
 
7
Version:   $Revision: 1.20 $
8
8
 
9
9
Copyright (c) Insight Software Consortium. All rights reserved.
10
10
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11
11
 
12
 
This software is distributed WITHOUT ANY WARRANTY; without even 
13
 
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
 
12
This software is distributed WITHOUT ANY WARRANTY; without even
 
13
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14
14
PURPOSE.  See the above copyright notices for more information.
15
15
 
16
16
=========================================================================*/
17
17
 
18
 
#ifndef _itkRayCastInterpolateImageFunction_txx
19
 
#define _itkRayCastInterpolateImageFunction_txx
 
18
#ifndef __itkRayCastInterpolateImageFunction_txx
 
19
#define __itkRayCastInterpolateImageFunction_txx
20
20
 
21
21
#include "itkRayCastInterpolateImageFunction.h"
22
22
 
27
27
// exposed to the user
28
28
namespace {
29
29
 
30
 
/** Helper class to maintain state when casting a ray.  This helper
31
 
 * class keeps the RayCastInterpolateImageFunction thread safe.
 
30
/** \class Helper class to maintain state when casting a ray.
 
31
 *  This helper class keeps the RayCastInterpolateImageFunction thread safe.
32
32
 */
33
33
template <class TInputImage, class TCoordRep = float>
34
34
class RayCastHelper
38
38
  itkStaticConstMacro(InputImageDimension, unsigned int,
39
39
                      TInputImage::ImageDimension);
40
40
 
41
 
  /** 
42
 
   * Type of the Transform Base class 
 
41
  /**
 
42
   * Type of the Transform Base class
43
43
   * The fixed image should be a 3D image
44
44
   */
45
45
  typedef itk::Transform<TCoordRep,3,3> TransformType;
54
54
  typedef itk::Vector<TCoordRep, 3>                  DirectionType;
55
55
  typedef itk::Point<TCoordRep, 3>                   PointType;
56
56
 
57
 
  typedef TInputImage InputImageType;
 
57
  typedef TInputImage                                InputImageType;
58
58
  typedef typename InputImageType::PixelType         PixelType;
59
59
  typedef typename InputImageType::IndexType         IndexType;
60
60
 
63
63
   */
64
64
  void SetImage(const InputImageType *input)
65
65
    {
66
 
      m_Image = input;
 
66
    m_Image = input;
67
67
    }
68
 
    
69
 
  /** 
 
68
 
 
69
  /**
70
70
   *  Initialise the ray using the position and direction of a line.
71
71
   *
72
72
   *  \param RayPosn       The position of the ray in 3D (mm).
75
75
   *  \return True if this is a valid ray.
76
76
   */
77
77
  bool SetRay(OutputPointType RayPosn, DirectionType RayDirn);
78
 
    
 
78
 
79
79
 
80
80
  /** \brief
81
81
   *  Integrate the interpolated intensities along the ray and
88
88
   *
89
89
   *  \param integral      The integrated intensities along the ray.
90
90
   *
91
 
   * \return True if a valid ray was specified.  
 
91
   * \return True if a valid ray was specified.
92
92
   */
93
 
  bool Integrate(double &integral) {
 
93
  bool Integrate(double &integral)
 
94
    {
94
95
    return IntegrateAboveThreshold(integral, 0);
95
 
  };
96
 
    
 
96
    };
 
97
 
97
98
 
98
99
  /** \brief
99
100
   * Integrate the interpolated intensities above a given threshold,
100
 
   * along the ray and return the result. 
 
101
   * along the ray and return the result.
101
102
   *
102
103
   * This routine can be called after instantiating the ray and
103
104
   * calling SetProjectionCoord2D() or Reset(). It may then be called
105
106
   * coordinates.
106
107
   *
107
108
   * \param integral      The integrated intensities along the ray.
108
 
   * \param threshold     The integration threshold [default value: 0] 
 
109
   * \param threshold     The integration threshold [default value: 0]
109
110
   *
110
 
   * \return True if a valid ray was specified.  
 
111
   * \return True if a valid ray was specified.
111
112
   */
112
113
  bool IntegrateAboveThreshold(double &integral, double threshold);
113
114
 
121
122
 
122
123
  /// Reset the iterator to the start of the ray.
123
124
  void Reset(void);
124
 
    
 
125
 
125
126
  /// Return the interpolated intensity of the current ray point.
126
127
  double GetCurrentIntensity(void) const;
127
 
    
 
128
 
128
129
  /// Return the ray point spacing in mm
129
130
  double GetRayPointSpacing(void) const {
130
131
    typename InputImageType::SpacingType spacing=this->m_Image->GetSpacing();
131
 
      
 
132
 
132
133
    if (m_ValidRay)
133
134
      return vcl_sqrt(m_VoxelIncrement[0]*spacing[0]*m_VoxelIncrement[0]*spacing[0]
134
135
                    + m_VoxelIncrement[1]*spacing[1]*m_VoxelIncrement[1]*spacing[1]
137
138
      return 0.;
138
139
  };
139
140
 
140
 
  /// Set the initial zero state of the object 
 
141
  /// Set the initial zero state of the object
141
142
  void ZeroState();
142
143
 
143
144
  /// Initialise the object
144
145
  void Initialise(void);
145
 
    
 
146
 
146
147
protected:
147
148
  /// Calculate the endpoint coordinats of the ray in voxels.
148
149
  void EndPointsInVoxels(void);
149
 
    
150
 
  /** 
 
150
 
 
151
  /**
151
152
   * Calculate the incremental direction vector in voxels, 'dVoxel',
152
 
   * required to traverse the ray. 
 
153
   * required to traverse the ray.
153
154
   */
154
155
  void CalcDirnVector(void);
155
 
    
156
 
  /** 
 
156
 
 
157
  /**
157
158
   * Reduce the length of the ray until both start and end
158
159
   * coordinates lie inside the volume.
159
160
   *
160
161
   * \return True if a valid ray has been, false otherwise.
161
162
   */
162
163
  bool AdjustRayLength(void);
163
 
    
164
 
  /** 
 
164
 
 
165
  /**
165
166
   *   Obtain pointers to the four voxels surrounding the point where the ray
166
 
   *   enters the volume. 
167
 
   */ 
 
167
   *   enters the volume.
 
168
   */
168
169
  void InitialiseVoxelPointers(void);
169
 
    
 
170
 
170
171
  /// Increment the voxel pointers surrounding the current point on the ray.
171
172
  void IncrementVoxelPointers(void);
172
 
    
 
173
 
173
174
  /// Record volume dimensions and resolution
174
175
  void RecordVolumeDimensions(void);
175
 
    
 
176
 
176
177
  /// Define the corners of the volume
177
178
  void DefineCorners(void);
178
 
    
 
179
 
179
180
  /** \brief
180
181
   * Calculate the planes which define the volume.
181
182
   *
185
186
   * slopes of the lines which go to make up the volume( defined as
186
187
   * lines in cube x,y,z dirn and then each of these lines has a slope
187
188
   * in the world x,y,z dirn [3]) and finally also to return the length
188
 
   * of the sides of the lines in mm. 
 
189
   * of the sides of the lines in mm.
189
190
   */
190
191
  void CalcPlanesAndCorners(void);
191
 
    
 
192
 
192
193
  /** \brief
193
194
   *  Calculate the ray intercepts with the volume.
194
195
   *
195
196
   *  See where the ray cuts the volume, check that truncation does not occur,
196
 
   *  if not, then start ray where it first intercepts the volume and set 
 
197
   *  if not, then start ray where it first intercepts the volume and set
197
198
   *  x_max to be where it leaves the volume.
198
199
   *
199
200
   *  \return True if a valid ray has been specified, false otherwise.
206
207
   *   intercepted.
207
208
   */
208
209
  typedef enum {
209
 
    UNDEFINED_DIRECTION=0,        //!< Undefined                         
 
210
    UNDEFINED_DIRECTION=0,        //!< Undefined
210
211
    TRANSVERSE_IN_X,              //!< x
211
212
    TRANSVERSE_IN_Y,              //!< y
212
213
    TRANSVERSE_IN_Z,              //!< z
217
218
  // efficiency. This inner class will go in/out of scope with every
218
219
  // call to Evaluate()
219
220
  const InputImageType *m_Image;
220
 
    
 
221
 
221
222
  /// Flag indicating whether the current ray is valid
222
223
  bool m_ValidRay;
223
 
    
 
224
 
224
225
  /** \brief
225
 
   * The start position of the ray in voxels. 
 
226
   * The start position of the ray in voxels.
226
227
   *
227
228
   * NB. Two of the components of this coordinate (i.e. those lying within
228
229
   * the planes of voxels being traversed) will be shifted by half a
230
231
   * to be determined by simply casting to 'int' and optionally adding 1.
231
232
   */
232
233
  double m_RayVoxelStartPosition[3];
233
 
    
 
234
 
234
235
  /** \brief
235
236
   * The end coordinate of the ray in voxels.
236
237
   *
240
241
   * to be determined by simply casting to 'int' and optionally adding 1.
241
242
   */
242
243
  double m_RayVoxelEndPosition[3];
243
 
    
244
 
    
 
244
 
 
245
 
245
246
  /** \brief
246
247
   * The current coordinate on the ray in voxels.
247
248
   *
251
252
   * to be determined by simply casting to 'int' and optionally adding 1.
252
253
   */
253
254
  double m_Position3Dvox[3];
254
 
    
 
255
 
255
256
  /** The incremental direction vector of the ray in voxels. */
256
257
  double m_VoxelIncrement[3];
257
 
    
 
258
 
258
259
  /// The direction in which the ray is incremented thorough the volume (x, y or z).
259
260
  TraversalDirection m_TraversalDirection;
260
 
    
 
261
 
261
262
  /// The total number of planes of voxels traversed by the ray.
262
263
  int m_TotalRayVoxelPlanes;
263
 
    
 
264
 
264
265
  /// The current number of planes of voxels traversed by the ray.
265
266
  int m_NumVoxelPlanesTraversed;
266
 
    
 
267
 
267
268
  /// Pointers to the current four voxels surrounding the ray's trajectory.
268
269
  const PixelType *m_RayIntersectionVoxels[4];
269
 
    
 
270
 
270
271
  /**
271
272
   * The voxel coordinate of the bottom-left voxel of the current
272
 
   * four voxels surrounding the ray's trajectory. 
 
273
   * four voxels surrounding the ray's trajectory.
273
274
   */
274
275
  int m_RayIntersectionVoxelIndex[3];
275
 
    
 
276
 
276
277
  /// The dimension in voxels of the 3D volume in along the x axis
277
278
  int m_NumberOfVoxelsInX;
278
279
  /// The dimension in voxels of the 3D volume in along the y axis
279
280
  int m_NumberOfVoxelsInY;
280
281
  /// The dimension in voxels of the 3D volume in along the z axis
281
282
  int m_NumberOfVoxelsInZ;
282
 
    
 
283
 
283
284
  /// Voxel dimension in x
284
285
  double m_VoxelDimensionInX;
285
286
  /// Voxel dimension in y
286
287
  double m_VoxelDimensionInY;
287
288
  /// Voxel dimension in z
288
289
  double m_VoxelDimensionInZ;
289
 
    
 
290
 
290
291
  /// The coordinate of the point at which the ray enters the volume in mm.
291
292
  double m_RayStartCoordInMM[3];
292
293
  /// The coordinate of the point at which the ray exits the volume in mm.
293
294
  double m_RayEndCoordInMM[3];
294
 
    
295
 
    
 
295
 
 
296
 
296
297
  /** \brief
297
 
      Planes which define the boundary of the volume in mm 
 
298
      Planes which define the boundary of the volume in mm
298
299
      (six planes and four parameters: Ax+By+Cz+D). */
299
300
  double m_BoundingPlane[6][4];
300
301
  /// The eight corners of the volume (x,y,z coordinates for each).
301
302
  double m_BoundingCorner[8][3];
302
 
    
 
303
 
303
304
  /// The position of the ray
304
305
  double m_CurrentRayPositionInMM[3];
305
 
    
 
306
 
306
307
  /// The direction of the ray
307
 
  double m_RayDirectionInMM[3];  
308
 
    
 
308
  double m_RayDirectionInMM[3];
 
309
 
309
310
};
310
311
 
311
 
 
312
 
 
313
312
/* -----------------------------------------------------------------------
314
313
   Initialise() - Initialise the object
315
314
   ----------------------------------------------------------------------- */
316
315
 
317
316
template<class TInputImage, class TCoordRep>
318
 
void 
 
317
void
319
318
RayCastHelper<TInputImage, TCoordRep>
320
 
::Initialise(void) 
 
319
::Initialise(void)
321
320
{
322
321
  // Save the dimensions of the volume and calculate the bounding box
323
322
  this->RecordVolumeDimensions();
333
332
   ----------------------------------------------------------------------- */
334
333
 
335
334
template<class TInputImage, class TCoordRep>
336
 
void 
 
335
void
337
336
RayCastHelper<TInputImage, TCoordRep>
338
 
::RecordVolumeDimensions(void) 
 
337
::RecordVolumeDimensions(void)
339
338
{
340
339
  typename InputImageType::SpacingType spacing=this->m_Image->GetSpacing();
341
340
  SizeType dim=this->m_Image->GetLargestPossibleRegion().GetSize();
355
354
   ----------------------------------------------------------------------- */
356
355
 
357
356
template<class TInputImage, class TCoordRep>
358
 
void 
 
357
void
359
358
RayCastHelper<TInputImage, TCoordRep>
360
 
::DefineCorners(void) 
 
359
::DefineCorners(void)
361
360
{
362
361
  // Define corner positions as if at the origin
363
362
 
364
 
  m_BoundingCorner[0][0] 
365
 
    = m_BoundingCorner[1][0] 
366
 
    = m_BoundingCorner[2][0] 
367
 
    = m_BoundingCorner[3][0] = 0;
368
 
 
369
 
  m_BoundingCorner[4][0] 
370
 
    = m_BoundingCorner[5][0] 
371
 
    = m_BoundingCorner[6][0] 
372
 
    = m_BoundingCorner[7][0] 
373
 
    = m_VoxelDimensionInX*m_NumberOfVoxelsInX;
374
 
  
375
 
  m_BoundingCorner[1][1] 
376
 
    = m_BoundingCorner[3][1] 
377
 
    = m_BoundingCorner[5][1] 
378
 
    = m_BoundingCorner[7][1] 
379
 
    = m_VoxelDimensionInY*m_NumberOfVoxelsInY;
380
 
 
381
 
  m_BoundingCorner[0][1] 
382
 
    = m_BoundingCorner[2][1] 
383
 
    = m_BoundingCorner[4][1] 
384
 
    = m_BoundingCorner[6][1] = 0;
385
 
 
386
 
  m_BoundingCorner[0][2] 
387
 
    = m_BoundingCorner[1][2] 
388
 
    = m_BoundingCorner[4][2] 
389
 
    = m_BoundingCorner[5][2] 
390
 
    = m_VoxelDimensionInZ*m_NumberOfVoxelsInZ;
391
 
 
392
 
  m_BoundingCorner[2][2] 
393
 
    = m_BoundingCorner[3][2] 
394
 
    = m_BoundingCorner[6][2] 
395
 
    = m_BoundingCorner[7][2] = 0;
 
363
  m_BoundingCorner[0][0] =
 
364
    m_BoundingCorner[1][0] =
 
365
    m_BoundingCorner[2][0] =
 
366
    m_BoundingCorner[3][0] = 0;
 
367
 
 
368
  m_BoundingCorner[4][0] =
 
369
    m_BoundingCorner[5][0] =
 
370
    m_BoundingCorner[6][0] =
 
371
    m_BoundingCorner[7][0] = m_VoxelDimensionInX*m_NumberOfVoxelsInX;
 
372
 
 
373
  m_BoundingCorner[1][1] =
 
374
    m_BoundingCorner[3][1] =
 
375
    m_BoundingCorner[5][1] =
 
376
    m_BoundingCorner[7][1] = m_VoxelDimensionInY*m_NumberOfVoxelsInY;
 
377
 
 
378
  m_BoundingCorner[0][1] =
 
379
    m_BoundingCorner[2][1] =
 
380
    m_BoundingCorner[4][1] =
 
381
    m_BoundingCorner[6][1] = 0;
 
382
 
 
383
  m_BoundingCorner[0][2] =
 
384
    m_BoundingCorner[1][2] =
 
385
    m_BoundingCorner[4][2] =
 
386
    m_BoundingCorner[5][2] =
 
387
    m_VoxelDimensionInZ*m_NumberOfVoxelsInZ;
 
388
 
 
389
  m_BoundingCorner[2][2] =
 
390
    m_BoundingCorner[3][2] =
 
391
    m_BoundingCorner[6][2] =
 
392
    m_BoundingCorner[7][2] = 0;
396
393
 
397
394
}
398
395
 
399
 
 
400
 
 
401
396
/* -----------------------------------------------------------------------
402
397
   CalcPlanesAndCorners() - Calculate the planes and corners of the volume.
403
398
   ----------------------------------------------------------------------- */
404
399
 
405
400
template<class TInputImage, class TCoordRep>
406
 
void 
 
401
void
407
402
RayCastHelper<TInputImage, TCoordRep>
408
 
::CalcPlanesAndCorners(void) 
 
403
::CalcPlanesAndCorners(void)
409
404
{
410
405
  int j;
411
406
 
412
407
 
413
408
  // find the equations of the planes
414
409
 
415
 
  int c1, c2, c3;
416
 
        
417
 
  for (j=0; j<6; j++) 
 
410
  int c1=0, c2=0, c3=0;
 
411
 
 
412
  for (j=0; j<6; j++)
418
413
    {                                // loop around for planes
419
414
    switch (j)
420
415
      {                // which corners to take
421
416
      case 0:
422
 
        c1=1; c2=2; c3=3;;
 
417
        c1=1; c2=2; c3=3;
423
418
        break;
424
419
      case 1:
425
420
        c1=4; c2=5; c3=6;
438
433
        break;
439
434
      }
440
435
 
441
 
    
 
436
 
442
437
    double line1x, line1y, line1z;
443
438
    double line2x, line2y, line2z;
444
439
 
445
440
    // lines from one corner to another in x,y,z dirns
446
441
    line1x = m_BoundingCorner[c1][0] - m_BoundingCorner[c2][0];
447
442
    line2x = m_BoundingCorner[c1][0] - m_BoundingCorner[c3][0];
448
 
    
 
443
 
449
444
    line1y = m_BoundingCorner[c1][1] - m_BoundingCorner[c2][1];
450
445
    line2y = m_BoundingCorner[c1][1] - m_BoundingCorner[c3][1];
451
 
    
 
446
 
452
447
    line1z = m_BoundingCorner[c1][2] - m_BoundingCorner[c2][2];
453
448
    line2z = m_BoundingCorner[c1][2] - m_BoundingCorner[c3][2];
454
 
    
 
449
 
455
450
    double A, B, C, D;
456
 
        
 
451
 
457
452
    // take cross product
458
453
    A = line1y*line2z - line2y*line1z;
459
454
    B = line2x*line1z - line1x*line2z;
460
455
    C = line1x*line2y - line2x*line1y;
461
456
 
462
457
    // find constant
463
 
    D = -(   A*m_BoundingCorner[c1][0] 
464
 
             + B*m_BoundingCorner[c1][1] 
 
458
    D = -(   A*m_BoundingCorner[c1][0]
 
459
             + B*m_BoundingCorner[c1][1]
465
460
             + C*m_BoundingCorner[c1][2] );
466
461
 
467
462
    // initialise plane value and normalise
469
464
    m_BoundingPlane[j][1] = B/vcl_sqrt(A*A + B*B + C*C);
470
465
    m_BoundingPlane[j][2] = C/vcl_sqrt(A*A + B*B + C*C);
471
466
    m_BoundingPlane[j][3] = D/vcl_sqrt(A*A + B*B + C*C);
472
 
    
473
 
    if ( (A*A + B*B + C*C) == 0 ) 
 
467
 
 
468
    if ( (A*A + B*B + C*C) == 0 )
474
469
      {
475
470
      itk::ExceptionObject err(__FILE__, __LINE__);
476
471
      err.SetLocation( ITK_LOCATION );
488
483
   ----------------------------------------------------------------------- */
489
484
 
490
485
template<class TInputImage, class TCoordRep>
491
 
bool 
 
486
bool
492
487
RayCastHelper<TInputImage, TCoordRep>
493
 
::CalcRayIntercepts() 
 
488
::CalcRayIntercepts()
494
489
{
495
490
  double maxInterDist, interDist;
496
491
  double cornerVect[4][3];
499
494
  double ax, ay, az, bx, by, bz;
500
495
  double cubeInter[6][3];
501
496
  double denom;
502
 
        
503
 
  int i,j, k; 
 
497
 
 
498
  int i,j, k;
504
499
  int NoSides = 6;  // =6 to allow truncation: =4 to remove truncated rays
505
500
 
506
501
  // Calculate intercept of ray with planes
508
503
  double interceptx[6], intercepty[6], interceptz[6];
509
504
  double d[6];
510
505
 
511
 
  for( j=0; j<NoSides; j++) 
 
506
  for( j=0; j<NoSides; j++)
512
507
    {
513
 
      
514
 
    denom = (  m_BoundingPlane[j][0]*m_RayDirectionInMM[0] 
515
 
               + m_BoundingPlane[j][1]*m_RayDirectionInMM[1] 
 
508
 
 
509
    denom = (  m_BoundingPlane[j][0]*m_RayDirectionInMM[0]
 
510
               + m_BoundingPlane[j][1]*m_RayDirectionInMM[1]
516
511
               + m_BoundingPlane[j][2]*m_RayDirectionInMM[2]);
517
 
      
 
512
 
518
513
    if( (long)(denom*100) != 0 )
519
514
      {
520
 
      d[j] = -(   m_BoundingPlane[j][3] 
 
515
      d[j] = -(   m_BoundingPlane[j][3]
521
516
                  + m_BoundingPlane[j][0]*m_CurrentRayPositionInMM[0]
522
 
                  + m_BoundingPlane[j][1]*m_CurrentRayPositionInMM[1] 
 
517
                  + m_BoundingPlane[j][1]*m_CurrentRayPositionInMM[1]
523
518
                  + m_BoundingPlane[j][2]*m_CurrentRayPositionInMM[2] ) / denom;
524
 
      
 
519
 
525
520
      interceptx[j] = m_CurrentRayPositionInMM[0] + d[j]*m_RayDirectionInMM[0];
526
521
      intercepty[j] = m_CurrentRayPositionInMM[1] + d[j]*m_RayDirectionInMM[1];
527
522
      interceptz[j] = m_CurrentRayPositionInMM[2] + d[j]*m_RayDirectionInMM[2];
528
 
      
 
523
 
529
524
      noInterFlag[j] = 1;  //OK
530
525
      }
531
526
    else
533
528
      noInterFlag[j] = 0;  //NOT OK
534
529
      }
535
530
    }
536
 
  
 
531
 
537
532
 
538
533
  nSidesCrossed = 0;
539
 
  for( j=0; j<NoSides; j++ ) 
 
534
  for( j=0; j<NoSides; j++ )
540
535
    {
541
 
      
 
536
 
542
537
    // Work out which corners to use
543
 
      
544
 
    if( j==0 ) 
 
538
 
 
539
    if( j==0 )
545
540
      {
546
541
      c[0] = 0; c[1] = 1; c[2] = 3; c[3] = 2;
547
542
      }
548
 
    else if( j==1 ) 
 
543
    else if( j==1 )
549
544
      {
550
545
      c[0] = 4; c[1] = 5; c[2] = 7; c[3] = 6;
551
546
      }
552
 
    else if( j==2 ) 
 
547
    else if( j==2 )
553
548
      {
554
549
      c[0] = 1; c[1] = 5; c[2] = 7; c[3] = 3;
555
550
      }
556
 
    else if( j==3 ) 
 
551
    else if( j==3 )
557
552
      {
558
553
      c[0] = 0; c[1] = 2; c[2] = 6; c[3] = 4;
559
554
      }
560
 
    else if( j==4 ) 
 
555
    else if( j==4 )
561
556
      { //TOP
562
557
      c[0] = 0; c[1] = 1; c[2] = 5; c[3] = 4;
563
558
      }
564
 
    else if( j==5 ) 
 
559
    else if( j==5 )
565
560
      { //BOTTOM
566
561
      c[0] = 2; c[1] = 3; c[2] = 7; c[3] = 6;
567
562
      }
568
563
 
569
564
    // Calculate vectors from corner of ct volume to intercept.
570
 
    for( i=0; i<4; i++ ) 
 
565
    for( i=0; i<4; i++ )
571
566
      {
572
 
      if( noInterFlag[j]==1 ) 
 
567
      if( noInterFlag[j]==1 )
573
568
        {
574
569
        cornerVect[i][0] = m_BoundingCorner[c[i]][0] - interceptx[j];
575
570
        cornerVect[i][1] = m_BoundingCorner[c[i]][1] - intercepty[j];
576
571
        cornerVect[i][2] = m_BoundingCorner[c[i]][2] - interceptz[j];
577
572
        }
578
 
      else if( noInterFlag[j]==0 ) 
 
573
      else if( noInterFlag[j]==0 )
579
574
        {
580
575
        cornerVect[i][0] = 0;
581
576
        cornerVect[i][1] = 0;
585
580
      }
586
581
 
587
582
    // Do cross product with these vectors
588
 
    for( i=0; i<4; i++ ) 
 
583
    for( i=0; i<4; i++ )
589
584
      {
590
585
      if( i==3 )
591
586
        {
617
612
    // if not, then the ray went through this plane
618
613
 
619
614
    crossFlag=0;
620
 
    for( i=0; i<3; i++ ) 
 
615
    for( i=0; i<3; i++ )
621
616
      {
622
 
      if( (   cross[0][i]<=0 
623
 
              && cross[1][i]<=0 
624
 
              && cross[2][i]<=0 
 
617
      if( (   cross[0][i]<=0
 
618
              && cross[1][i]<=0
 
619
              && cross[2][i]<=0
625
620
              && cross[3][i]<=0)
626
621
 
627
 
          || (   cross[0][i]>=0 
628
 
                 && cross[1][i]>=0 
629
 
                 && cross[2][i]>=0 
 
622
          || (   cross[0][i]>=0
 
623
                 && cross[1][i]>=0
 
624
                 && cross[2][i]>=0
630
625
                 && cross[3][i]>=0) )
631
626
        {
632
627
        crossFlag++;
634
629
      }
635
630
 
636
631
 
637
 
    if( crossFlag==3 && noInterFlag[j]==1 ) 
 
632
    if( crossFlag==3 && noInterFlag[j]==1 )
638
633
      {
639
634
      cubeInter[nSidesCrossed][0] = interceptx[j];
640
635
      cubeInter[nSidesCrossed][1] = intercepty[j];
658
653
    }
659
654
 
660
655
  // If 'nSidesCrossed' is larger than 2, this means that the ray goes through
661
 
  // a corner of the volume and due to rounding errors, the ray is 
 
656
  // a corner of the volume and due to rounding errors, the ray is
662
657
  // deemed to go through more than two planes.  To obtain the correct
663
658
  // start and end positions we choose the two intercept values which
664
659
  // are furthest from each other.
665
660
 
666
 
        
667
 
  if( nSidesCrossed >= 3 ) 
 
661
 
 
662
  if( nSidesCrossed >= 3 )
668
663
    {
669
664
    maxInterDist = 0;
670
665
    for( j=0; j<nSidesCrossed-1; j++ )
673
668
        {
674
669
        interDist = 0;
675
670
        for( i=0; i<3; i++ )
676
 
          { 
 
671
          {
677
672
          interDist += (cubeInter[j][i] - cubeInter[k][i])*
678
673
            (cubeInter[j][i] - cubeInter[k][i]);
679
674
          }
684
679
          m_RayStartCoordInMM[0] = cubeInter[j][0];
685
680
          m_RayStartCoordInMM[1] = cubeInter[j][1];
686
681
          m_RayStartCoordInMM[2] = cubeInter[j][2];
687
 
          
 
682
 
688
683
          m_RayEndCoordInMM[0] = cubeInter[k][0];
689
684
          m_RayEndCoordInMM[1] = cubeInter[k][1];
690
685
          m_RayEndCoordInMM[2] = cubeInter[k][2];
710
705
   ----------------------------------------------------------------------- */
711
706
 
712
707
template<class TInputImage, class TCoordRep>
713
 
bool 
 
708
bool
714
709
RayCastHelper<TInputImage, TCoordRep>
715
 
::SetRay(OutputPointType RayPosn, DirectionType RayDirn) 
 
710
::SetRay(OutputPointType RayPosn, DirectionType RayDirn)
716
711
{
717
712
 
718
713
  // Store the position and direction of the ray
723
718
  m_NumberOfVoxelsInX = dim[0];
724
719
  m_NumberOfVoxelsInY = dim[1];
725
720
  m_NumberOfVoxelsInZ = dim[2];
726
 
                      
 
721
 
727
722
  m_VoxelDimensionInX = spacing[0];
728
723
  m_VoxelDimensionInY = spacing[1];
729
724
  m_VoxelDimensionInZ = spacing[2];
730
725
 
731
 
  m_CurrentRayPositionInMM[0] = 
732
 
    RayPosn[0] + 0.5*m_VoxelDimensionInX*(double)m_NumberOfVoxelsInX; 
733
 
 
734
 
  m_CurrentRayPositionInMM[1] = 
735
 
    RayPosn[1] + 0.5*m_VoxelDimensionInY*(double)m_NumberOfVoxelsInY; 
736
 
 
737
 
  m_CurrentRayPositionInMM[2] = 
738
 
    RayPosn[2] + 0.5*m_VoxelDimensionInZ*(double)m_NumberOfVoxelsInZ; 
739
 
                                        
 
726
  m_CurrentRayPositionInMM[0] =
 
727
    RayPosn[0] + 0.5*m_VoxelDimensionInX*(double)m_NumberOfVoxelsInX;
 
728
 
 
729
  m_CurrentRayPositionInMM[1] =
 
730
    RayPosn[1] + 0.5*m_VoxelDimensionInY*(double)m_NumberOfVoxelsInY;
 
731
 
 
732
  m_CurrentRayPositionInMM[2] =
 
733
    RayPosn[2] + 0.5*m_VoxelDimensionInZ*(double)m_NumberOfVoxelsInZ;
 
734
 
740
735
  m_RayDirectionInMM[0] = RayDirn[0];
741
736
  m_RayDirectionInMM[1] = RayDirn[1];
742
737
  m_RayDirectionInMM[2] = RayDirn[2];
745
740
 
746
741
  m_ValidRay = this->CalcRayIntercepts();
747
742
 
748
 
  if (! m_ValidRay) 
 
743
  if (! m_ValidRay)
749
744
    {
750
745
    Reset();
751
746
    return false;
784
779
   ----------------------------------------------------------------------- */
785
780
 
786
781
template<class TInputImage, class TCoordRep>
787
 
void 
 
782
void
788
783
RayCastHelper<TInputImage, TCoordRep>
789
 
::EndPointsInVoxels(void) 
 
784
::EndPointsInVoxels(void)
790
785
{
791
786
  m_RayVoxelStartPosition[0] = m_RayStartCoordInMM[0]/m_VoxelDimensionInX;
792
787
  m_RayVoxelStartPosition[1] = m_RayStartCoordInMM[1]/m_VoxelDimensionInY;
804
799
   ----------------------------------------------------------------------- */
805
800
 
806
801
template<class TInputImage, class TCoordRep>
807
 
void 
 
802
void
808
803
RayCastHelper<TInputImage, TCoordRep>
809
 
::CalcDirnVector(void) 
 
804
::CalcDirnVector(void)
810
805
{
811
806
  double xNum, yNum, zNum;
812
807
 
821
816
 
822
817
  // Iterate in X direction
823
818
 
824
 
  if( (xNum >= yNum) && (xNum >= zNum) ) 
 
819
  if( (xNum >= yNum) && (xNum >= zNum) )
825
820
    {
826
 
    if( m_RayVoxelStartPosition[0] < m_RayVoxelEndPosition[0] ) 
827
 
      { 
 
821
    if( m_RayVoxelStartPosition[0] < m_RayVoxelEndPosition[0] )
 
822
      {
828
823
      m_VoxelIncrement[0] = 1;
829
824
 
830
 
      m_VoxelIncrement[1] 
831
 
        = (m_RayVoxelStartPosition[1] 
832
 
           - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0] 
 
825
      m_VoxelIncrement[1]
 
826
        = (m_RayVoxelStartPosition[1]
 
827
           - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
833
828
                                        - m_RayVoxelEndPosition[0]);
834
 
          
835
 
      m_VoxelIncrement[2] 
836
 
        = (m_RayVoxelStartPosition[2] 
837
 
           - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0] 
 
829
 
 
830
      m_VoxelIncrement[2]
 
831
        = (m_RayVoxelStartPosition[2]
 
832
           - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
838
833
                                        - m_RayVoxelEndPosition[0]);
839
834
      }
840
 
    else 
 
835
    else
841
836
      {
842
837
      m_VoxelIncrement[0] = -1;
843
838
 
844
 
      m_VoxelIncrement[1] 
845
 
        = -(m_RayVoxelStartPosition[1] 
846
 
            - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0] 
 
839
      m_VoxelIncrement[1]
 
840
        = -(m_RayVoxelStartPosition[1]
 
841
            - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
847
842
                                         - m_RayVoxelEndPosition[0]);
848
843
 
849
 
      m_VoxelIncrement[2] 
 
844
      m_VoxelIncrement[2]
850
845
        = -(m_RayVoxelStartPosition[2]
851
 
            - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0] 
 
846
            - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
852
847
                                         - m_RayVoxelEndPosition[0]);
853
848
      }
854
849
 
855
 
    // This section is to alter the start position in order to 
 
850
    // This section is to alter the start position in order to
856
851
    // place the center of the voxels in there correct positions,
857
852
    // rather than placing them at the corner of voxels which is
858
853
    // what happens if this is not carried out.  The reason why
859
854
    // x has no -0.5 is because this is the direction we are going
860
 
    // to iterate in and therefore we wish to go from center to 
 
855
    // to iterate in and therefore we wish to go from center to
861
856
    // center rather than finding the surrounding voxels.
862
 
    
863
 
    m_RayVoxelStartPosition[1] 
864
 
      += ( (int)m_RayVoxelStartPosition[0]
865
 
           - m_RayVoxelStartPosition[0])*m_VoxelIncrement[1]*m_VoxelIncrement[0] 
 
857
 
 
858
    m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[0]
 
859
        - m_RayVoxelStartPosition[0])*m_VoxelIncrement[1]*m_VoxelIncrement[0]
866
860
      + 0.5*m_VoxelIncrement[1] - 0.5;
867
861
 
868
 
    m_RayVoxelStartPosition[2] 
869
 
      += ( (int)m_RayVoxelStartPosition[0] 
 
862
    m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[0]
870
863
           - m_RayVoxelStartPosition[0])*m_VoxelIncrement[2]*m_VoxelIncrement[0]
871
864
      + 0.5*m_VoxelIncrement[2] - 0.5;
872
865
 
873
 
    m_RayVoxelStartPosition[0] 
874
 
      = (int)m_RayVoxelStartPosition[0] + 0.5*m_VoxelIncrement[0];
 
866
    m_RayVoxelStartPosition[0] = (int)m_RayVoxelStartPosition[0] + 0.5*m_VoxelIncrement[0];
875
867
 
876
868
    m_TotalRayVoxelPlanes = (int)xNum;
877
869
 
880
872
 
881
873
  // Iterate in Y direction
882
874
 
883
 
  else if( (yNum >= xNum) && (yNum >= zNum) ) 
 
875
  else if( (yNum >= xNum) && (yNum >= zNum) )
884
876
    {
885
877
 
886
 
    if( m_RayVoxelStartPosition[1] < m_RayVoxelEndPosition[1] ) 
 
878
    if( m_RayVoxelStartPosition[1] < m_RayVoxelEndPosition[1] )
887
879
      {
888
880
      m_VoxelIncrement[1] = 1;
889
881
 
890
 
      m_VoxelIncrement[0] 
891
 
        = (m_RayVoxelStartPosition[0] 
892
 
           - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1] 
 
882
      m_VoxelIncrement[0]
 
883
        = (m_RayVoxelStartPosition[0]
 
884
           - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
893
885
                                        - m_RayVoxelEndPosition[1]);
894
886
 
895
 
      m_VoxelIncrement[2] 
896
 
        = (m_RayVoxelStartPosition[2] 
897
 
           - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1] 
 
887
      m_VoxelIncrement[2]
 
888
        = (m_RayVoxelStartPosition[2]
 
889
           - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
898
890
                                        - m_RayVoxelEndPosition[1]);
899
891
      }
900
 
    else 
 
892
    else
901
893
      {
902
894
      m_VoxelIncrement[1] = -1;
903
895
 
904
 
      m_VoxelIncrement[0] 
905
 
        = -(m_RayVoxelStartPosition[0] 
906
 
            - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1] 
 
896
      m_VoxelIncrement[0]
 
897
        = -(m_RayVoxelStartPosition[0]
 
898
            - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
907
899
                                         - m_RayVoxelEndPosition[1]);
908
900
 
909
 
      m_VoxelIncrement[2] 
910
 
        = -(m_RayVoxelStartPosition[2] 
911
 
            - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1] 
 
901
      m_VoxelIncrement[2]
 
902
        = -(m_RayVoxelStartPosition[2]
 
903
            - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
912
904
                                         - m_RayVoxelEndPosition[1]);
913
905
      }
914
906
 
915
 
    m_RayVoxelStartPosition[0] 
916
 
      += ( (int)m_RayVoxelStartPosition[1] 
917
 
           - m_RayVoxelStartPosition[1])*m_VoxelIncrement[0]*m_VoxelIncrement[1] 
 
907
    m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[1]
 
908
                                    - m_RayVoxelStartPosition[1])*m_VoxelIncrement[0]*m_VoxelIncrement[1]
918
909
      + 0.5*m_VoxelIncrement[0] - 0.5;
919
910
 
920
 
    m_RayVoxelStartPosition[2] 
921
 
      += ( (int)m_RayVoxelStartPosition[1] 
922
 
           - m_RayVoxelStartPosition[1])*m_VoxelIncrement[2]*m_VoxelIncrement[1] 
 
911
    m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[1]
 
912
                                    - m_RayVoxelStartPosition[1])*m_VoxelIncrement[2]*m_VoxelIncrement[1]
923
913
      + 0.5*m_VoxelIncrement[2] - 0.5;
924
914
 
925
 
    m_RayVoxelStartPosition[1] 
926
 
      = (int)m_RayVoxelStartPosition[1] + 0.5*m_VoxelIncrement[1];
927
 
    
 
915
    m_RayVoxelStartPosition[1] = (int)m_RayVoxelStartPosition[1] + 0.5*m_VoxelIncrement[1];
 
916
 
928
917
    m_TotalRayVoxelPlanes = (int)yNum;
929
918
 
930
919
    m_TraversalDirection = TRANSVERSE_IN_Y;
932
921
 
933
922
  // Iterate in Z direction
934
923
 
935
 
  else 
 
924
  else
936
925
    {
937
926
 
938
 
    if( m_RayVoxelStartPosition[2] < m_RayVoxelEndPosition[2] ) 
 
927
    if( m_RayVoxelStartPosition[2] < m_RayVoxelEndPosition[2] )
939
928
      {
940
929
      m_VoxelIncrement[2] = 1;
941
 
          
 
930
 
942
931
      m_VoxelIncrement[0]
943
 
        = (m_RayVoxelStartPosition[0] 
944
 
           - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2] 
 
932
        = (m_RayVoxelStartPosition[0]
 
933
           - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
945
934
                                        - m_RayVoxelEndPosition[2]);
946
935
 
947
 
      m_VoxelIncrement[1] 
948
 
        = (m_RayVoxelStartPosition[1] 
949
 
           - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2] 
 
936
      m_VoxelIncrement[1]
 
937
        = (m_RayVoxelStartPosition[1]
 
938
           - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
950
939
                                        - m_RayVoxelEndPosition[2]);
951
940
      }
952
 
    else 
 
941
    else
953
942
      {
954
943
      m_VoxelIncrement[2] = -1;
955
944
 
956
 
      m_VoxelIncrement[0] 
957
 
        = -(m_RayVoxelStartPosition[0] 
958
 
            - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2] 
 
945
      m_VoxelIncrement[0]
 
946
        = -(m_RayVoxelStartPosition[0]
 
947
            - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
959
948
                                         - m_RayVoxelEndPosition[2]);
960
949
 
961
 
      m_VoxelIncrement[1] 
962
 
        = -(m_RayVoxelStartPosition[1] 
963
 
            - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2] 
 
950
      m_VoxelIncrement[1]
 
951
        = -(m_RayVoxelStartPosition[1]
 
952
            - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
964
953
                                         - m_RayVoxelEndPosition[2]);
965
954
      }
966
955
 
967
 
    m_RayVoxelStartPosition[0] 
968
 
      += ( (int)m_RayVoxelStartPosition[2] 
969
 
           - m_RayVoxelStartPosition[2])*m_VoxelIncrement[0]*m_VoxelIncrement[2] 
 
956
    m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[2]
 
957
                                    - m_RayVoxelStartPosition[2])*m_VoxelIncrement[0]*m_VoxelIncrement[2]
970
958
      + 0.5*m_VoxelIncrement[0] - 0.5;
971
959
 
972
 
    m_RayVoxelStartPosition[1] 
973
 
      += ( (int)m_RayVoxelStartPosition[2] 
974
 
           - m_RayVoxelStartPosition[2])*m_VoxelIncrement[1]*m_VoxelIncrement[2] 
 
960
    m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[2]
 
961
                                    - m_RayVoxelStartPosition[2])*m_VoxelIncrement[1]*m_VoxelIncrement[2]
975
962
      + 0.5*m_VoxelIncrement[1] - 0.5;
976
963
 
977
 
    m_RayVoxelStartPosition[2] 
978
 
      = (int)m_RayVoxelStartPosition[2] + 0.5*m_VoxelIncrement[2];
979
 
    
 
964
    m_RayVoxelStartPosition[2] = (int)m_RayVoxelStartPosition[2] + 0.5*m_VoxelIncrement[2];
 
965
 
980
966
    m_TotalRayVoxelPlanes = (int)zNum;
981
967
 
982
968
    m_TraversalDirection = TRANSVERSE_IN_Z;
989
975
   ----------------------------------------------------------------------- */
990
976
 
991
977
template<class TInputImage, class TCoordRep>
992
 
bool 
 
978
bool
993
979
RayCastHelper<TInputImage, TCoordRep>
994
 
::AdjustRayLength(void) 
 
980
::AdjustRayLength(void)
995
981
{
996
982
  bool startOK, endOK;
997
983
 
998
984
  int Istart[3];
999
985
  int Idirn[3];
1000
986
 
1001
 
  if (m_TraversalDirection == TRANSVERSE_IN_X) 
1002
 
    {
1003
 
    Idirn[0] = 0;    
1004
 
    Idirn[1] = 1;    
1005
 
    Idirn[2] = 1;
1006
 
    }
1007
 
  else if (m_TraversalDirection == TRANSVERSE_IN_Y) 
1008
 
    {
1009
 
    Idirn[0] = 1;    
1010
 
    Idirn[1] = 0;    
1011
 
    Idirn[2] = 1;
1012
 
    }
1013
 
  else if (m_TraversalDirection == TRANSVERSE_IN_Z) 
1014
 
    {
1015
 
    Idirn[0] = 1;    
1016
 
    Idirn[1] = 1;    
 
987
  if (m_TraversalDirection == TRANSVERSE_IN_X)
 
988
    {
 
989
    Idirn[0] = 0;
 
990
    Idirn[1] = 1;
 
991
    Idirn[2] = 1;
 
992
    }
 
993
  else if (m_TraversalDirection == TRANSVERSE_IN_Y)
 
994
    {
 
995
    Idirn[0] = 1;
 
996
    Idirn[1] = 0;
 
997
    Idirn[2] = 1;
 
998
    }
 
999
  else if (m_TraversalDirection == TRANSVERSE_IN_Z)
 
1000
    {
 
1001
    Idirn[0] = 1;
 
1002
    Idirn[1] = 1;
1017
1003
    Idirn[2] = 0;
1018
1004
    }
1019
 
  else 
 
1005
  else
1020
1006
    {
1021
1007
    itk::ExceptionObject err(__FILE__, __LINE__);
1022
1008
    err.SetLocation( ITK_LOCATION );
1027
1013
    }
1028
1014
 
1029
1015
 
1030
 
  do 
 
1016
  do
1031
1017
    {
1032
1018
 
1033
1019
    startOK = false;
1039
1025
 
1040
1026
    if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
1041
1027
        (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
1042
 
        (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) ) 
 
1028
        (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) )
1043
1029
      {
1044
1030
 
1045
1031
      startOK = true;
1046
1032
      }
1047
 
    else 
 
1033
    else
1048
1034
      {
1049
1035
      m_RayVoxelStartPosition[0] += m_VoxelIncrement[0];
1050
1036
      m_RayVoxelStartPosition[1] += m_VoxelIncrement[1];
1053
1039
      m_TotalRayVoxelPlanes--;
1054
1040
      }
1055
1041
 
1056
 
    Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0] 
 
1042
    Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0]
1057
1043
                            + m_TotalRayVoxelPlanes*m_VoxelIncrement[0]);
1058
1044
 
1059
 
    Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1] 
 
1045
    Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1]
1060
1046
                            + m_TotalRayVoxelPlanes*m_VoxelIncrement[1]);
1061
1047
 
1062
 
    Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2] 
 
1048
    Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2]
1063
1049
                            + m_TotalRayVoxelPlanes*m_VoxelIncrement[2]);
1064
 
  
 
1050
 
1065
1051
    if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
1066
1052
        (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
1067
 
        (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) ) 
 
1053
        (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) )
1068
1054
      {
1069
1055
 
1070
1056
      endOK = true;
1071
1057
      }
1072
 
    else 
 
1058
    else
1073
1059
      {
1074
1060
      m_TotalRayVoxelPlanes--;
1075
 
      }        
 
1061
      }
1076
1062
 
1077
1063
    } while ( (! (startOK && endOK)) && (m_TotalRayVoxelPlanes > 1) );
1078
1064
 
1086
1072
   ----------------------------------------------------------------------- */
1087
1073
 
1088
1074
template<class TInputImage, class TCoordRep>
1089
 
void 
 
1075
void
1090
1076
RayCastHelper<TInputImage, TCoordRep>
1091
 
::Reset(void) 
 
1077
::Reset(void)
1092
1078
{
1093
1079
  int i;
1094
1080
 
1095
1081
  m_NumVoxelPlanesTraversed = -1;
1096
 
    
 
1082
 
1097
1083
  // If this is a valid ray...
1098
1084
 
1099
 
  if (m_ValidRay) 
 
1085
  if (m_ValidRay)
1100
1086
    {
1101
1087
    for (i=0; i<3; i++)
1102
1088
      {
1106
1092
    }
1107
1093
 
1108
1094
  // otherwise set parameters to zero
1109
 
  
1110
 
  else 
 
1095
 
 
1096
  else
1111
1097
    {
1112
1098
    for (i=0; i<3; i++)
1113
1099
      {
1116
1102
    for (i=0; i<3; i++)
1117
1103
      {
1118
1104
      m_RayVoxelEndPosition[i] = 0.;
1119
 
      }    
 
1105
      }
1120
1106
    for (i=0; i<3; i++)
1121
1107
      {
1122
1108
      m_VoxelIncrement[i] = 0.;
1123
 
      }      
 
1109
      }
1124
1110
    m_TraversalDirection = UNDEFINED_DIRECTION;
1125
 
    
 
1111
 
1126
1112
    m_TotalRayVoxelPlanes = 0;
1127
 
    
 
1113
 
1128
1114
    for (i=0; i<4; i++)
1129
1115
      {
1130
1116
      m_RayIntersectionVoxels[i] = 0;
1142
1128
   ----------------------------------------------------------------------- */
1143
1129
 
1144
1130
template<class TInputImage, class TCoordRep>
1145
 
void 
 
1131
void
1146
1132
RayCastHelper<TInputImage, TCoordRep>
1147
 
::InitialiseVoxelPointers(void) 
 
1133
::InitialiseVoxelPointers(void)
1148
1134
{
1149
1135
  IndexType index;
1150
1136
 
1158
1144
  m_RayIntersectionVoxelIndex[1] = Iy;
1159
1145
  m_RayIntersectionVoxelIndex[2] = Iz;
1160
1146
 
1161
 
  switch( m_TraversalDirection ) 
1162
 
    
1163
 
    {
1164
 
    case TRANSVERSE_IN_X: 
1165
 
    {
1166
 
 
1167
 
    if( (Ix >= 0) && (Ix     < m_NumberOfVoxelsInX) && 
1168
 
        (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) && 
1169
 
        (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ)) 
1170
 
      {
1171
 
 
1172
 
      index[0]=Ix; index[1]=Iy; index[2]=Iz; 
1173
 
      m_RayIntersectionVoxels[0]
1174
 
        = this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index);
1175
 
 
1176
 
      index[0]=Ix; index[1]=Iy+1; index[2]=Iz; 
1177
 
      m_RayIntersectionVoxels[1] 
1178
 
        = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1179
 
 
1180
 
      index[0]=Ix; index[1]=Iy; index[2]=Iz+1; 
1181
 
      m_RayIntersectionVoxels[2] 
1182
 
        = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1183
 
 
1184
 
      index[0]=Ix; index[1]=Iy+1; index[2]=Iz+1; 
1185
 
      m_RayIntersectionVoxels[3] 
1186
 
        = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1187
 
      }
1188
 
    else
1189
 
      m_RayIntersectionVoxels[0] 
1190
 
        = m_RayIntersectionVoxels[1] 
1191
 
        = m_RayIntersectionVoxels[2] 
1192
 
        = m_RayIntersectionVoxels[3] = NULL;
1193
 
 
1194
 
    break;
1195
 
    }
1196
 
 
1197
 
    case TRANSVERSE_IN_Y: 
1198
 
    {
1199
 
 
1200
 
    if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX) && 
1201
 
        (Iy >= 0) && (Iy     < m_NumberOfVoxelsInY) && 
1202
 
        (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ)) 
1203
 
      {
1204
 
 
1205
 
      index[0]=Ix; index[1]=Iy; index[2]=Iz; 
1206
 
      m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
1207
 
                                     + this->m_Image->ComputeOffset(index) );
1208
 
 
1209
 
      index[0]=Ix+1; index[1]=Iy; index[2]=Iz; 
1210
 
      m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
1211
 
                                     + this->m_Image->ComputeOffset(index) );
1212
 
 
1213
 
      index[0]=Ix; index[1]=Iy; index[2]=Iz+1; 
1214
 
      m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
1215
 
                                     + this->m_Image->ComputeOffset(index) );
1216
 
 
1217
 
      index[0]=Ix+1; index[1]=Iy; index[2]=Iz+1; 
1218
 
      m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
1219
 
                                     + this->m_Image->ComputeOffset(index) );
1220
 
      }
1221
 
    else
1222
 
      m_RayIntersectionVoxels[0] 
1223
 
        = m_RayIntersectionVoxels[1] 
1224
 
        = m_RayIntersectionVoxels[2] 
1225
 
        = m_RayIntersectionVoxels[3] = NULL;        
1226
 
 
1227
 
    break;
1228
 
    }
1229
 
 
1230
 
    case TRANSVERSE_IN_Z: 
1231
 
    {
1232
 
 
1233
 
    if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX)   && 
1234
 
        (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) && 
1235
 
        (Iz >= 0) && (Iz     < m_NumberOfVoxelsInZ)) 
1236
 
      {
1237
 
 
1238
 
      index[0]=Ix; index[1]=Iy; index[2]=Iz; 
1239
 
      m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
1240
 
                                     + this->m_Image->ComputeOffset(index) );
1241
 
 
1242
 
      index[0]=Ix+1; index[1]=Iy; index[2]=Iz; 
1243
 
      m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
1244
 
                                     + this->m_Image->ComputeOffset(index) );
1245
 
 
1246
 
      index[0]=Ix; index[1]=Iy+1; index[2]=Iz; 
1247
 
      m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
1248
 
                                     + this->m_Image->ComputeOffset(index) );
1249
 
 
1250
 
      index[0]=Ix+1; index[1]=Iy+1; index[2]=Iz; 
1251
 
      m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
1252
 
                                     + this->m_Image->ComputeOffset(index) );
1253
 
 
1254
 
      }
1255
 
    else
1256
 
      m_RayIntersectionVoxels[0] 
1257
 
        = m_RayIntersectionVoxels[1]
1258
 
        = m_RayIntersectionVoxels[2]
1259
 
        = m_RayIntersectionVoxels[3] = NULL;
1260
 
 
1261
 
    break;
1262
 
    }
1263
 
 
1264
 
    default: 
1265
 
    {
1266
 
    itk::ExceptionObject err(__FILE__, __LINE__);
1267
 
    err.SetLocation( ITK_LOCATION );
1268
 
    err.SetDescription( "The ray traversal direction is unset "
1269
 
                        "- InitialiseVoxelPointers().");
1270
 
    throw err;
1271
 
    return;
1272
 
    }
1273
 
    }  
 
1147
  switch( m_TraversalDirection )
 
1148
    {
 
1149
    case TRANSVERSE_IN_X:
 
1150
      {
 
1151
 
 
1152
      if( (Ix >= 0) && (Ix     < m_NumberOfVoxelsInX) &&
 
1153
          (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
 
1154
          (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ))
 
1155
        {
 
1156
        
 
1157
        index[0]=Ix; index[1]=Iy; index[2]=Iz;
 
1158
        m_RayIntersectionVoxels[0]
 
1159
          = this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index);
 
1160
        
 
1161
        index[0]=Ix; index[1]=Iy+1; index[2]=Iz;
 
1162
        m_RayIntersectionVoxels[1]
 
1163
          = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
 
1164
        
 
1165
        index[0]=Ix; index[1]=Iy; index[2]=Iz+1;
 
1166
        m_RayIntersectionVoxels[2]
 
1167
          = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
 
1168
        
 
1169
        index[0]=Ix; index[1]=Iy+1; index[2]=Iz+1;
 
1170
        m_RayIntersectionVoxels[3]
 
1171
          = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
 
1172
        }
 
1173
      else
 
1174
        {
 
1175
        m_RayIntersectionVoxels[0] =
 
1176
          m_RayIntersectionVoxels[1] =
 
1177
          m_RayIntersectionVoxels[2] =
 
1178
          m_RayIntersectionVoxels[3] = NULL;
 
1179
        }
 
1180
      break;
 
1181
      }
 
1182
 
 
1183
    case TRANSVERSE_IN_Y:
 
1184
      {
 
1185
 
 
1186
      if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX) &&
 
1187
          (Iy >= 0) && (Iy     < m_NumberOfVoxelsInY) &&
 
1188
          (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ))
 
1189
        {
 
1190
        
 
1191
        index[0]=Ix; index[1]=Iy; index[2]=Iz;
 
1192
        m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
 
1193
                                       + this->m_Image->ComputeOffset(index) );
 
1194
        
 
1195
        index[0]=Ix+1; index[1]=Iy; index[2]=Iz;
 
1196
        m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
 
1197
                                       + this->m_Image->ComputeOffset(index) );
 
1198
        
 
1199
        index[0]=Ix; index[1]=Iy; index[2]=Iz+1;
 
1200
        m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
 
1201
                                       + this->m_Image->ComputeOffset(index) );
 
1202
        
 
1203
        index[0]=Ix+1; index[1]=Iy; index[2]=Iz+1;
 
1204
        m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
 
1205
                                       + this->m_Image->ComputeOffset(index) );
 
1206
        }
 
1207
      else
 
1208
        {
 
1209
        m_RayIntersectionVoxels[0]
 
1210
        = m_RayIntersectionVoxels[1]
 
1211
        = m_RayIntersectionVoxels[2]
 
1212
        = m_RayIntersectionVoxels[3] = NULL;
 
1213
        }
 
1214
      break;
 
1215
      }
 
1216
 
 
1217
    case TRANSVERSE_IN_Z:
 
1218
      {
 
1219
 
 
1220
      if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX)   &&
 
1221
          (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
 
1222
          (Iz >= 0) && (Iz     < m_NumberOfVoxelsInZ))
 
1223
        {
 
1224
        
 
1225
        index[0]=Ix; index[1]=Iy; index[2]=Iz;
 
1226
        m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
 
1227
                                       + this->m_Image->ComputeOffset(index) );
 
1228
        
 
1229
        index[0]=Ix+1; index[1]=Iy; index[2]=Iz;
 
1230
        m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
 
1231
                                       + this->m_Image->ComputeOffset(index) );
 
1232
        
 
1233
        index[0]=Ix; index[1]=Iy+1; index[2]=Iz;
 
1234
        m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
 
1235
                                       + this->m_Image->ComputeOffset(index) );
 
1236
        
 
1237
        index[0]=Ix+1; index[1]=Iy+1; index[2]=Iz;
 
1238
        m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
 
1239
                                       + this->m_Image->ComputeOffset(index) );
 
1240
        
 
1241
        }
 
1242
      else
 
1243
        {
 
1244
        m_RayIntersectionVoxels[0]
 
1245
        = m_RayIntersectionVoxels[1]
 
1246
        = m_RayIntersectionVoxels[2]
 
1247
        = m_RayIntersectionVoxels[3] = NULL;
 
1248
        }
 
1249
      break;
 
1250
      }
 
1251
 
 
1252
    default:
 
1253
      {
 
1254
      itk::ExceptionObject err(__FILE__, __LINE__);
 
1255
      err.SetLocation( ITK_LOCATION );
 
1256
      err.SetDescription( "The ray traversal direction is unset "
 
1257
                          "- InitialiseVoxelPointers().");
 
1258
      throw err;
 
1259
      return;
 
1260
      }
 
1261
    }
1274
1262
}
1275
1263
 
1276
 
 
1277
 
 
1278
 
 
1279
1264
/* -----------------------------------------------------------------------
1280
1265
   IncrementVoxelPointers() - Increment the voxel pointers
1281
1266
   ----------------------------------------------------------------------- */
1282
1267
 
1283
1268
template<class TInputImage, class TCoordRep>
1284
 
void 
 
1269
void
1285
1270
RayCastHelper<TInputImage, TCoordRep>
1286
 
::IncrementVoxelPointers(void) 
 
1271
::IncrementVoxelPointers(void)
1287
1272
{
1288
1273
  double xBefore = m_Position3Dvox[0];
1289
1274
  double yBefore = m_Position3Dvox[1];
1301
1286
  m_RayIntersectionVoxelIndex[1] += dy;
1302
1287
  m_RayIntersectionVoxelIndex[2] += dz;
1303
1288
 
1304
 
  int totalRayVoxelPlanes 
 
1289
  int totalRayVoxelPlanes
1305
1290
    = dx + dy*m_NumberOfVoxelsInX + dz*m_NumberOfVoxelsInX*m_NumberOfVoxelsInY;
1306
1291
 
1307
1292
  m_RayIntersectionVoxels[0] += totalRayVoxelPlanes;
1316
1301
   ----------------------------------------------------------------------- */
1317
1302
 
1318
1303
template<class TInputImage, class TCoordRep>
1319
 
double 
 
1304
double
1320
1305
RayCastHelper<TInputImage, TCoordRep>
1321
1306
::GetCurrentIntensity(void) const
1322
1307
{
1331
1316
  b = (double) (*m_RayIntersectionVoxels[1] - a);
1332
1317
  c = (double) (*m_RayIntersectionVoxels[2] - a);
1333
1318
  d = (double) (*m_RayIntersectionVoxels[3] - a - b - c);
1334
 
  
1335
 
  switch( m_TraversalDirection ) 
1336
 
    {
1337
 
    case TRANSVERSE_IN_X: {
1338
 
 
1339
 
    y = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
1340
 
    z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
1341
 
    break;                                            
1342
 
    }                                                    
1343
 
                                                    
1344
 
    case TRANSVERSE_IN_Y: {                            
1345
 
                                                          
1346
 
    y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
1347
 
    z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
1348
 
    break;                                            
1349
 
    }                                                    
1350
 
                                                        
1351
 
    case TRANSVERSE_IN_Z: {                            
1352
 
                                                          
1353
 
    y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
1354
 
    z = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
1355
 
    break;
1356
 
    }
1357
 
 
1358
 
    default: 
1359
 
    {
1360
 
    itk::ExceptionObject err(__FILE__, __LINE__);
1361
 
    err.SetLocation( ITK_LOCATION );
1362
 
    err.SetDescription( "The ray traversal direction is unset "
1363
 
                        "- GetCurrentIntensity().");
1364
 
    throw err;
1365
 
    return 0;
1366
 
    }
1367
 
    }
1368
 
  
 
1319
 
 
1320
  switch( m_TraversalDirection )
 
1321
    {
 
1322
    case TRANSVERSE_IN_X:
 
1323
      {
 
1324
      y = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
 
1325
      z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
 
1326
      break;
 
1327
      }
 
1328
    case TRANSVERSE_IN_Y:
 
1329
      {
 
1330
      y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
 
1331
      z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
 
1332
      break;
 
1333
      }
 
1334
    case TRANSVERSE_IN_Z:
 
1335
      {
 
1336
      y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
 
1337
      z = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
 
1338
      break;
 
1339
      }
 
1340
    default:
 
1341
      {
 
1342
      itk::ExceptionObject err(__FILE__, __LINE__);
 
1343
      err.SetLocation( ITK_LOCATION );
 
1344
      err.SetDescription( "The ray traversal direction is unset "
 
1345
                          "- GetCurrentIntensity().");
 
1346
      throw err;
 
1347
      return 0;
 
1348
      }
 
1349
    }
 
1350
 
1369
1351
  return a + b*y + c*z + d*y*z;
1370
1352
}
1371
1353
 
1372
 
 
1373
 
 
1374
 
 
1375
1354
/* -----------------------------------------------------------------------
1376
1355
   IncrementIntensities() - Increment the intensities of the current ray point
1377
1356
   ----------------------------------------------------------------------- */
1378
1357
 
1379
1358
template<class TInputImage, class TCoordRep>
1380
 
void 
 
1359
void
1381
1360
RayCastHelper<TInputImage, TCoordRep>
1382
1361
::IncrementIntensities(double increment)
1383
1362
{
1401
1380
   ----------------------------------------------------------------------- */
1402
1381
 
1403
1382
template<class TInputImage, class TCoordRep>
1404
 
bool 
 
1383
bool
1405
1384
RayCastHelper<TInputImage, TCoordRep>
1406
 
::IntegrateAboveThreshold(double &integral, double threshold) 
 
1385
::IntegrateAboveThreshold(double &integral, double threshold)
1407
1386
{
1408
1387
  double intensity;
1409
1388
//  double posn3D_x, posn3D_y, posn3D_z;
1411
1390
  integral = 0.;
1412
1391
 
1413
1392
  // Check if this is a valid ray
1414
 
  
 
1393
 
1415
1394
  if (! m_ValidRay)
1416
1395
    {
1417
1396
    return false;
1418
1397
    }
1419
 
  /* Step along the ray as quickly as possible 
 
1398
  /* Step along the ray as quickly as possible
1420
1399
     integrating the interpolated intensities. */
1421
1400
 
1422
 
  for (m_NumVoxelPlanesTraversed=0; 
1423
 
       m_NumVoxelPlanesTraversed<m_TotalRayVoxelPlanes; 
1424
 
       m_NumVoxelPlanesTraversed++) 
 
1401
  for (m_NumVoxelPlanesTraversed=0;
 
1402
       m_NumVoxelPlanesTraversed<m_TotalRayVoxelPlanes;
 
1403
       m_NumVoxelPlanesTraversed++)
1425
1404
    {
1426
1405
    intensity = this->GetCurrentIntensity();
1427
1406
 
1441
1420
  return true;
1442
1421
}
1443
1422
 
1444
 
 
1445
 
 
1446
1423
/* -----------------------------------------------------------------------
1447
1424
   ZeroState() - Set the default (zero) state of the object
1448
1425
   ----------------------------------------------------------------------- */
1449
1426
 
1450
1427
template<class TInputImage, class TCoordRep>
1451
 
void  
 
1428
void
1452
1429
RayCastHelper<TInputImage, TCoordRep>
1453
 
::ZeroState()  
 
1430
::ZeroState()
1454
1431
{
1455
1432
  int i;
1456
1433
 
1462
1439
 
1463
1440
  m_VoxelDimensionInX = 0;
1464
1441
  m_VoxelDimensionInY = 0;
1465
 
  m_VoxelDimensionInZ = 0;  
 
1442
  m_VoxelDimensionInZ = 0;
1466
1443
 
1467
1444
  for (i=0; i<3; i++)
1468
1445
    {
1469
 
    m_CurrentRayPositionInMM[i] = 0.; 
 
1446
    m_CurrentRayPositionInMM[i] = 0.;
1470
1447
    }
1471
1448
  for (i=0; i<3; i++)
1472
1449
    {
1488
1465
 
1489
1466
  m_TotalRayVoxelPlanes = 0;
1490
1467
  m_NumVoxelPlanesTraversed = -1;
1491
 
  
 
1468
 
1492
1469
  for (i=0; i<4; i++)
1493
1470
    {
1494
1471
    m_RayIntersectionVoxels[i] = 0;
1513
1490
 *
1514
1491
 **************************************************************************/
1515
1492
 
1516
 
 
1517
 
 
1518
1493
/* -----------------------------------------------------------------------
1519
1494
   Constructor
1520
1495
   ----------------------------------------------------------------------- */
1521
1496
 
1522
1497
template<class TInputImage, class TCoordRep>
1523
1498
RayCastInterpolateImageFunction< TInputImage, TCoordRep >
1524
 
::RayCastInterpolateImageFunction() 
 
1499
::RayCastInterpolateImageFunction()
1525
1500
{
1526
1501
  m_Threshold = 0.;
1527
1502
 
1546
1521
  os << indent << "FocalPoint: " << m_FocalPoint << std::endl;
1547
1522
  os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
1548
1523
  os << indent << "Interpolator: " << m_Interpolator.GetPointer() << std::endl;
1549
 
  
 
1524
 
1550
1525
}
1551
1526
 
1552
1527
/* -----------------------------------------------------------------------
1561
1536
{
1562
1537
  double integral = 0;
1563
1538
 
1564
 
  OutputPointType transformedFocalPoint 
 
1539
  OutputPointType transformedFocalPoint
1565
1540
    = m_Transform->TransformPoint( m_FocalPoint );
1566
1541
 
1567
1542
  DirectionType direction = transformedFocalPoint - point;
1585
1560
{
1586
1561
  OutputPointType point;
1587
1562
  this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
1588
 
  
 
1563
 
1589
1564
  return this->Evaluate( point );
1590
1565
}
1591
1566