1 /*===================================================================
3 The Medical Imaging Interaction Toolkit (MITK)
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
13 See LICENSE.txt or http://www.mitk.org for details.
15 ===================================================================*/
17 #ifndef _itkRegularizedIVIMReconstructionSingleIteration_txx
18 #define _itkRegularizedIVIMReconstructionSingleIteration_txx
21 #include "itkConstShapedNeighborhoodIterator.h"
22 #include "itkNeighborhoodInnerProduct.h"
23 #include "itkImageRegionIterator.h"
24 #include "itkNeighborhoodAlgorithm.h"
25 #include "itkZeroFluxNeumannBoundaryCondition.h"
26 #include "itkOffset.h"
27 #include "itkProgressReporter.h"
28 #include "itkRegularizedIVIMLocalVariationImageFilter.h"
34 #define IVIM_FOO -100000
42 template <class TInputPixel, class TOutputPixel, class TRefPixelType>
43 RegularizedIVIMReconstructionSingleIteration<TInputPixel, TOutputPixel, TRefPixelType>
44 ::RegularizedIVIMReconstructionSingleIteration()
47 m_LocalVariation = LocalVariationImageType::New();
51 * generate requested region
53 template <class TInputPixel, class TOutputPixel, class TRefPixelType>
55 RegularizedIVIMReconstructionSingleIteration<TInputPixel, TOutputPixel, TRefPixelType>
56 ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError)
58 // call the superclass' implementation of this method
59 Superclass::GenerateInputRequestedRegion();
61 // get pointers to the input and output
62 typename Superclass::InputImagePointer inputPtr =
63 const_cast< InputImageType * >( this->GetInput() );
64 typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
66 if ( !inputPtr || !outputPtr )
71 // get a copy of the input requested region (should equal the output
73 typename InputImageType::RegionType inputRequestedRegion;
74 inputRequestedRegion = inputPtr->GetRequestedRegion();
76 // pad the input requested region by 1
77 inputRequestedRegion.PadByRadius( 1 );
79 // crop the input requested region at the input's largest possible region
80 if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
82 inputPtr->SetRequestedRegion( inputRequestedRegion );
87 // Couldn't crop the region (requested region is outside the largest
88 // possible region). Throw an exception.
90 // store what we tried to request (prior to trying to crop)
91 inputPtr->SetRequestedRegion( inputRequestedRegion );
94 InvalidRequestedRegionError e(__FILE__, __LINE__);
95 e.SetLocation(ITK_LOCATION);
96 e.SetDescription("Requested region outside possible region.");
97 e.SetDataObject(inputPtr);
106 template <class TInputPixel, class TOutputPixel, class TRefPixelType>
108 RegularizedIVIMReconstructionSingleIteration<TInputPixel, TOutputPixel, TRefPixelType>
109 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
113 typename OutputImageType::Pointer output = this->GetOutput();
114 typename InputImageType::ConstPointer input = this->GetInput();
116 // Find the data-set boundary "faces"
117 itk::Size<InputImageDimension> size;
118 for( int i=0; i<InputImageDimension; i++)
121 NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> bC;
122 typename NeighborhoodAlgorithm::
123 ImageBoundaryFacesCalculator<InputImageType>::FaceListType
124 faceList = bC(input, outputRegionForThread, size);
126 NeighborhoodAlgorithm::
127 ImageBoundaryFacesCalculator<LocalVariationImageType> lv_bC;
128 typename NeighborhoodAlgorithm::
129 ImageBoundaryFacesCalculator<LocalVariationImageType>::FaceListType
130 lv_faceList = lv_bC(m_LocalVariation, outputRegionForThread, size);
132 ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;
133 ZeroFluxNeumannBoundaryCondition<LocalVariationImageType> lv_nbc;
134 std::vector<double> ws;
135 std::vector<double> hs;
137 typename NeighborhoodAlgorithm::
138 ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator
139 lv_fit=lv_faceList.begin();
141 // Process each of the boundary faces. These are N-d regions which border
142 // the edge of the buffer.
143 for ( typename NeighborhoodAlgorithm::
144 ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator
145 fit=faceList.begin(); fit != faceList.end(); ++fit)
148 // iterators over output, input, original and local variation image
149 ImageRegionIterator<OutputImageType> output_image_it =
150 ImageRegionIterator<OutputImageType>(output, *fit);
151 ImageRegionConstIterator<InputImageType> input_image_it =
152 ImageRegionConstIterator<InputImageType>(input, *fit);
153 ImageRegionConstIterator<RefImageType> orig_image_it =
154 ImageRegionConstIterator<RefImageType>(m_OriginalImage, *fit);
155 ImageRegionConstIterator<LocalVariationImageType> loc_var_image_it =
156 ImageRegionConstIterator<LocalVariationImageType>(
157 m_LocalVariation, *fit);
159 // neighborhood in input image
160 ConstShapedNeighborhoodIterator<InputImageType>
161 input_image_neighbors_it(size, input, *fit);
162 typename ConstShapedNeighborhoodIterator<InputImageType>::
164 input_image_neighbors_it.OverrideBoundaryCondition(&nbc);
165 input_image_neighbors_it.ClearActiveList();
166 for(int i=0; i<InputImageDimension; i++)
170 input_image_neighbors_it.ActivateOffset(offset);
172 input_image_neighbors_it.ActivateOffset(offset);
174 input_image_neighbors_it.GoToBegin();
176 // neighborhood in local variation image
177 ConstShapedNeighborhoodIterator<LocalVariationImageType>
178 loc_var_image_neighbors_it(size, m_LocalVariation, *lv_fit);
179 loc_var_image_neighbors_it.OverrideBoundaryCondition(&lv_nbc);
180 loc_var_image_neighbors_it.ClearActiveList();
181 for(int i=0; i<InputImageDimension; i++)
185 loc_var_image_neighbors_it.ActivateOffset(offset);
187 loc_var_image_neighbors_it.ActivateOffset(offset);
189 loc_var_image_neighbors_it.GoToBegin();
191 const unsigned int neighborhoodSize = InputImageDimension*2;
192 ws.resize(neighborhoodSize);
194 while ( ! output_image_it.IsAtEnd() )
197 //MITK_INFO << "looking at voxel " << output_image_it.GetIndex();
199 // 1 / ||nabla_alpha(u)||_a
200 double locvar_alpha_inv = 1.0/loc_var_image_it.Get();
201 //MITK_INFO << "locvar: " << loc_var_image_it.Get();
202 // compute w_alphabetas
205 typename ConstShapedNeighborhoodIterator<LocalVariationImageType>::
206 ConstIterator loc_var_neighbors_it;
207 for (loc_var_neighbors_it = loc_var_image_neighbors_it.Begin();
208 ! loc_var_neighbors_it.IsAtEnd();
209 loc_var_neighbors_it++)
212 // 1 / ||nabla_alpha(u)||_a + 1 / ||nabla_beta(u)||_a
214 locvar_alpha_inv + (1.0/(double)loc_var_neighbors_it.Get());
217 //MITK_INFO << "nb var: " << count << ": " << loc_var_neighbors_it.Get();
219 //MITK_INFO << "wsum: " << wsum;
221 // h_alphaalpha * u_alpha^zero
222 typename RefImageType::PixelType orig = orig_image_it.Get();
223 // vnl_vector<double> estim(orig.GetSize());
224 // vnl_matrix<double> diff(orig.GetSize(),1);
225 // vnl_matrix<double> estimdash(orig.GetSize(),2);
226 vnl_vector_fixed<double,2> step;
227 step[0] = 0; step[1]=0;
228 for(size_t ind=0; ind<m_BValues.size(); ind++)
230 //MITK_INFO << "refval: " << orig[ind];
231 double estim = (1-input_image_it.Get()[0])*exp(-m_BValues[ind]*input_image_it.Get()[1]);
232 double estimdash1 = exp(-m_BValues[ind]*input_image_it.Get()[1]);
233 double estimdash2 = (-1.0) * (1.0-input_image_it.Get()[0]) * m_BValues[ind] * exp(-m_BValues[ind]*input_image_it.Get()[1]);
234 //MITK_INFO << "estimdash1: " << estimdash1 << "; estimdash2: " << estimdash2 << "; diff: " << (double) orig[ind] - estim;
235 if(orig[ind] != IVIM_FOO)
237 step[0] += ((double) orig[ind] - estim) * estimdash2;
238 step[1] += ((double) orig[ind] - estim) * estimdash1;
242 step[1] *= m_Lambda / (m_Lambda+wsum);
244 // add the different h_alphabeta * u_beta
246 typename ConstShapedNeighborhoodIterator<InputImageType>::
247 ConstIterator input_neighbors_it;
248 for (input_neighbors_it = input_image_neighbors_it.Begin();
249 ! input_neighbors_it.IsAtEnd();
250 input_neighbors_it++)
252 step[1] += (input_neighbors_it.Get()[1] - input_image_it.Get()[1]) * (ws[count++] / (m_Lambda+wsum));
255 //MITK_INFO << "stepfinal: " << step[0] << "; " << step[1];
259 out[0] = input_image_it.Get()[0] + .001*step[0];
260 out[1] = input_image_it.Get()[1] + .00001*step[1];
261 output_image_it.Set( out );
263 //MITK_INFO << "(" << input_image_it.Get()[0] << " ; " << input_image_it.Get()[1] << ") => (" << out[0] << " ; " << out[1] << ")";
264 // increment iterators
269 ++input_image_neighbors_it;
270 ++loc_var_image_neighbors_it;
279 * first calculate local variation in the image
281 template <class TInputPixel, class TOutputPixel, class TRefPixelType>
283 RegularizedIVIMReconstructionSingleIteration<TInputPixel, TOutputPixel, TRefPixelType>
284 ::BeforeThreadedGenerateData()
286 typedef typename itk::RegularizedIVIMLocalVariationImageFilter
287 <InputImageType,LocalVariationImageType> FilterType;
288 typename FilterType::Pointer filter = FilterType::New();
289 filter->SetInput(this->GetInput(0));
290 filter->SetNumberOfThreads(this->GetNumberOfThreads());
292 this->m_LocalVariation = filter->GetOutput();
296 * Standard "PrintSelf" method
298 template <class TInputPixel, class TOutputPixel, class TRefPixelType>
300 RegularizedIVIMReconstructionSingleIteration<TInputPixel, TOutputPixel, TRefPixelType>
305 Superclass::PrintSelf( os, indent );
308 } // end namespace itk