1
/*=========================================================================
3
Program: Insight Segmentation & Registration Toolkit
4
Module: $RCSfile: itkProjectionImageFilter.txx,v $
6
Date: $Date: 2008-02-08 04:52:21 $
7
Version: $Revision: 1.1 $
9
Copyright (c) Insight Software Consortium. All rights reserved.
10
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
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.
16
=========================================================================*/
17
#ifndef __itkProjectionImageFilter_txx
18
#define __itkProjectionImageFilter_txx
20
#include "itkProjectionImageFilter.h"
21
#include "itkImageRegionIterator.h"
22
#include "itkImageRegionConstIterator.h"
23
#include "itkImageRegionConstIteratorWithIndex.h"
24
#include "itkProgressReporter.h"
25
#include "itkImageLinearConstIteratorWithIndex.h"
34
template <class TInputImage, class TOutputImage, class TAccumulator>
35
ProjectionImageFilter<TInputImage,TOutputImage,TAccumulator>
36
::ProjectionImageFilter()
38
this->SetNumberOfRequiredInputs( 1 );
39
m_ProjectionDimension=InputImageDimension-1;
43
template <class TInputImage, class TOutputImage, class TAccumulator>
45
ProjectionImageFilter<TInputImage,TOutputImage,TAccumulator>
46
::GenerateOutputInformation()
48
itkDebugMacro("GenerateOutputInformation Start");
50
if(m_ProjectionDimension>=TInputImage::ImageDimension)
52
itkExceptionMacro(<<"Invalid ProjectionDimension. ProjectionDimension is "
53
<< m_ProjectionDimension
54
<< " but input ImageDimension is "
55
<< TInputImage::ImageDimension);
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;
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() );
73
inputIndex = input->GetLargestPossibleRegion().GetIndex();
74
inputSize = input->GetLargestPossibleRegion().GetSize();
75
inSpacing = input->GetSpacing();
76
inOrigin = input->GetOrigin();
78
// Set the LargestPossibleRegion of the output.
79
// Reduce the size of the accumulated dimension.
81
if( static_cast< unsigned int >( InputImageDimension ) ==
82
static_cast< unsigned int >( OutputImageDimension ) )
84
for(unsigned int i = 0; i<InputImageDimension; i++)
86
if (i != m_ProjectionDimension)
88
outputSize[i] = inputSize[i];
89
outputIndex[i] = inputIndex[i];
90
outSpacing[i] = inSpacing[i];
91
outOrigin[i] = inOrigin[i];
97
outSpacing[i] = inSpacing[i]*inputSize[i];
98
outOrigin[i] = inOrigin[i] + (i-1)*inSpacing[i]/2;
104
// Then OutputImageDimension = InputImageDimension - 1
105
for(unsigned int i = 0; i<OutputImageDimension; i++)
107
if (i != m_ProjectionDimension)
109
outputSize[i] = inputSize[i];
110
outputIndex[i] = inputIndex[i];
111
outSpacing[i] = inSpacing[i];
112
outOrigin[i] = inOrigin[i];
116
outputSize[i] = inputSize[InputImageDimension - 1];
117
outputIndex[i] = inputIndex[InputImageDimension - 1];
118
outSpacing[i] = inSpacing[InputImageDimension - 1];
119
outOrigin[i] = inOrigin[InputImageDimension - 1];
124
outputRegion.SetSize(outputSize);
125
outputRegion.SetIndex(outputIndex);
126
output->SetOrigin(outOrigin);
127
output->SetSpacing(outSpacing);
128
output->SetLargestPossibleRegion(outputRegion);
130
itkDebugMacro("GenerateOutputInformation End");
134
template <class TInputImage, class TOutputImage, class TAccumulator>
136
ProjectionImageFilter<TInputImage,TOutputImage,TAccumulator>
137
::GenerateInputRequestedRegion()
139
itkDebugMacro("GenerateInputRequestedRegion Start");
141
if(m_ProjectionDimension>=TInputImage::ImageDimension)
143
itkExceptionMacro(<<"Invalid ProjectionDimension "
144
<< m_ProjectionDimension
145
<< " but ImageDimension is "
146
<< TInputImage::ImageDimension);
149
Superclass::GenerateInputRequestedRegion();
151
if ( this->GetInput() )
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;
161
outputIndex = this->GetOutput()->GetRequestedRegion().GetIndex();
162
outputSize = this->GetOutput()->GetRequestedRegion().GetSize();
163
inputLargSize = this->GetInput()->GetLargestPossibleRegion().GetSize();
164
inputLargIndex = this->GetInput()->GetLargestPossibleRegion().GetIndex();
166
if( static_cast< unsigned int >( InputImageDimension ) ==
167
static_cast< unsigned int >( OutputImageDimension ) )
169
for(unsigned int i=0; i<TInputImage::ImageDimension; i++)
171
if(i!=m_ProjectionDimension)
173
inputSize[i] = outputSize[i];
174
inputIndex[i] = outputIndex[i];
178
inputSize[i]=inputLargSize[i];
179
inputIndex[i]=inputLargIndex[i];
185
for(unsigned int i=0; i<OutputImageDimension; i++)
187
if(i!=m_ProjectionDimension)
189
inputSize[i] = outputSize[i];
190
inputIndex[i] = outputIndex[i];
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];
200
inputSize[m_ProjectionDimension] =
201
inputLargSize[m_ProjectionDimension];
202
inputIndex[m_ProjectionDimension] =
203
inputLargIndex[m_ProjectionDimension];
206
RequestedRegion.SetSize(inputSize);
207
RequestedRegion.SetIndex(inputIndex);
208
InputImagePointer input = const_cast< TInputImage * > ( this->GetInput() );
209
input->SetRequestedRegion (RequestedRegion);
212
itkDebugMacro("GenerateInputRequestedRegion End");
217
* GenerateData Performs the accumulation
219
template <class TInputImage, class TOutputImage, class TAccumulator>
221
ProjectionImageFilter<TInputImage,TOutputImage,TAccumulator>
222
::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
225
if(m_ProjectionDimension>=TInputImage::ImageDimension)
227
itkExceptionMacro(<<"Invalid ProjectionDimension "
228
<< m_ProjectionDimension
229
<< " but ImageDimension is "
230
<< TInputImage::ImageDimension);
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());
238
typedef typename TOutputImage::PixelType OutputPixelType;
240
// get some values, just to be easier to manipulate
241
typename Superclass::InputImageConstPointer inputImage = this->GetInput();
243
typename TInputImage::RegionType inputRegion =
244
inputImage->GetLargestPossibleRegion();
246
typename TInputImage::SizeType inputSize = inputRegion.GetSize();
247
typename TInputImage::IndexType inputIndex = inputRegion.GetIndex();
249
typename TOutputImage::Pointer outputImage = this->GetOutput();
250
typename TOutputImage::RegionType outputRegion =
251
outputImage->GetLargestPossibleRegion();
253
typename TOutputImage::SizeType outputSizeForThread =
254
outputRegionForThread.GetSize();
255
typename TOutputImage::IndexType outputIndexForThread =
256
outputRegionForThread.GetIndex();
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;
263
if( static_cast< unsigned int >( InputImageDimension ) ==
264
static_cast< unsigned int >( OutputImageDimension ) )
266
for( unsigned int i=0; i< InputImageDimension; i++ )
268
if( i != m_ProjectionDimension )
270
inputSizeForThread[i] = outputSizeForThread[i];
271
inputIndexForThread[i] = outputIndexForThread[i];
277
for(unsigned int i=0; i<OutputImageDimension; i++)
279
if(i!=m_ProjectionDimension)
281
inputSizeForThread[i] = outputSizeForThread[i];
282
inputIndexForThread[i] = outputIndexForThread[i];
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];
292
inputSizeForThread[m_ProjectionDimension] =
293
inputSize[m_ProjectionDimension];
294
inputIndexForThread[m_ProjectionDimension] =
295
inputIndex[m_ProjectionDimension];
297
inputRegionForThread.SetSize( inputSizeForThread );
298
inputRegionForThread.SetIndex( inputIndexForThread );
300
unsigned long projectionSize = inputSize[m_ProjectionDimension];
302
// create the iterators for input and output image
303
typedef ImageLinearConstIteratorWithIndex<TInputImage> InputIteratorType;
304
InputIteratorType iIt( inputImage, inputRegionForThread );
305
iIt.SetDirection( m_ProjectionDimension );
308
// instantiate the accumulator
309
AccumulatorType accumulator = this->NewAccumulator( projectionSize );
311
// ok, everything is ready... lets the linear iterator do its job !
312
while( !iIt.IsAtEnd() )
314
// init the accumulator before a new set of pixels
315
accumulator.Initialize();
317
while( !iIt.IsAtEndOfLine() )
319
accumulator( iIt.Get() );
323
// move the ouput iterator and set the output value
324
typename TOutputImage::IndexType oIdx;
325
typename TInputImage::IndexType iIdx = iIt.GetIndex();
327
if( static_cast< unsigned int >( InputImageDimension ) ==
328
static_cast< unsigned int >( OutputImageDimension ) )
330
for( unsigned int i=0; i< InputImageDimension; i++ )
332
if( i != m_ProjectionDimension )
344
for( unsigned int i=0; i< OutputImageDimension; i++ )
346
if( i != m_ProjectionDimension )
352
oIdx[i] = iIdx[InputImageDimension - 1];
357
outputImage->SetPixel( oIdx,
358
static_cast<OutputPixelType>( accumulator.GetValue() ) );
360
// one more line done !
361
progress.CompletedPixel();
363
// continue with the next one
370
template <class TInputImage, class TOutputImage, class TAccumulator>
372
ProjectionImageFilter<TInputImage,TOutputImage,TAccumulator>::
373
NewAccumulator( unsigned long size ) const
375
return TAccumulator( size );
379
template <class TInputImage, class TOutputImage, class TAccumulator>
381
ProjectionImageFilter<TInputImage,TOutputImage,TAccumulator>::
382
PrintSelf(std::ostream& os, Indent indent) const
384
Superclass::PrintSelf(os,indent);
386
os << indent << "ProjectionDimension: " << m_ProjectionDimension << std::endl;
390
} // end namespace itk