1 /*=========================================================================
2  *
3  *  Copyright Insight Software Consortium
4  *
5  *  Licensed under the Apache License, Version 2.0 (the "License");
6  *  you may not use this file except in compliance with the License.
7  *  You may obtain a copy of the License at
8  *
9  *         http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  *  Unless required by applicable law or agreed to in writing, software
12  *  distributed under the License is distributed on an "AS IS" BASIS,
13  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  *  See the License for the specific language governing permissions and
15  *  limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef itkExtractImageFilter_hxx
19 #define itkExtractImageFilter_hxx
20 
21 #include "itkExtractImageFilter.h"
22 #include "itkImageAlgorithm.h"
23 #include "itkObjectFactory.h"
24 #include "itkProgressReporter.h"
25 
26 namespace itk
27 {
28 /**
29  *
30  */
31 template< typename TInputImage, typename TOutputImage >
32 ExtractImageFilter< TInputImage, TOutputImage >
ExtractImageFilter()33 ::ExtractImageFilter():
34   m_DirectionCollapseStrategy(DIRECTIONCOLLAPSETOUNKOWN)
35 {
36   Superclass::InPlaceOff();
37   this->DynamicMultiThreadingOn();
38 }
39 
40 /**
41  *
42  */
43 template< typename TInputImage, typename TOutputImage >
44 void
45 ExtractImageFilter< TInputImage, TOutputImage >
PrintSelf(std::ostream & os,Indent indent) const46 ::PrintSelf(std::ostream & os, Indent indent) const
47 {
48   Superclass::PrintSelf(os, indent);
49 
50   os << indent << "ExtractionRegion: " << m_ExtractionRegion << std::endl;
51   os << indent << "OutputImageRegion: " << m_OutputImageRegion << std::endl;
52   os << indent << "DirectionCollapseStrategy: " << m_DirectionCollapseStrategy << std::endl;
53 }
54 
55 template< typename TInputImage, typename TOutputImage >
56 void
57 ExtractImageFilter< TInputImage, TOutputImage >
CallCopyOutputRegionToInputRegion(InputImageRegionType & destRegion,const OutputImageRegionType & srcRegion)58 ::CallCopyOutputRegionToInputRegion(InputImageRegionType & destRegion,
59                                     const OutputImageRegionType & srcRegion)
60 {
61   ExtractImageFilterRegionCopierType extractImageRegionCopier;
62 
63   extractImageRegionCopier(destRegion, srcRegion, m_ExtractionRegion);
64 }
65 
66 template< typename TInputImage, typename TOutputImage >
67 void
68 ExtractImageFilter< TInputImage, TOutputImage >
SetExtractionRegion(InputImageRegionType extractRegion)69 ::SetExtractionRegion(InputImageRegionType extractRegion)
70 {
71   static_assert(InputImageDimension >= OutputImageDimension, "InputImageDimension must be greater than OutputImageDimension");
72   m_ExtractionRegion = extractRegion;
73 
74   unsigned int         nonzeroSizeCount = 0;
75   InputImageSizeType   inputSize = extractRegion.GetSize();
76   OutputImageSizeType  outputSize;
77   outputSize.Fill(0);
78   OutputImageIndexType outputIndex;
79   outputIndex.Fill(0);
80 
81   /**
82    * check to see if the number of non-zero entries in the extraction region
83    * matches the number of dimensions in the output image.
84    */
85   for ( unsigned int i = 0; i < InputImageDimension; ++i )
86     {
87     if ( inputSize[i] )
88       {
89       outputSize[nonzeroSizeCount] = inputSize[i];
90       outputIndex[nonzeroSizeCount] = extractRegion.GetIndex()[i];
91       nonzeroSizeCount++;
92       }
93     }
94 
95   if ( nonzeroSizeCount != OutputImageDimension )
96     {
97     itkExceptionMacro("The number of zero sized dimensions in the input image Extraction Region\n"
98                       << "is not consistent with the dimensionality of the output image.\n"
99                       << "Expected the extraction region size (" << extractRegion.GetSize()
100                       << ") to contain " << InputImageDimension-OutputImageDimension
101                       << " zero sized dimensions to collapse." );
102     }
103 
104   m_OutputImageRegion.SetSize(outputSize);
105   m_OutputImageRegion.SetIndex(outputIndex);
106   this->Modified();
107 }
108 
109 /**
110  * ExtractImageFilter can produce an image which is a different resolution
111  * than its input image.  As such, ExtractImageFilter needs to provide an
112  * implementation for GenerateOutputInformation() in order to inform
113  * the pipeline execution model.  The original documentation of this
114  * method is below.
115  *
116  * \sa ProcessObject::GenerateOutputInformaton()
117  */
118 template< typename TInputImage, typename TOutputImage >
119 void
120 ExtractImageFilter< TInputImage, TOutputImage >
GenerateOutputInformation()121 ::GenerateOutputInformation()
122 {
123   // do not call the superclass' implementation of this method since
124   // this filter allows the input and the output to be of different dimensions
125 
126   // get pointers to the input and output
127   typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
128   typename Superclass::InputImageConstPointer inputPtr  = this->GetInput();
129 
130   if ( !outputPtr || !inputPtr )
131     {
132     return;
133     }
134 
135   // Set the output image size to the same value as the extraction region.
136   outputPtr->SetLargestPossibleRegion(m_OutputImageRegion);
137 
138   // Set the output spacing and origin
139   const ImageBase< InputImageDimension > *phyData;
140 
141   phyData =
142     dynamic_cast< const ImageBase< InputImageDimension > * >( this->GetInput() );
143 
144   if ( phyData )
145     {
146     // Copy what we can from the image from spacing and origin of the input
147     // This logic needs to be augmented with logic that select which
148     // dimensions to copy
149 
150     const typename InputImageType::SpacingType &
151     inputSpacing = inputPtr->GetSpacing();
152     const typename InputImageType::DirectionType &
153     inputDirection = inputPtr->GetDirection();
154     const typename InputImageType::PointType &
155     inputOrigin = inputPtr->GetOrigin();
156 
157     typename OutputImageType::SpacingType outputSpacing;
158     typename OutputImageType::DirectionType outputDirection;
159     typename OutputImageType::PointType outputOrigin;
160     outputOrigin.Fill(0.0);
161 
162     if ( static_cast< unsigned int >( OutputImageDimension ) >
163          static_cast< unsigned int >( InputImageDimension ) )
164       {
165       // copy the input to the output and fill the rest of the
166       // output with zeros.
167       for ( unsigned int i = 0; i < InputImageDimension; ++i )
168         {
169         outputSpacing[i] = inputSpacing[i];
170         outputOrigin[i] = inputOrigin[i];
171         for ( unsigned int dim = 0; dim < InputImageDimension; ++dim )
172           {
173           outputDirection[i][dim] = inputDirection[i][dim];
174           }
175         }
176       for (unsigned int i=InputImageDimension; i < OutputImageDimension; ++i )
177         {
178         outputSpacing[i] = 1.0;
179         outputOrigin[i] = 0.0;
180         for ( unsigned int dim = 0; dim < InputImageDimension; ++dim )
181           {
182           outputDirection[i][dim] = 0.0;
183           }
184         outputDirection[i][i] = 1.0;
185         }
186       }
187     else
188       {
189       // copy the non-collapsed part of the input spacing and origing to the
190       // output
191       outputDirection.SetIdentity();
192       int nonZeroCount = 0;
193       for ( unsigned int i = 0; i < InputImageDimension; ++i )
194         {
195         if ( m_ExtractionRegion.GetSize()[i] )
196           {
197           outputSpacing[nonZeroCount] = inputSpacing[i];
198           outputOrigin[nonZeroCount] = inputOrigin[i];
199           int nonZeroCount2 = 0;
200           for ( unsigned int dim = 0; dim < InputImageDimension; ++dim )
201             {
202             if ( m_ExtractionRegion.GetSize()[dim] )
203               {
204               outputDirection[nonZeroCount][nonZeroCount2] =
205                 inputDirection[i][dim];
206               ++nonZeroCount2;
207               }
208             }
209           nonZeroCount++;
210           }
211         }
212       }
213     // if the filter changes from a higher to a lower dimension, or
214     // if, after rebuilding the direction cosines, there's a zero
215     // length cosine vector, reset the directions to identity
216     // or throw an exception, depending on the collapse strategy.
217     if( static_cast<int>(InputImageDimension) != static_cast<int>(OutputImageDimension) )
218       {
219       switch(m_DirectionCollapseStrategy)
220         {
221       case DIRECTIONCOLLAPSETOIDENTITY:
222           {
223           outputDirection.SetIdentity();
224           }
225         break;
226       case DIRECTIONCOLLAPSETOSUBMATRIX:
227           {
228           if ( vnl_determinant( outputDirection.GetVnlMatrix() ) == 0.0 )
229             {
230             itkExceptionMacro( << "Invalid submatrix extracted for collapsed direction." );
231             }
232           }
233         break;
234       case DIRECTIONCOLLAPSETOGUESS:
235           {
236           if ( vnl_determinant( outputDirection.GetVnlMatrix() ) == 0.0 )
237             {
238             outputDirection.SetIdentity();
239             }
240           }
241         break;
242       case DIRECTIONCOLLAPSETOUNKOWN:
243       default:
244           {
245           itkExceptionMacro( << "It is required that the strategy for collapsing the direction matrix be explicitly specified. "
246             << "Set with either myfilter->SetDirectionCollapseToIdentity() or myfilter->SetDirectionCollapseToSubmatrix() "
247             << typeid( ImageBase< InputImageDimension > * ).name() );
248           }
249         }
250       }
251     // set the spacing and origin
252     outputPtr->SetSpacing(outputSpacing);
253     outputPtr->SetDirection(outputDirection);
254     outputPtr->SetOrigin(outputOrigin);
255     outputPtr->SetNumberOfComponentsPerPixel(
256       inputPtr->GetNumberOfComponentsPerPixel() );
257     }
258   else
259     {
260     // pointer could not be cast back down
261     itkExceptionMacro( << "itk::ExtractImageFilter::GenerateOutputInformation "
262                        << "cannot cast input to "
263                        << typeid( ImageBase< InputImageDimension > * ).name() );
264     }
265 }
266 
267 template< typename TInputImage, typename TOutputImage >
268 void
269 ExtractImageFilter< TInputImage, TOutputImage >
GenerateData()270 ::GenerateData()
271 {
272 
273   // InPlace::AllocateOutputs set the running in place ivar.
274   // This method will be called again, by GenerateData, but there is
275   // no harm done.
276   this->AllocateOutputs();
277 
278   // The input matched the output, nothing to do.
279   if ( this->GetRunningInPlace() )
280     {
281     OutputImageType *outputPtr = this->GetOutput();
282 
283     // the in-place grafting copies the meta data, this needs to be
284     // set back.
285     outputPtr->SetLargestPossibleRegion(m_OutputImageRegion);
286 
287     this->UpdateProgress( 1.0 );
288     return;
289     }
290 
291   this->Superclass::GenerateData();
292 }
293 
294 /**
295  * ExtractImageFilter can be implemented as a multithreaded filter.
296  * Therefore, this implementation provides a DynamicThreadedGenerateData()
297  * routine which is called for each processing thread. The output
298  * image data is allocated automatically by the superclass prior to
299  * calling DynamicThreadedGenerateData().  DynamicThreadedGenerateData can only
300  * write to the portion of the output image specified by the
301  * parameter "outputRegionForThread"
302  *
303  * \sa ImageToImageFilter::ThreadedGenerateData(),
304  *     ImageToImageFilter::GenerateData()
305  */
306 template< typename TInputImage, typename TOutputImage >
307 void
308 ExtractImageFilter< TInputImage, TOutputImage >
DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread)309 ::DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread)
310 {
311   itkDebugMacro(<< "Actually executing");
312 
313   const InputImageType *inputPtr = this->GetInput();
314   OutputImageType      *outputPtr = this->GetOutput();
315 
316 
317   // Define the portion of the input to walk for this thread
318   InputImageRegionType inputRegionForThread;
319   this->CallCopyOutputRegionToInputRegion(inputRegionForThread, outputRegionForThread);
320 
321   // copy the input pixel to the output
322   ImageAlgorithm::Copy( inputPtr, outputPtr, inputRegionForThread, outputRegionForThread );
323 }
324 } // end namespace itk
325 
326 #endif
327