← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Code/Common/itkImageBase.h

  • Committer: Bazaar Package Importer
  • Author(s): Steve M. Robbins
  • Date: 2008-12-19 20:16:49 UTC
  • mfrom: (1.2.1 upstream) (4.1.1 sid)
  • Revision ID: james.westby@ubuntu.com-20081219201649-drt97guwl2ryt0cn

* New upstream version.
  - patches/nifti-versioning.patch: Remove.  Applied upstream.
  - control:
  - rules: Update version numbers, package names.

* control: Build-depend on uuid-dev (gdcm uses it).

* copyright: Update download URL.

* rules: Adhere to parallel=N in DEB_BUILD_OPTIONS by setting MAKEFLAGS.

* compat: Set to 7.
* control: Update build-dep on debhelper to version >= 7.

* CMakeCache.txt.debian: Set CMAKE_BUILD_TYPE to "RELEASE" so that we
  build with -O3 (not -O2), necessary to optimize the templated code.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Code/Common/itkImageBase.h
3
3
  Program:   Insight Segmentation & Registration Toolkit
4
4
  Module:    $RCSfile: itkImageBase.h,v $
5
5
  Language:  C++
6
 
  Date:      $Date: 2007-10-21 18:03:05 $
7
 
  Version:   $Revision: 1.69 $
 
6
  Date:      $Date: 2008-11-01 15:28:00 $
 
7
  Version:   $Revision: 1.75 $
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.
12
12
  Portions of this code are covered under the VTK copyright.
13
13
  See VTKCopyright.txt or http://www.kitware.com/VTKCopyright.htm for details.
14
14
 
15
 
     This software is distributed WITHOUT ANY WARRANTY; without even 
16
 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
 
15
     This software is distributed WITHOUT ANY WARRANTY; without even
 
16
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
17
17
     PURPOSE.  See the above copyright notices for more information.
18
18
 
19
19
=========================================================================*/
31
31
#include "itkFixedArray.h"
32
32
#include "itkPoint.h"
33
33
#include "itkMatrix.h"
 
34
#include "itkContinuousIndex.h"
34
35
#include "itkImageHelper.h"
35
36
#include <vnl/vnl_matrix_fixed.txx>
36
37
 
37
38
#include "itkImageRegion.h"
38
39
 
 
40
#ifdef ITK_USE_TEMPLATE_META_PROGRAMMING_LOOP_UNROLLING
 
41
#include "itkImageTransformHelper.h"
 
42
#endif
 
43
 
39
44
namespace itk
40
45
{
41
46
 
49
54
struct GetImageDimension
50
55
{
51
56
  itkStaticConstMacro(ImageDimension, unsigned int, TImage::ImageDimension);
52
 
}; 
53
 
  
 
57
};
 
58
 
54
59
/** \class ImageBase
55
60
 * \brief Base class for templated image classes.
56
61
 *
62
67
 * subclasses of ImageBase, namely Image and ImageAdaptor.
63
68
 *
64
69
 * There are three sets of meta-data describing an image. These are "Region"
65
 
 * objects that define a portion of an image via a starting index for the 
66
 
 * image array and a size. The ivar LargestPossibleRegion defines the size 
 
70
 * objects that define a portion of an image via a starting index for the
 
71
 * image array and a size. The ivar LargestPossibleRegion defines the size
67
72
 * and starting index of the image dataset. The entire image dataset, however,
68
73
 * may not be resident in memory. The region of the image that is resident in
69
74
 * memory is defined by the "BufferedRegion". The Buffer is a contiguous block
70
75
 * of memory.  The third set of meta-data defines a region of interest, called
71
76
 * the "RequestedRegion". The RequestedRegion is used by the pipeline
72
 
 * execution model to define what a filter is requested to produce. 
 
77
 * execution model to define what a filter is requested to produce.
73
78
 *
74
79
 * [RegionIndex, RegionSize] C [BufferIndex, BufferSize]
75
80
 *                           C [ImageIndex, ImageSize]
83
88
{
84
89
public:
85
90
  /** Standard typedefs. */
86
 
  typedef ImageBase           Self;
87
 
  typedef DataObject  Superclass;
88
 
  typedef SmartPointer<Self>  Pointer;
89
 
  typedef SmartPointer<const Self>  ConstPointer;
90
 
  
 
91
  typedef ImageBase                     Self;
 
92
  typedef DataObject                    Superclass;
 
93
  typedef SmartPointer<Self>            Pointer;
 
94
  typedef SmartPointer<const Self>      ConstPointer;
 
95
 
91
96
  /** Method for creation through the object factory. */
92
97
  itkNewMacro(Self);
93
98
 
101
106
  itkStaticConstMacro(ImageDimension, unsigned int, VImageDimension );
102
107
 
103
108
  /** Index typedef support. An index is used to access pixel values. */
104
 
  typedef Index<VImageDimension>  IndexType;
105
 
  typedef typename IndexType::IndexValueType  IndexValueType;
106
 
  
 
109
  typedef Index<VImageDimension>                  IndexType;
 
110
  typedef typename IndexType::IndexValueType      IndexValueType;
 
111
 
107
112
  /** Offset typedef support. An offset represent relative position
108
113
   * between indices. */
109
 
  typedef Offset<VImageDimension>  OffsetType;
110
 
  typedef typename OffsetType::OffsetValueType OffsetValueType;
 
114
  typedef Offset<VImageDimension>                 OffsetType;
 
115
  typedef typename OffsetType::OffsetValueType    OffsetValueType;
111
116
 
112
117
  /** Size typedef support. A size is used to define region bounds. */
113
 
  typedef Size<VImageDimension>  SizeType;
114
 
  typedef typename SizeType::SizeValueType SizeValueType;
115
 
    
 
118
  typedef Size<VImageDimension>                   SizeType;
 
119
  typedef typename SizeType::SizeValueType        SizeValueType;
 
120
 
116
121
  /** Region typedef support. A region is used to specify a subset of an image. */
117
 
  typedef ImageRegion<VImageDimension>  RegionType;
 
122
  typedef ImageRegion<VImageDimension>            RegionType;
118
123
 
119
124
  /** Spacing typedef support.  Spacing holds the size of a pixel.  The
120
125
   * spacing is the geometric distance between image samples. ITK only
121
126
   * supports positive spacing value: negative values may cause
122
 
   * undesirable results.*/
123
 
  typedef Vector<double, VImageDimension> SpacingType;
 
127
   * undesirable results.  */
 
128
  typedef Vector<double, VImageDimension>         SpacingType;
124
129
 
125
130
  /** Origin typedef support.  The origin is the geometric coordinates
126
131
   * of the index (0,0). */
135
140
  void Initialize();
136
141
 
137
142
  /** Image dimension. The dimension of an image is fixed at construction. */
138
 
  static unsigned int GetImageDimension() 
 
143
  static unsigned int GetImageDimension()
139
144
    { return VImageDimension; }
140
145
 
141
146
  /** Set the origin of the image. The origin is the geometric
148
153
 
149
154
  /** Set the direction cosines of the image. The direction cosines
150
155
   * are vectors that point from one pixel to the next.
151
 
   * 
 
156
   *
152
157
   * One row of the matrix indicates the direction cosines of the unit vector
153
158
   * that is parallel to the lines of the image grid corresponding to that
154
 
   * dimension. For example, and image with Direction matrix 
 
159
   * dimension. For example, and image with Direction matrix
155
160
   *
156
161
   *    0.866   0.500
157
162
   *   -0.500   0.866
179
184
   * For ImageBase and Image, the default direction is identity. */
180
185
  itkGetConstReferenceMacro(Direction, DirectionType);
181
186
 
182
 
  /** Set the spacing (size of a pixel) of the image. The
183
 
   * spacing is the geometric distance between image samples.
184
 
   * It is stored internally as double, but may be set from
185
 
   * float. \sa GetSpacing() */
186
 
  itkSetMacro(Spacing, SpacingType);
187
 
  virtual void SetSpacing( const double spacing[VImageDimension] );
188
 
  virtual void SetSpacing( const float spacing[VImageDimension] );
189
 
 
190
187
  /** Get the spacing (size of a pixel) `of the image. The
191
188
   * spacing is the geometric distance between image samples.
192
189
   * The value returned is a pointer to a double array.
195
192
 
196
193
  /** Get the origin of the image. The origin is the geometric
197
194
   * coordinates of the index (0,0).  The value returned is a pointer
198
 
   * to a double array.  For ImageBase and Image, the default origin is 
 
195
   * to a double array.  For ImageBase and Image, the default origin is
199
196
   * 0. */
200
197
  itkGetConstReferenceMacro(Origin, PointType);
201
198
 
217
214
    { return m_LargestPossibleRegion;};
218
215
 
219
216
  /** Set the region object that defines the size and starting index
220
 
   * of the region of the image currently loaded in memory. 
 
217
   * of the region of the image currently loaded in memory.
221
218
   * \sa ImageRegion, SetLargestPossibleRegion(), SetRequestedRegion() */
222
219
  virtual void SetBufferedRegion(const RegionType &region);
223
220
 
224
221
  /** Get the region object that defines the size and starting index
225
 
   * of the region of the image currently loaded in memory. 
 
222
   * of the region of the image currently loaded in memory.
226
223
   * \sa ImageRegion, SetLargestPossibleRegion(), SetRequestedRegion() */
227
224
  virtual const RegionType& GetBufferedRegion() const
228
225
  { return m_BufferedRegion;};
229
 
  
 
226
 
230
227
  /** Set the region object that defines the size and starting index
231
228
   * for the region of the image requested (i.e., the region of the
232
229
   * image to be operated on by a filter). Setting the RequestedRegion
237
234
  virtual void SetRequestedRegion(const RegionType &region);
238
235
 
239
236
  /** Set the requested region from this data object to match the requested
240
 
   * region of the data object passed in as a parameter.  This method 
 
237
   * region of the data object passed in as a parameter.  This method
241
238
   * implements the API from DataObject. The data object parameter must be
242
239
   * castable to an ImageBase. Setting the RequestedRegion does not cause
243
240
   * the object to be modified. This method is called internally by
263
260
   * some data accessing algorithms. The entries in the offset table
264
261
   * are only valid after the BufferedRegion is set. */
265
262
  const OffsetValueType *GetOffsetTable() const { return m_OffsetTable; };
266
 
  
 
263
 
267
264
  /** Compute an offset from the beginning of the buffer for a pixel
268
265
   * at the specified index. The index is not checked as to whether it
269
266
   * is inside the current buffer, so the computed offset could
272
269
   * prior to calling ComputeOffset. */
273
270
 
274
271
 
275
 
#if 1
 
272
#ifdef ITK_USE_TEMPLATE_META_PROGRAMMING_LOOP_UNROLLING
276
273
  inline OffsetValueType ComputeOffset(const IndexType &ind) const
277
 
  {
 
274
    {
278
275
    OffsetValueType offset = 0;
279
276
    ImageHelper<VImageDimension,VImageDimension>::ComputeOffset(this->GetBufferedRegion().GetIndex(),
280
277
                                                                ind,
281
278
                                                                m_OffsetTable,
282
279
                                                                offset);
283
280
    return offset;
284
 
  }
 
281
    }
285
282
#else
286
283
  OffsetValueType ComputeOffset(const IndexType &ind) const
287
 
  {
 
284
    {
288
285
    // need to add bounds checking for the region/buffer?
289
286
    OffsetValueType offset=0;
290
287
    const IndexType &bufferedRegionIndex = this->GetBufferedRegion().GetIndex();
291
 
  
 
288
 
292
289
    // data is arranged as [][][][slice][row][col]
293
290
    // with Index[0] = col, Index[1] = row, Index[2] = slice
294
291
    for (int i=VImageDimension-1; i > 0; i--)
297
294
      }
298
295
    offset += (ind[0] - bufferedRegionIndex[0]);
299
296
 
300
 
    return offset;    
301
 
  }
 
297
    return offset;
 
298
    }
302
299
#endif
303
300
  /** Compute the index of the pixel at a specified offset from the
304
301
   * beginning of the buffered region. Bounds checking is not
307
304
   * should be between 0 and the number of pixels in the
308
305
   * BufferedRegion (the latter can be found using
309
306
   * ImageRegion::GetNumberOfPixels()). */
310
 
#if 1
 
307
#ifdef ITK_USE_TEMPLATE_META_PROGRAMMING_LOOP_UNROLLING
311
308
  inline IndexType ComputeIndex(OffsetValueType offset) const
312
 
  {
 
309
    {
313
310
    IndexType index;
314
311
    const IndexType &bufferedRegionIndex = this->GetBufferedRegion().GetIndex();
315
312
    ImageHelper<VImageDimension,VImageDimension>::ComputeIndex(bufferedRegionIndex,
317
314
                                                               m_OffsetTable,
318
315
                                                               index);
319
316
    return index;
320
 
  }
 
317
    }
321
318
#else
322
319
  IndexType ComputeIndex(OffsetValueType offset) const
323
 
  {
 
320
    {
324
321
    IndexType index;
325
322
    const IndexType &bufferedRegionIndex = this->GetBufferedRegion().GetIndex();
326
 
    
 
323
 
327
324
    for (int i=VImageDimension-1; i > 0; i--)
328
325
      {
329
326
      index[i] = static_cast<IndexValueType>(offset / m_OffsetTable[i]);
332
329
      }
333
330
    index[0] = bufferedRegionIndex[0] + static_cast<IndexValueType>(offset);
334
331
 
335
 
    return index;    
336
 
  }
337
 
#endif
 
332
    return index;
 
333
    }
 
334
#endif
 
335
 
 
336
  /** Set the spacing (size of a pixel) of the image. The
 
337
   * spacing is the geometric distance between image samples.
 
338
   * It is stored internally as double, but may be set from
 
339
   * float. These methods also pre-compute the Index to Physical
 
340
   * point transforms of the image.
 
341
   * \sa GetSpacing() */
 
342
  virtual void SetSpacing (const SpacingType & spacing);
 
343
  virtual void SetSpacing (const double spacing[VImageDimension]);
 
344
  virtual void SetSpacing (const float spacing[VImageDimension]);
 
345
 
 
346
 
 
347
  /** Get the index (discrete) from a physical point.
 
348
   * Floating point index results are truncated to integers.
 
349
   * Returns true if the resulting index is within the image, false otherwise
 
350
   * \sa Transform */
 
351
#ifdef ITK_USE_TEMPLATE_META_PROGRAMMING_LOOP_UNROLLING
 
352
  template<class TCoordRep>
 
353
  bool TransformPhysicalPointToIndex(
 
354
    const Point<TCoordRep, VImageDimension>& point,
 
355
    IndexType & index ) const
 
356
    {
 
357
      ImageTransformHelper<VImageDimension,VImageDimension-1,VImageDimension-1>::TransformPhysicalPointToIndex(
 
358
        this->m_PhysicalPointToIndex, this->m_Origin, point, index);
 
359
 
 
360
    // Now, check to see if the index is within allowed bounds
 
361
    const bool isInside = this->GetLargestPossibleRegion().IsInside( index );
 
362
    return isInside;
 
363
    }
 
364
#else
 
365
  template<class TCoordRep>
 
366
  bool TransformPhysicalPointToIndex(
 
367
            const Point<TCoordRep, VImageDimension>& point,
 
368
            IndexType & index                                ) const
 
369
    {
 
370
    for (unsigned int i = 0; i < VImageDimension; i++)
 
371
      {
 
372
      TCoordRep sum = NumericTraits<TCoordRep>::Zero;
 
373
      for (unsigned int j = 0; j < VImageDimension; j++)
 
374
        {
 
375
        sum += this->m_PhysicalPointToIndex[i][j] * (point[j] - this->m_Origin[j]);
 
376
        }
 
377
      index[i] = static_cast< IndexValueType>( sum );
 
378
      }
 
379
 
 
380
    // Now, check to see if the index is within allowed bounds
 
381
    const bool isInside = this->GetLargestPossibleRegion().IsInside( index );
 
382
 
 
383
    return isInside;
 
384
    }
 
385
#endif
 
386
 
 
387
  /** \brief Get the continuous index from a physical point
 
388
   *
 
389
   * Returns true if the resulting index is within the image, false otherwise.
 
390
   * \sa Transform */
 
391
  template<class TCoordRep>
 
392
  bool TransformPhysicalPointToContinuousIndex(
 
393
              const Point<TCoordRep, VImageDimension>& point,
 
394
              ContinuousIndex<TCoordRep, VImageDimension>& index   ) const
 
395
    {
 
396
    Vector<double, VImageDimension> cvector;
 
397
 
 
398
    for( unsigned int k = 0; k < VImageDimension; k++ )
 
399
      {
 
400
      cvector[k] = point[k] - this->m_Origin[k];
 
401
      }
 
402
    cvector = m_PhysicalPointToIndex * cvector;
 
403
    for( unsigned int i = 0; i < VImageDimension; i++ )
 
404
      {
 
405
      index[i] = static_cast<TCoordRep>(cvector[i]);
 
406
      }
 
407
 
 
408
    // Now, check to see if the index is within allowed bounds
 
409
    const bool isInside = this->GetLargestPossibleRegion().IsInside( index );
 
410
 
 
411
    return isInside;
 
412
    }
 
413
 
 
414
 
 
415
  /** Get a physical point (in the space which
 
416
   * the origin and spacing infomation comes from)
 
417
   * from a continuous index (in the index space)
 
418
   * \sa Transform */
 
419
  template<class TCoordRep>
 
420
  void TransformContinuousIndexToPhysicalPoint(
 
421
            const ContinuousIndex<TCoordRep, VImageDimension>& index,
 
422
            Point<TCoordRep, VImageDimension>& point        ) const
 
423
    {
 
424
    for( unsigned int r=0; r<VImageDimension; r++)
 
425
      {
 
426
      TCoordRep sum = NumericTraits<TCoordRep>::Zero;
 
427
      for( unsigned int c=0; c<VImageDimension; c++ )
 
428
        {
 
429
        sum += this->m_IndexToPhysicalPoint(r,c) * index[c];
 
430
        }
 
431
      point[r] = sum + this->m_Origin[r];
 
432
      }
 
433
    }
 
434
 
 
435
  /** Get a physical point (in the space which
 
436
   * the origin and spacing infomation comes from)
 
437
   * from a discrete index (in the index space)
 
438
   *
 
439
   * \sa Transform */
 
440
#ifdef ITK_USE_TEMPLATE_META_PROGRAMMING_LOOP_UNROLLING
 
441
  template<class TCoordRep>
 
442
  void TransformIndexToPhysicalPoint(
 
443
                      const IndexType & index,
 
444
                      Point<TCoordRep, VImageDimension>& point ) const
 
445
    {
 
446
      ImageTransformHelper<VImageDimension,VImageDimension-1,VImageDimension-1>::TransformIndexToPhysicalPoint(
 
447
        this->m_IndexToPhysicalPoint, this->m_Origin, index, point);
 
448
    }
 
449
#else
 
450
  template<class TCoordRep>
 
451
  void TransformIndexToPhysicalPoint(
 
452
                      const IndexType & index,
 
453
                      Point<TCoordRep, VImageDimension>& point ) const
 
454
    {
 
455
    for (unsigned int i = 0; i < VImageDimension; i++)
 
456
      {
 
457
      point[i] = this->m_Origin[i];
 
458
      for (unsigned int j = 0; j < VImageDimension; j++)
 
459
        {
 
460
        point[i] += m_IndexToPhysicalPoint[i][j] * index[j];
 
461
        }
 
462
      }
 
463
    }
 
464
#endif
 
465
 
 
466
 
 
467
  /** Get a physical point (in the space which
 
468
   * the origin and spacing infomation comes from)
 
469
   * from a discrete index (in the index space)
 
470
   *
 
471
   * \sa Transform */
 
472
 
 
473
  /** Take a vector or covariant vector that has been computed in the
 
474
   * coordinate system parallel to the image grid and rotate it by the
 
475
   * direction cosines in order to get it in terms of the coordinate system of
 
476
   * the image acquisition device.  This implementation in the OrientedImage
 
477
   * multiply the array (vector or covariant vector) by the matrix of Direction
 
478
   * Cosines. The arguments of the method are of type FixedArray to make
 
479
   * possible to use this method with both Vector and CovariantVector.
 
480
   * The Method is implemented differently in the itk::Image.
 
481
   *
 
482
   * \sa Image
 
483
   */ 
 
484
  template<class TCoordRep>
 
485
  void TransformLocalVectorToPhysicalVector(
 
486
    const FixedArray<TCoordRep, VImageDimension> & inputGradient,
 
487
          FixedArray<TCoordRep, VImageDimension> & outputGradient ) const
 
488
    {
 
489
    //
 
490
    // This temporary implementation should be replaced with Template MetaProgramming.
 
491
    // 
 
492
#ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION
 
493
    const DirectionType & direction = this->GetDirection();
 
494
    for (unsigned int i = 0; i < VImageDimension; i++)
 
495
      {
 
496
      typedef typename NumericTraits<TCoordRep>::AccumulateType CoordSumType;
 
497
      CoordSumType sum = NumericTraits<CoordSumType>::Zero;
 
498
      for (unsigned int j = 0; j < VImageDimension; j++)
 
499
        {
 
500
        sum += direction[i][j] * inputGradient[j];
 
501
        }
 
502
      outputGradient[i] = static_cast<TCoordRep>( sum );
 
503
      }
 
504
#else
 
505
    for (unsigned int i = 0; i < VImageDimension; i++)
 
506
      {
 
507
      outputGradient[i] = inputGradient[i];
 
508
      }
 
509
#endif
 
510
    }
338
511
 
339
512
  /** Copy information from the specified data set.  This method is
340
513
   * part of the pipeline execution model. By default, a ProcessObject
393
566
   * throws a InvalidRequestedRegionError exception is the requested
394
567
   * region is not within the LargestPossibleRegion. */
395
568
  virtual bool VerifyRequestedRegion();
396
 
  
 
569
 
397
570
  /** INTERNAL This method is used internally by filters to copy meta-data from
398
571
   * the output to the input. Users should not have a need to use this method.
399
572
   *
406
579
   * \endcode
407
580
   *
408
581
   * \sa ImageBase, VectorImage
409
 
   * 
410
 
   * Returns/Sets the number of components in the image. Note that for all 
 
582
   *
 
583
   * Returns/Sets the number of components in the image. Note that for all
411
584
   * images this is 1. Even for Image< RGBPixel< T >, 3 >.
412
585
   * This is >= 1 only for time-series images such as itk::VectorImage. */
413
586
  virtual unsigned int GetNumberOfComponentsPerPixel() const;
414
 
  virtual void SetNumberOfComponentsPerPixel( unsigned int ); 
 
587
  virtual void SetNumberOfComponentsPerPixel( unsigned int );
415
588
 
416
589
protected:
417
590
  ImageBase();
424
597
   * the BufferedRegion is set. */
425
598
  void ComputeOffsetTable();
426
599
 
 
600
  /** Compute helper matrices used to transform Index coordinates to
 
601
   * PhysicalPoint coordinates and back. This method is virtual and will be
 
602
   * overloaded in derived classes in order to provide backward compatibility
 
603
   * behavior in classes that did not used to take image orientation into
 
604
   * account.  */ 
 
605
  virtual void ComputeIndexToPhysicalPointMatrices();
 
606
 
427
607
protected:
428
608
  /** Origin and spacing of physical coordinates. This variables are
429
609
   * protected for efficiency.  They are referenced frequently by
430
610
   * inner loop calculations. */
431
 
  SpacingType  m_Spacing;
432
 
  PointType   m_Origin;
433
 
  DirectionType m_Direction;
 
611
  SpacingType         m_Spacing;
 
612
  PointType           m_Origin;
 
613
  DirectionType       m_Direction;
 
614
 
 
615
  /** Matrices intended to help with the conversion of Index coordinates
 
616
   *  to PhysicalPoint coordinates */
 
617
  DirectionType       m_IndexToPhysicalPoint;
 
618
  DirectionType       m_PhysicalPointToIndex;
434
619
 
435
620
private:
436
621
  ImageBase(const Self&); //purposely not implemented
441
626
  RegionType          m_LargestPossibleRegion;
442
627
  RegionType          m_RequestedRegion;
443
628
  RegionType          m_BufferedRegion;
 
629
 
444
630
};
445
631
 
446
632
} // end namespace itk