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

« back to all changes in this revision

Viewing changes to Code/BasicFilters/itkProjectionImageFilter.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
* 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:
 
1
/*=========================================================================
 
2
 
 
3
  Program:   Insight Segmentation & Registration Toolkit
 
4
  Module:    $RCSfile: itkProjectionImageFilter.txx,v $
 
5
  Language:  C++
 
6
  Date:      $Date: 2008-02-08 04:52:21 $
 
7
  Version:   $Revision: 1.1 $
 
8
 
 
9
  Copyright (c) Insight Software Consortium. All rights reserved.
 
10
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
 
11
 
 
12
     This software is distributed WITHOUT ANY WARRANTY; without even
 
13
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 
14
     PURPOSE.  See the above copyright notices for more information.
 
15
 
 
16
=========================================================================*/
 
17
#ifndef __itkProjectionImageFilter_txx
 
18
#define __itkProjectionImageFilter_txx
 
19
 
 
20
#include "itkProjectionImageFilter.h"
 
21
#include "itkImageRegionIterator.h"
 
22
#include "itkImageRegionConstIterator.h"
 
23
#include "itkImageRegionConstIteratorWithIndex.h"
 
24
#include "itkProgressReporter.h"
 
25
#include "itkImageLinearConstIteratorWithIndex.h"
 
26
 
 
27
 
 
28
namespace itk
 
29
{
 
30
 
 
31
/**
 
32
 * Constructor
 
33
 */
 
34
template <class TInputImage, class TOutputImage, class TAccumulator>
 
35
ProjectionImageFilter<TInputImage,TOutputImage,TAccumulator>
 
36
::ProjectionImageFilter()
 
37
{
 
38
  this->SetNumberOfRequiredInputs( 1 );
 
39
  m_ProjectionDimension=InputImageDimension-1;
 
40
}
 
41
 
 
42
 
 
43
template <class TInputImage, class TOutputImage, class TAccumulator>
 
44
void
 
45
ProjectionImageFilter<TInputImage,TOutputImage,TAccumulator>
 
46
::GenerateOutputInformation()
 
47
{
 
48
  itkDebugMacro("GenerateOutputInformation Start");
 
49
 
 
50
  if(m_ProjectionDimension>=TInputImage::ImageDimension)
 
51
    {
 
52
    itkExceptionMacro(<<"Invalid ProjectionDimension. ProjectionDimension is "
 
53
                      << m_ProjectionDimension
 
54
                      << " but input ImageDimension is "
 
55
                      << TInputImage::ImageDimension);
 
56
    }
 
57
 
 
58
  typename TOutputImage::RegionType outputRegion;
 
59
  typename TInputImage::IndexType inputIndex;
 
60
  typename TInputImage::SizeType  inputSize;
 
61
  typename TOutputImage::SizeType  outputSize;
 
62
  typename TOutputImage::IndexType outputIndex;
 
63
  typename TInputImage::SpacingType inSpacing;
 
64
  typename TInputImage::PointType inOrigin;
 
65
  typename TOutputImage::SpacingType outSpacing;
 
66
  typename TOutputImage::PointType outOrigin;
 
67
 
 
68
  // Get pointers to the input and output
 
69
  typename Superclass::OutputImagePointer output = this->GetOutput();
 
70
  typename Superclass::InputImagePointer input = 
 
71
    const_cast< TInputImage * >( this->GetInput() );
 
72
 
 
73
  inputIndex = input->GetLargestPossibleRegion().GetIndex();
 
74
  inputSize = input->GetLargestPossibleRegion().GetSize();
 
75
  inSpacing = input->GetSpacing();
 
76
  inOrigin = input->GetOrigin();
 
77
 
 
78
  // Set the LargestPossibleRegion of the output.
 
79
  // Reduce the size of the accumulated dimension.
 
80
 
 
81
  if( static_cast< unsigned int >( InputImageDimension ) == 
 
82
      static_cast< unsigned int >( OutputImageDimension )    )
 
83
    {
 
84
    for(unsigned int i = 0; i<InputImageDimension; i++)
 
85
      {
 
86
      if (i != m_ProjectionDimension)
 
87
        {
 
88
        outputSize[i]  = inputSize[i];
 
89
        outputIndex[i] = inputIndex[i];
 
90
        outSpacing[i] = inSpacing[i];
 
91
        outOrigin[i]  = inOrigin[i];
 
92
        }
 
93
      else
 
94
        {
 
95
        outputSize[i]  = 1;
 
96
        outputIndex[i] = 0;
 
97
        outSpacing[i] = inSpacing[i]*inputSize[i];
 
98
        outOrigin[i]  = inOrigin[i] + (i-1)*inSpacing[i]/2;
 
99
        }
 
100
      }
 
101
    }
 
102
  else
 
103
    {
 
104
    // Then OutputImageDimension = InputImageDimension - 1
 
105
    for(unsigned int i = 0; i<OutputImageDimension; i++)
 
106
      {
 
107
      if (i != m_ProjectionDimension)
 
108
        {
 
109
        outputSize[i]  = inputSize[i];
 
110
        outputIndex[i] = inputIndex[i];
 
111
        outSpacing[i] = inSpacing[i];
 
112
        outOrigin[i]  = inOrigin[i];
 
113
        }
 
114
      else
 
115
        {
 
116
        outputSize[i]  = inputSize[InputImageDimension - 1];
 
117
        outputIndex[i] = inputIndex[InputImageDimension - 1];
 
118
        outSpacing[i] = inSpacing[InputImageDimension - 1];
 
119
        outOrigin[i]  = inOrigin[InputImageDimension - 1];
 
120
        }
 
121
      }
 
122
    }
 
123
 
 
124
  outputRegion.SetSize(outputSize);
 
125
  outputRegion.SetIndex(outputIndex);
 
126
  output->SetOrigin(outOrigin);
 
127
  output->SetSpacing(outSpacing);
 
128
  output->SetLargestPossibleRegion(outputRegion);
 
129
 
 
130
  itkDebugMacro("GenerateOutputInformation End");
 
131
}
 
132
 
 
133
 
 
134
template <class TInputImage, class  TOutputImage, class TAccumulator>
 
135
void
 
136
ProjectionImageFilter<TInputImage,TOutputImage,TAccumulator>
 
137
::GenerateInputRequestedRegion()
 
138
{
 
139
  itkDebugMacro("GenerateInputRequestedRegion Start");
 
140
 
 
141
  if(m_ProjectionDimension>=TInputImage::ImageDimension)
 
142
    {
 
143
    itkExceptionMacro(<<"Invalid ProjectionDimension "
 
144
                      << m_ProjectionDimension
 
145
                      << " but ImageDimension is "
 
146
                      << TInputImage::ImageDimension);
 
147
    }
 
148
 
 
149
  Superclass::GenerateInputRequestedRegion();
 
150
 
 
151
  if ( this->GetInput() )
 
152
    {
 
153
    typename TInputImage::RegionType RequestedRegion;
 
154
    typename TInputImage::SizeType  inputSize;
 
155
    typename TInputImage::IndexType inputIndex;
 
156
    typename TInputImage::SizeType  inputLargSize;
 
157
    typename TInputImage::IndexType inputLargIndex;
 
158
    typename TOutputImage::SizeType  outputSize;
 
159
    typename TOutputImage::IndexType outputIndex;
 
160
 
 
161
    outputIndex = this->GetOutput()->GetRequestedRegion().GetIndex();
 
162
    outputSize = this->GetOutput()->GetRequestedRegion().GetSize();
 
163
    inputLargSize = this->GetInput()->GetLargestPossibleRegion().GetSize();
 
164
    inputLargIndex = this->GetInput()->GetLargestPossibleRegion().GetIndex();
 
165
 
 
166
    if( static_cast< unsigned int >( InputImageDimension ) == 
 
167
        static_cast< unsigned int >( OutputImageDimension )    )
 
168
      {
 
169
      for(unsigned int i=0; i<TInputImage::ImageDimension; i++)
 
170
        {
 
171
        if(i!=m_ProjectionDimension)
 
172
          {
 
173
          inputSize[i] = outputSize[i];
 
174
          inputIndex[i] = outputIndex[i];
 
175
          }
 
176
        else
 
177
          {
 
178
          inputSize[i]=inputLargSize[i];
 
179
          inputIndex[i]=inputLargIndex[i];
 
180
          }
 
181
        }
 
182
      }
 
183
    else
 
184
      {
 
185
      for(unsigned int i=0; i<OutputImageDimension; i++)
 
186
        {
 
187
        if(i!=m_ProjectionDimension)
 
188
          {
 
189
          inputSize[i] = outputSize[i];
 
190
          inputIndex[i] = outputIndex[i];
 
191
          }
 
192
        else
 
193
          {
 
194
          // the size of the output image on this dimension is the size
 
195
          // of the input image on the removed dimension
 
196
          inputSize[InputImageDimension - 1] = outputSize[i];
 
197
          inputIndex[InputImageDimension - 1] = outputIndex[i];
 
198
          }
 
199
        }
 
200
        inputSize[m_ProjectionDimension] = 
 
201
                     inputLargSize[m_ProjectionDimension];
 
202
        inputIndex[m_ProjectionDimension] = 
 
203
                     inputLargIndex[m_ProjectionDimension];
 
204
      }
 
205
 
 
206
    RequestedRegion.SetSize(inputSize);
 
207
    RequestedRegion.SetIndex(inputIndex);
 
208
    InputImagePointer input = const_cast< TInputImage * > ( this->GetInput() );
 
209
    input->SetRequestedRegion (RequestedRegion);
 
210
    }
 
211
 
 
212
  itkDebugMacro("GenerateInputRequestedRegion End");
 
213
}
 
214
 
 
215
 
 
216
/**
 
217
 * GenerateData Performs the accumulation
 
218
 */
 
219
template <class TInputImage, class TOutputImage, class TAccumulator>
 
220
void
 
221
ProjectionImageFilter<TInputImage,TOutputImage,TAccumulator>
 
222
::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, 
 
223
                       int threadId )
 
224
{
 
225
  if(m_ProjectionDimension>=TInputImage::ImageDimension)
 
226
    {
 
227
    itkExceptionMacro(<<"Invalid ProjectionDimension "
 
228
                      << m_ProjectionDimension
 
229
                      << " but ImageDimension is "
 
230
                      << TInputImage::ImageDimension);
 
231
    }
 
232
 
 
233
  // use the output image to report the progress: there is no need to 
 
234
  // call CompletedPixel() for all input pixel
 
235
  ProgressReporter progress(this, threadId, 
 
236
                            outputRegionForThread.GetNumberOfPixels());
 
237
 
 
238
  typedef typename TOutputImage::PixelType OutputPixelType;
 
239
  
 
240
  // get some values, just to be easier to manipulate
 
241
  typename Superclass::InputImageConstPointer  inputImage = this->GetInput();
 
242
 
 
243
  typename TInputImage::RegionType inputRegion = 
 
244
                       inputImage->GetLargestPossibleRegion();
 
245
 
 
246
  typename TInputImage::SizeType inputSize = inputRegion.GetSize();
 
247
  typename TInputImage::IndexType inputIndex = inputRegion.GetIndex();
 
248
 
 
249
  typename TOutputImage::Pointer outputImage = this->GetOutput();
 
250
  typename TOutputImage::RegionType outputRegion = 
 
251
                       outputImage->GetLargestPossibleRegion();
 
252
 
 
253
  typename TOutputImage::SizeType outputSizeForThread = 
 
254
                       outputRegionForThread.GetSize();
 
255
  typename TOutputImage::IndexType outputIndexForThread = 
 
256
                       outputRegionForThread.GetIndex();
 
257
 
 
258
  // compute the input region for this thread
 
259
  typename TInputImage::RegionType inputRegionForThread = inputRegion;
 
260
  typename TInputImage::SizeType inputSizeForThread = inputSize;
 
261
  typename TInputImage::IndexType inputIndexForThread = inputIndex;
 
262
  
 
263
  if( static_cast< unsigned int >( InputImageDimension ) == 
 
264
      static_cast< unsigned int >( OutputImageDimension )    )
 
265
    {
 
266
    for( unsigned int i=0; i< InputImageDimension; i++ )
 
267
      {
 
268
      if( i != m_ProjectionDimension )
 
269
        {
 
270
        inputSizeForThread[i] = outputSizeForThread[i];
 
271
        inputIndexForThread[i] = outputIndexForThread[i];
 
272
        }
 
273
      }
 
274
    }
 
275
  else
 
276
    {
 
277
    for(unsigned int i=0; i<OutputImageDimension; i++)
 
278
      {
 
279
      if(i!=m_ProjectionDimension)
 
280
        {
 
281
        inputSizeForThread[i] = outputSizeForThread[i];
 
282
        inputIndexForThread[i] = outputIndexForThread[i];
 
283
        }
 
284
      else
 
285
        {
 
286
        // the size of the output image on this dimension is the size
 
287
        // of the input image on the removed dimension
 
288
        inputSizeForThread[InputImageDimension - 1] = outputSizeForThread[i];
 
289
        inputIndexForThread[InputImageDimension - 1] = outputIndexForThread[i];
 
290
        }
 
291
      }
 
292
      inputSizeForThread[m_ProjectionDimension] = 
 
293
                                 inputSize[m_ProjectionDimension];
 
294
      inputIndexForThread[m_ProjectionDimension] = 
 
295
                                 inputIndex[m_ProjectionDimension];
 
296
    }
 
297
  inputRegionForThread.SetSize( inputSizeForThread );
 
298
  inputRegionForThread.SetIndex( inputIndexForThread );
 
299
 
 
300
  unsigned long projectionSize = inputSize[m_ProjectionDimension];
 
301
 
 
302
  // create the iterators for input and output image
 
303
  typedef ImageLinearConstIteratorWithIndex<TInputImage> InputIteratorType;
 
304
  InputIteratorType iIt( inputImage, inputRegionForThread );
 
305
  iIt.SetDirection( m_ProjectionDimension );
 
306
  iIt.GoToBegin();
 
307
 
 
308
  // instantiate the accumulator
 
309
  AccumulatorType accumulator = this->NewAccumulator( projectionSize );
 
310
 
 
311
  // ok, everything is ready... lets the linear iterator do its job !
 
312
  while( !iIt.IsAtEnd() )
 
313
    {
 
314
    // init the accumulator before a new set of pixels
 
315
    accumulator.Initialize();
 
316
 
 
317
    while( !iIt.IsAtEndOfLine() )
 
318
      {
 
319
      accumulator( iIt.Get() );
 
320
      ++iIt;
 
321
      }
 
322
 
 
323
    // move the ouput iterator and set the output value
 
324
    typename TOutputImage::IndexType oIdx;
 
325
    typename TInputImage::IndexType iIdx = iIt.GetIndex();
 
326
 
 
327
    if( static_cast< unsigned int >( InputImageDimension ) == 
 
328
        static_cast< unsigned int >( OutputImageDimension )    )
 
329
      {
 
330
      for( unsigned int i=0; i< InputImageDimension; i++ )
 
331
        {
 
332
        if( i != m_ProjectionDimension )
 
333
          { 
 
334
          oIdx[i] = iIdx[i]; 
 
335
          }
 
336
        else
 
337
          { 
 
338
          oIdx[i] = 0; 
 
339
          }
 
340
        }
 
341
      }
 
342
    else
 
343
      {
 
344
      for( unsigned int i=0; i< OutputImageDimension; i++ )
 
345
        {
 
346
        if( i != m_ProjectionDimension )
 
347
          {
 
348
          oIdx[i] = iIdx[i]; 
 
349
          }
 
350
        else
 
351
          { 
 
352
          oIdx[i] = iIdx[InputImageDimension - 1]; 
 
353
          }
 
354
        }
 
355
      }
 
356
 
 
357
    outputImage->SetPixel( oIdx, 
 
358
                        static_cast<OutputPixelType>( accumulator.GetValue() ) );
 
359
 
 
360
    // one more line done !
 
361
    progress.CompletedPixel();
 
362
 
 
363
    // continue with the next one
 
364
    iIt.NextLine();
 
365
    }
 
366
 
 
367
}
 
368
 
 
369
 
 
370
template <class TInputImage, class TOutputImage, class TAccumulator>
 
371
TAccumulator
 
372
ProjectionImageFilter<TInputImage,TOutputImage,TAccumulator>::
 
373
NewAccumulator( unsigned long size ) const
 
374
{
 
375
  return TAccumulator( size );
 
376
}
 
377
 
 
378
 
 
379
template <class TInputImage, class TOutputImage, class TAccumulator>
 
380
void
 
381
ProjectionImageFilter<TInputImage,TOutputImage,TAccumulator>::
 
382
PrintSelf(std::ostream& os, Indent indent) const
 
383
{
 
384
  Superclass::PrintSelf(os,indent);
 
385
 
 
386
  os << indent << "ProjectionDimension: " << m_ProjectionDimension << std::endl;
 
387
}
 
388
 
 
389
 
 
390
} // end namespace itk
 
391
 
 
392
 
 
393
#endif