1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://www.centreleonberard.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================*/
18 #ifndef __clitkExplosionControlledThresholdConnectedImageFilter_txx
19 #define __clitkExplosionControlledThresholdConnectedImageFilter_txx
21 #include "clitkExplosionControlledThresholdConnectedImageFilter.h"
22 #include "itkBinaryThresholdImageFunction.h"
23 #include "itkFloodFilledImageFunctionConditionalIterator.h"
24 #include "itkShapedFloodFilledImageFunctionConditionalIterator.h"
25 #include "itkProgressReporter.h"
30 //--------------------------------------------------------------------
34 template <class TInputImage, class TOutputImage>
35 ExplosionControlledThresholdConnectedImageFilter<TInputImage, TOutputImage>
36 ::ExplosionControlledThresholdConnectedImageFilter()
38 m_Lower = itk::NumericTraits<InputImagePixelType>::NonpositiveMin();
39 m_Upper = itk::NumericTraits<InputImagePixelType>::max();
40 m_FinalLower = itk::NumericTraits<InputImagePixelType>::NonpositiveMin();
41 m_FinalUpper = itk::NumericTraits<InputImagePixelType>::max();
42 m_AdaptLowerBorder=false;
43 m_AdaptUpperBorder=false;
44 m_MaximumUpperThreshold=itk::NumericTraits<InputImagePixelType>::max();
45 m_MinimumLowerThreshold=itk::NumericTraits<InputImagePixelType>::NonpositiveMin();
47 m_ReplaceValue = itk::NumericTraits<OutputImagePixelType>::One;
48 m_MinimumThresholdStepSize=itk::NumericTraits<InputImagePixelType>::One;
49 m_ThresholdStepSize=64;
51 m_FullyConnected=false;
54 //--------------------------------------------------------------------
57 //--------------------------------------------------------------------
58 template <class TInputImage, class TOutputImage>
60 ExplosionControlledThresholdConnectedImageFilter<TInputImage, TOutputImage>
63 if( this->m_Seeds.size() > 0 )
65 this->m_Seeds.clear();
69 //--------------------------------------------------------------------
72 //--------------------------------------------------------------------
73 template <class TInputImage, class TOutputImage>
75 ExplosionControlledThresholdConnectedImageFilter<TInputImage, TOutputImage>
76 ::SetSeed(const IndexType & seed)
79 this->AddSeed ( seed );
81 //--------------------------------------------------------------------
84 //--------------------------------------------------------------------
85 template <class TInputImage, class TOutputImage>
87 ExplosionControlledThresholdConnectedImageFilter<TInputImage, TOutputImage>
88 ::AddSeed ( const IndexType & seed )
90 this->m_Seeds.push_back ( seed );
93 //--------------------------------------------------------------------
96 //--------------------------------------------------------------------
98 * Standard PrintSelf method.
100 template <class TInputImage, class TOutputImage>
102 ExplosionControlledThresholdConnectedImageFilter<TInputImage, TOutputImage>
103 ::PrintSelf(std::ostream& os, itk::Indent indent) const
105 this->Superclass::PrintSelf(os, indent);
106 os << indent << "Upper: "
107 << static_cast<typename itk::NumericTraits<InputImagePixelType>::PrintType>(m_Upper)
109 os << indent << "Lower: "
110 << static_cast<typename itk::NumericTraits<InputImagePixelType>::PrintType>(m_Lower)
112 os << indent << "ReplaceValue: "
113 << static_cast<double>(m_Multiplier)
115 os << indent << "ReplaceValue: "
116 << static_cast<typename itk::NumericTraits<OutputImagePixelType>::PrintType>(m_ReplaceValue)
119 //--------------------------------------------------------------------
122 //--------------------------------------------------------------------
123 template <class TInputImage, class TOutputImage>
125 ExplosionControlledThresholdConnectedImageFilter<TInputImage,TOutputImage>
126 ::GenerateInputRequestedRegion()
128 Superclass::GenerateInputRequestedRegion();
129 if ( this->GetInput() )
131 InputImagePointer image =
132 const_cast< InputImageType * >( this->GetInput() );
133 image->SetRequestedRegionToLargestPossibleRegion();
136 //--------------------------------------------------------------------
139 //--------------------------------------------------------------------
140 template <class TInputImage, class TOutputImage>
142 ExplosionControlledThresholdConnectedImageFilter<TInputImage,TOutputImage>
143 ::EnlargeOutputRequestedRegion(itk::DataObject *output)
145 Superclass::EnlargeOutputRequestedRegion(output);
146 output->SetRequestedRegionToLargestPossibleRegion();
148 //--------------------------------------------------------------------
151 //--------------------------------------------------------------------
152 template <class TInputImage, class TOutputImage>
154 ExplosionControlledThresholdConnectedImageFilter<TInputImage,TOutputImage>
157 typename Superclass::InputImageConstPointer inputImage = this->GetInput();
158 typename Superclass::OutputImagePointer outputImage = this->GetOutput();
161 outputImage->SetBufferedRegion( outputImage->GetRequestedRegion() );
162 outputImage->Allocate();
163 outputImage->FillBuffer ( itk::NumericTraits<OutputImagePixelType>::Zero );
166 //--------------------------------------------------
167 // Get initial region size
168 //--------------------------------------------------
169 typedef itk::BinaryThresholdImageFunction<InputImageType> FunctionType;
170 //typedef itk::ShapedFloodFilledImageFunctionConditionalIterator<OutputImageType, FunctionType> IteratorType;
171 typedef itk::FloodFilledImageFunctionConditionalIterator<OutputImageType, FunctionType> IteratorType;
172 typename FunctionType::Pointer function = FunctionType::New();
173 function->SetInputImage ( inputImage );
174 function->ThresholdBetween ( m_Lower, m_Upper );
175 IteratorType it = IteratorType ( outputImage, function, m_Seeds );
176 //it.SetFullyConnected(m_FullyConnected);
177 unsigned int initialSize=0;
180 while( !it.IsAtEnd())
185 if (m_Verbose) std::cout<<"Initial region size using thresholds ["<<m_Lower<<", "<<m_Upper<<"], is "<<initialSize<<"..."<<std::endl;
188 //--------------------------------------------------
189 // Decrease lower threshold
190 //--------------------------------------------------
191 m_FinalLower=m_Lower;
192 if (m_AdaptLowerBorder)
195 InputImagePixelType currentLower=m_Lower;
196 unsigned int currentSize=initialSize;
197 unsigned int previousSize=initialSize;
199 // Lower the threshold till explosion
200 while ( (previousSize<m_MinimumSize) || ( (currentLower> m_MinimumLowerThreshold) && ( (unsigned int) currentSize <= (unsigned int )m_Multiplier* previousSize) ) )
203 currentLower-=m_ThresholdStepSize;
204 function->ThresholdBetween(currentLower, m_Upper);
206 // Get current region size
207 previousSize=currentSize;
209 IteratorType it = IteratorType ( outputImage, function, m_Seeds );
211 while( !it.IsAtEnd())
216 if (m_Verbose) std::cout<<"Decreasing lower threshold to "<<currentLower<<", size is "<<currentSize <<"(previously "<<previousSize<<") ..."<<std::endl;
220 // Explosion occured, increase lower theshold
221 if ( (double)currentSize > m_Multiplier*(double) previousSize)
223 // Raise lower threshold
224 if (m_Verbose) std::cout<<"Explosion detected, adapting lower threshold ..."<<std::endl;
225 currentLower+=m_ThresholdStepSize;
226 InputImagePixelType currentStepSize=m_ThresholdStepSize;
228 while (currentStepSize>m_MinimumThresholdStepSize)
231 currentLower-=currentStepSize;
233 // Get current region size
235 IteratorType it = IteratorType ( outputImage, function, m_Seeds );
237 while( !it.IsAtEnd())
243 if (m_Verbose) std::cout<<"Adapting lower threshold to "<<currentLower<<", size is "<<currentSize <<"(previously "<<previousSize<<") ..."<<std::endl;
245 // Explosion: go back
246 if ((double)currentSize > m_Multiplier* (double) previousSize)
247 currentLower+=currentStepSize;
249 // No explosion: adapt previous
250 else previousSize=currentSize;
254 m_FinalLower=currentLower;
255 if (m_Verbose) std::cout<<"Final lower threshold (precision="<<m_MinimumThresholdStepSize<<") is set to "<<m_FinalLower<<"..."<<std::endl;
258 if (m_Verbose) std::cout<<"No explosion before minimum lower threshold reached!"<<std::endl;
261 } // end update lower
264 //--------------------------------------------------
265 // Increase upper threshold
266 //--------------------------------------------------
267 m_FinalUpper=m_Upper;
268 if (m_AdaptUpperBorder)
271 InputImagePixelType currentUpper=m_Upper;
272 unsigned int currentSize=initialSize;
273 unsigned int previousSize=initialSize;
275 // Upper the threshold till explosion
276 while ((previousSize<m_MinimumSize) || ( (currentUpper< m_MaximumUpperThreshold) && ((double)currentSize <= m_Multiplier* (double) previousSize) ) )
279 currentUpper+=m_ThresholdStepSize;
280 function->ThresholdBetween(m_Lower,currentUpper);
282 // Get current region size
283 previousSize=currentSize;
285 IteratorType it = IteratorType ( outputImage, function, m_Seeds );
287 while( !it.IsAtEnd())
292 if (m_Verbose) std::cout<<"Increasing upper threshold to "<<currentUpper<<", size is "<<currentSize <<"(previously "<<previousSize<<") ..."<<std::endl;
296 // Explosion occured, decrease upper theshold
297 if ((double)currentSize > m_Multiplier* (double) previousSize)
299 // Lower upper threshold
300 if (m_Verbose) std::cout<<"Explosion detected, adapting upper threshold ..."<<std::endl;
301 currentUpper-=m_ThresholdStepSize;
302 InputImagePixelType currentStepSize=m_ThresholdStepSize;
304 while (currentStepSize>m_MinimumThresholdStepSize)
307 currentUpper+=currentStepSize;
309 // Get current region size
311 function->ThresholdBetween(m_Lower,currentUpper);
312 IteratorType it = IteratorType ( outputImage, function, m_Seeds );
314 while( !it.IsAtEnd())
319 if (m_Verbose) std::cout<<"Adapting upper threshold to "<<currentUpper<<", size is "<<currentSize <<"(previously "<<previousSize<<") ..."<<std::endl;
321 // Explosion: go back
322 if ((double)currentSize > m_Multiplier* (double) previousSize)
323 currentUpper-=currentStepSize;
325 // No explosion: adapt previous
326 else previousSize=currentSize;
330 m_FinalUpper=currentUpper;
331 if (m_Verbose) std::cout<<"Final upper threshold (precision="<<m_MinimumThresholdStepSize<<") is set to "<<m_FinalUpper<<"..."<<std::endl;
334 else { if (m_Verbose) std::cout<<"No explosion before maximum upper threshold reached!"<<std::endl; }
336 } // end update upper
340 //--------------------------------------------------
341 // Update the output with final thresholds
342 //--------------------------------------------------
343 function->ThresholdBetween(m_FinalLower, m_FinalUpper);
344 IteratorType it2 = IteratorType ( outputImage, function, m_Seeds );
346 while( !it2.IsAtEnd())
348 it2.Set(m_ReplaceValue);
353 //--------------------------------------------------------------------
356 } // end namespace clitk