itkTotalVariationSingleIterationImageFilter.txx
Go to the documentation of this file.
00001 /*========================================================================= 00002 00003 Program: Insight Segmentation & Registration Toolkit 00004 Language: C++ 00005 Date: $Date: 2006-01-11 19:43:31 $ 00006 Version: $Revision: x $ 00007 00008 Copyright (c) Insight Software Consortium. All rights reserved. 00009 See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. 00010 00011 This software is distributed WITHOUT ANY WARRANTY; without even 00012 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 00013 PURPOSE. See the above copyright notices for more information. 00014 00015 =========================================================================*/ 00016 00017 #ifndef _itkTotalVariationSingleIterationImageFilter_txx 00018 #define _itkTotalVariationSingleIterationImageFilter_txx 00019 00020 #include "itkTotalVariationSingleIterationImageFilter.h" 00021 00022 // itk includes 00023 #include "itkConstShapedNeighborhoodIterator.h" 00024 #include "itkNeighborhoodInnerProduct.h" 00025 #include "itkImageRegionIterator.h" 00026 #include "itkNeighborhoodAlgorithm.h" 00027 #include "itkZeroFluxNeumannBoundaryCondition.h" 00028 #include "itkOffset.h" 00029 #include "itkProgressReporter.h" 00030 #include "itkLocalVariationImageFilter.h" 00031 00032 // other includes 00033 #include <vector> 00034 #include <algorithm> 00035 00036 namespace itk 00037 { 00038 00042 template <class TInputImage, class TOutputImage> 00043 TotalVariationSingleIterationImageFilter<TInputImage, TOutputImage> 00044 ::TotalVariationSingleIterationImageFilter() 00045 { 00046 m_Lambda = 1.0; 00047 m_LocalVariation = LocalVariationImageType::New(); 00048 } 00049 00053 template <class TInputImage, class TOutputImage> 00054 void 00055 TotalVariationSingleIterationImageFilter<TInputImage, TOutputImage> 00056 ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError) 00057 { 00058 // call the superclass' implementation of this method 00059 Superclass::GenerateInputRequestedRegion(); 00060 00061 // get pointers to the input and output 00062 typename Superclass::InputImagePointer inputPtr = 00063 const_cast< TInputImage * >( this->GetInput() ); 00064 typename Superclass::OutputImagePointer outputPtr = this->GetOutput(); 00065 00066 if ( !inputPtr || !outputPtr ) 00067 { 00068 return; 00069 } 00070 00071 // get a copy of the input requested region (should equal the output 00072 // requested region) 00073 typename TInputImage::RegionType inputRequestedRegion; 00074 inputRequestedRegion = inputPtr->GetRequestedRegion(); 00075 00076 // pad the input requested region by 1 00077 inputRequestedRegion.PadByRadius( 1 ); 00078 00079 // crop the input requested region at the input's largest possible region 00080 if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) ) 00081 { 00082 inputPtr->SetRequestedRegion( inputRequestedRegion ); 00083 return; 00084 } 00085 else 00086 { 00087 // Couldn't crop the region (requested region is outside the largest 00088 // possible region). Throw an exception. 00089 00090 // store what we tried to request (prior to trying to crop) 00091 inputPtr->SetRequestedRegion( inputRequestedRegion ); 00092 00093 // build an exception 00094 InvalidRequestedRegionError e(__FILE__, __LINE__); 00095 e.SetLocation(ITK_LOCATION); 00096 e.SetDescription("Requested region outside possible region."); 00097 e.SetDataObject(inputPtr); 00098 throw e; 00099 } 00100 } 00101 00102 00106 template< class TInputImage, class TOutputImage> 00107 void 00108 TotalVariationSingleIterationImageFilter< TInputImage, TOutputImage> 00109 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, 00110 int threadId) 00111 { 00112 00113 typename OutputImageType::Pointer output = this->GetOutput(); 00114 typename InputImageType::ConstPointer input = this->GetInput(); 00115 00116 // Find the data-set boundary "faces" 00117 itk::Size<InputImageDimension> size; 00118 for( int i=0; i<InputImageDimension; i++) 00119 size[i] = 1; 00120 00121 NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> bC; 00122 typename NeighborhoodAlgorithm:: 00123 ImageBoundaryFacesCalculator<InputImageType>::FaceListType 00124 faceList = bC(input, outputRegionForThread, size); 00125 00126 NeighborhoodAlgorithm:: 00127 ImageBoundaryFacesCalculator<LocalVariationImageType> lv_bC; 00128 typename NeighborhoodAlgorithm:: 00129 ImageBoundaryFacesCalculator<LocalVariationImageType>::FaceListType 00130 lv_faceList = lv_bC(m_LocalVariation, outputRegionForThread, size); 00131 00132 // support progress methods/callbacks 00133 ProgressReporter progress( 00134 this, threadId, outputRegionForThread.GetNumberOfPixels()); 00135 00136 ZeroFluxNeumannBoundaryCondition<InputImageType> nbc; 00137 ZeroFluxNeumannBoundaryCondition<LocalVariationImageType> lv_nbc; 00138 std::vector<double> ws; 00139 std::vector<double> hs; 00140 00141 typename NeighborhoodAlgorithm:: 00142 ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator 00143 lv_fit=lv_faceList.begin(); 00144 00145 // Process each of the boundary faces. These are N-d regions which border 00146 // the edge of the buffer. 00147 for ( typename NeighborhoodAlgorithm:: 00148 ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator 00149 fit=faceList.begin(); fit != faceList.end(); ++fit) 00150 { 00151 00152 // iterators over output, input, original and local variation image 00153 ImageRegionIterator<OutputImageType> output_image_it = 00154 ImageRegionIterator<OutputImageType>(output, *fit); 00155 ImageRegionConstIterator<InputImageType> input_image_it = 00156 ImageRegionConstIterator<InputImageType>(input, *fit); 00157 ImageRegionConstIterator<InputImageType> orig_image_it = 00158 ImageRegionConstIterator<InputImageType>(m_OriginalImage, *fit); 00159 ImageRegionConstIterator<LocalVariationImageType> loc_var_image_it = 00160 ImageRegionConstIterator<LocalVariationImageType>( 00161 m_LocalVariation, *fit); 00162 00163 // neighborhood in input image 00164 ConstShapedNeighborhoodIterator<InputImageType> 00165 input_image_neighbors_it(size, input, *fit); 00166 typename ConstShapedNeighborhoodIterator<InputImageType>:: 00167 OffsetType offset; 00168 input_image_neighbors_it.OverrideBoundaryCondition(&nbc); 00169 input_image_neighbors_it.ClearActiveList(); 00170 for(int i=0; i<InputImageDimension; i++) 00171 { 00172 offset.Fill(0); 00173 offset[i] = -1; 00174 input_image_neighbors_it.ActivateOffset(offset); 00175 offset[i] = 1; 00176 input_image_neighbors_it.ActivateOffset(offset); 00177 } 00178 input_image_neighbors_it.GoToBegin(); 00179 00180 // neighborhood in local variation image 00181 ConstShapedNeighborhoodIterator<LocalVariationImageType> 00182 loc_var_image_neighbors_it(size, m_LocalVariation, *lv_fit); 00183 loc_var_image_neighbors_it.OverrideBoundaryCondition(&lv_nbc); 00184 loc_var_image_neighbors_it.ClearActiveList(); 00185 for(int i=0; i<InputImageDimension; i++) 00186 { 00187 offset.Fill(0); 00188 offset[i] = -1; 00189 loc_var_image_neighbors_it.ActivateOffset(offset); 00190 offset[i] = 1; 00191 loc_var_image_neighbors_it.ActivateOffset(offset); 00192 } 00193 loc_var_image_neighbors_it.GoToBegin(); 00194 00195 const unsigned int neighborhoodSize = InputImageDimension*2; 00196 ws.resize(neighborhoodSize); 00197 00198 while ( ! output_image_it.IsAtEnd() ) 00199 { 00200 00201 // 1 / ||nabla_alpha(u)||_a 00202 double locvar_alpha_inv = 1.0/loc_var_image_it.Get(); 00203 00204 // compute w_alphabetas 00205 int count = 0; 00206 double wsum = 0; 00207 typename ConstShapedNeighborhoodIterator<LocalVariationImageType>:: 00208 ConstIterator loc_var_neighbors_it; 00209 for (loc_var_neighbors_it = loc_var_image_neighbors_it.Begin(); 00210 ! loc_var_neighbors_it.IsAtEnd(); 00211 loc_var_neighbors_it++) 00212 { 00213 // w_alphabeta(u) = 00214 // 1 / ||nabla_alpha(u)||_a + 1 / ||nabla_beta(u)||_a 00215 ws[count] = 00216 locvar_alpha_inv + (1.0/(double)loc_var_neighbors_it.Get()); 00217 wsum += ws[count++]; 00218 } 00219 00220 // h_alphaalpha * u_alpha^zero 00221 typename OutputImageType::PixelType res = 00222 static_cast<typename OutputImageType::PixelType>( 00223 ((typename OutputImageType::PixelType) 00224 orig_image_it.Get()) * (m_Lambda / (m_Lambda+wsum))); 00225 00226 // add the different h_alphabeta * u_beta 00227 count = 0; 00228 typename ConstShapedNeighborhoodIterator<InputImageType>:: 00229 ConstIterator input_neighbors_it; 00230 for (input_neighbors_it = input_image_neighbors_it.Begin(); 00231 ! input_neighbors_it.IsAtEnd(); 00232 input_neighbors_it++) 00233 { 00234 res += input_neighbors_it.Get() * (ws[count++] / (m_Lambda+wsum)); 00235 } 00236 00237 // set output result 00238 output_image_it.Set( res ); 00239 00240 // increment iterators 00241 ++output_image_it; 00242 ++input_image_it; 00243 ++orig_image_it; 00244 ++loc_var_image_it; 00245 ++input_image_neighbors_it; 00246 ++loc_var_image_neighbors_it; 00247 00248 // report progress 00249 progress.CompletedPixel(); 00250 00251 } 00252 00253 ++lv_fit; 00254 } 00255 } 00256 00260 template< class TInputImage, class TOutputImage> 00261 void 00262 TotalVariationSingleIterationImageFilter< TInputImage, TOutputImage> 00263 ::BeforeThreadedGenerateData() 00264 { 00265 typedef typename itk::LocalVariationImageFilter 00266 <TInputImage,LocalVariationImageType> FilterType; 00267 typename FilterType::Pointer filter = FilterType::New(); 00268 filter->SetInput(this->GetInput(0)); 00269 filter->SetNumberOfThreads(this->GetNumberOfThreads()); 00270 filter->Update(); 00271 this->m_LocalVariation = filter->GetOutput(); 00272 } 00273 00277 template <class TInputImage, class TOutput> 00278 void 00279 TotalVariationSingleIterationImageFilter<TInputImage, TOutput> 00280 ::PrintSelf( 00281 std::ostream& os, 00282 Indent indent) const 00283 { 00284 Superclass::PrintSelf( os, indent ); 00285 } 00286 00287 } // end namespace itk 00288 00289 #endif
