Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
mitkDWIHeadMotionCorrectionFilter.cpp
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 
17 
18 #ifndef MITKDWIHEADMOTIONCORRECTIONFILTER_CPP
19 #define MITKDWIHEADMOTIONCORRECTIONFILTER_CPP
20 
22 
23 #include "itkSplitDWImageFilter.h"
25 
26 #include "mitkImageTimeSelector.h"
27 
29 //#include "mitkRegistrationMethodITK4.h"
31 
33 #include <mitkImageWriteAccessor.h>
34 #include <mitkProperties.h>
35 
36 #include <vector>
37 
38 #include "mitkIOUtil.h"
39 #include <itkImage.h>
40 
42  : m_CurrentStep(0),
43  m_Steps(100),
44  m_IsInValidState(true),
45  m_AbortRegistration(false),
46  m_AverageUnweighted(true)
47 {
48 
49 }
50 
51 
53 {
55 
56  InputImageType* input = const_cast<InputImageType*>(this->GetInput(0));
57 
59  mitk::CastToItkImage(input, itkVectorImagePointer);
60 
61  typedef mitk::PyramidImageRegistrationMethod RegistrationMethod;
62  // typedef mitk::RegistrationMethodITKV4 RegistrationMethod
63  unsigned int numberOfSteps = itkVectorImagePointer->GetNumberOfComponentsPerPixel () ;
64  m_Steps = numberOfSteps;
65  //
66  // (1) Extract the b-zero images to a 3d+t image, register them by NCorr metric and
67  // rigid registration : they will then be used are reference image for registering
68  // the gradient images
69  //
72  b0_extractor->SetInput( itkVectorImagePointer );
73  b0_extractor->SetDirections( static_cast<mitk::GradientDirectionsProperty*>( input->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer() );
74  b0_extractor->Update();
75 
77  b0Image->InitializeByItk( b0_extractor->GetOutput() );
78  b0Image->SetImportChannel( b0_extractor->GetOutput()->GetBufferPointer(),
80 
81  // (2.1) Use the extractor to access the extracted b0 volumes
84 
85  t_selector->SetInput( b0Image );
86  t_selector->SetTimeNr(0);
87  t_selector->Update();
88 
89  // first unweighted image as reference space for the registration
90  mitk::Image::Pointer b0referenceImage = t_selector->GetOutput();
91  const unsigned int numberOfb0Images = b0Image->GetTimeSteps();
92 
93  // register b-zeros only if the flag to average is set, use the first one otherwise
94  if( m_AverageUnweighted && numberOfb0Images > 1)
95  {
96 
98  registrationMethod->SetFixedImage( b0referenceImage );
99  registrationMethod->SetTransformToRigid();
100 
101  // the unweighted images are of same modality
102  registrationMethod->SetCrossModalityOff();
103 
104  // use the advanced (windowed sinc) interpolation
105  registrationMethod->SetUseAdvancedInterpolation(true);
106 
107  // Initialize the temporary output image
108  mitk::Image::Pointer registeredB0Image = b0Image->Clone();
109 
112 
113  t_selector2->SetInput( b0Image );
114 
115  for( unsigned int i=1; i<numberOfb0Images; i++)
116  {
117 
118  m_CurrentStep = i + 1;
119  if( m_AbortRegistration == true) {
120  m_IsInValidState = false;
121  mitkThrow() << "Stopped by user.";
122  };
123 
124  t_selector2->SetTimeNr(i);
125  t_selector2->Update();
126 
127  registrationMethod->SetMovingImage( t_selector2->GetOutput() );
128 
129  try
130  {
131  MITK_INFO << " === (" << i <<"/"<< numberOfb0Images-1 << ") :: Starting registration";
132  registrationMethod->Update();
133  }
134  catch( const itk::ExceptionObject& e)
135  {
136  m_IsInValidState = false;
137  mitkThrow() << "Failed to register the b0 images, the PyramidRegistration threw an exception: \n" << e.what();
138  }
139 
140  // import volume to the inter-results
141  mitk::ImageWriteAccessor imac(registrationMethod->GetResampledMovingImage());
142  registeredB0Image->SetImportVolume( imac.GetData(),
144  }
145 
146  // use the accumulateImageFilter as provided by the ItkAccumulateFilter method in the header file
147  AccessFixedDimensionByItk_1(registeredB0Image, ItkAccumulateFilter, (4), b0referenceImage );
148 
149  }
150 
151 
152  //
153  // (2) Split the diffusion image into a 3d+t regular image, extract only the weighted images
154  //
156  split_filter->SetInput (itkVectorImagePointer );
157  split_filter->SetExtractAllAboveThreshold(8, static_cast<mitk::BValueMapProperty*>(input->GetProperty(mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME.c_str()).GetPointer() )->GetBValueMap() );
158 
159  try
160  {
161  split_filter->Update();
162  }
163  catch( const itk::ExceptionObject &e)
164  {
165  m_IsInValidState = false;
166  mitkThrow() << " Caught exception from SplitImageFilter : " << e.what();
167  }
168 
169  mitk::Image::Pointer splittedImage = mitk::Image::New();
170  splittedImage->InitializeByItk( split_filter->GetOutput() );
171  splittedImage->SetImportChannel( split_filter->GetOutput()->GetBufferPointer(),
173 
174 
175  //
176  // (3) Use again the time-selector to access the components separately in order
177  // to perform the registration of Image -> unweighted reference
178  //
179 
180  RegistrationMethod::Pointer weightedRegistrationMethod
182 
183  weightedRegistrationMethod->SetTransformToAffine();
184  weightedRegistrationMethod->SetCrossModalityOn();
185 
186  //
187  // - (3.1) Set the reference image
188  // - a single b0 image
189  // - average over the registered b0 images if multiple present
190  //
191  weightedRegistrationMethod->SetFixedImage( b0referenceImage );
192  // use the advanced (windowed sinc) interpolation
193  weightedRegistrationMethod->SetUseAdvancedInterpolation(true);
194  weightedRegistrationMethod->SetVerboseOn();
195 
196  //
197  // - (3.2) Register all timesteps in the splitted image onto the first reference
198  //
199  unsigned int maxImageIdx = splittedImage->GetTimeSteps();
200  mitk::TimeGeometry* tsg = splittedImage->GetTimeGeometry();
201  mitk::ProportionalTimeGeometry* ptg = dynamic_cast<ProportionalTimeGeometry*>(tsg);
202  ptg->Expand(maxImageIdx+1);
203  ptg->SetTimeStepGeometry( ptg->GetGeometryForTimeStep(0), maxImageIdx );
204 
205 
206  mitk::Image::Pointer registeredWeighted = mitk::Image::New();
207  registeredWeighted->Initialize( splittedImage->GetPixelType(0), *tsg );
208 
209  // insert the first unweighted reference as the first volume
210  // in own scope to release the accessor asap after copy
211  {
212  mitk::ImageWriteAccessor imac(b0referenceImage);
213  registeredWeighted->SetImportVolume( imac.GetData(),
215  }
216 
217 
218  // mitk::Image::Pointer registeredWeighted = splittedImage->Clone();
219  // this time start at 0, we have only gradient images in the 3d+t file
220  // the reference image comes form an other image
221  mitk::ImageTimeSelector::Pointer t_selector_w =
223 
224  t_selector_w->SetInput( splittedImage );
225 
226  // store the rotation parts of the transformations in a vector
227  typedef RegistrationMethod::TransformMatrixType MatrixType;
228  std::vector< MatrixType > estimated_transforms;
229 
230 
231  for( unsigned int i=0; i<maxImageIdx; i++)
232  {
233  m_CurrentStep = numberOfb0Images + i + 1;
234  if( m_AbortRegistration == true) {
235  m_IsInValidState = false;
236  mitkThrow() << "Stopped by user.";
237  };
238 
239  t_selector_w->SetTimeNr(i);
240  t_selector_w->Update();
241 
242  weightedRegistrationMethod->SetMovingImage( t_selector_w->GetOutput() );
243 
244  try
245  {
246  MITK_INFO << " === (" << i+1 <<"/"<< maxImageIdx << ") :: Starting registration";
247  weightedRegistrationMethod->Update();
248  }
249  catch( const itk::ExceptionObject& e)
250  {
251  m_IsInValidState = false;
252  mitkThrow() << "Failed to register the b0 images, the PyramidRegistration threw an exception: \n" << e.what();
253  }
254 
255  // allow expansion
256  mitk::ImageWriteAccessor imac(weightedRegistrationMethod->GetResampledMovingImage());
257  registeredWeighted->SetImportVolume( imac.GetData(),
258  i+1, 0, mitk::Image::CopyMemory);
259 
260  estimated_transforms.push_back( weightedRegistrationMethod->GetLastRotationMatrix() );
261  }
262 
263 
264  //
265  // (4) Cast the resulting image back to an diffusion weighted image
266  //
271  bzero_vector.fill(0);
272 
273  // compose the direction vector
274  // - no direction for the first image
275  // - correct ordering of the directions based on the index list
276  gradients_new->push_back( bzero_vector );
277 
278  SplitFilterType::IndexListType index_list = split_filter->GetIndexList();
279  SplitFilterType::IndexListType::const_iterator lIter = index_list.begin();
280 
281  while( lIter != index_list.end() )
282  {
283  gradients_new->push_back( gradients->at( *lIter ) );
284  ++lIter;
285  }
286 
287  registeredWeighted->SetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( gradients_new.GetPointer() ));
288  registeredWeighted->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( static_cast<mitk::FloatProperty*>(input->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ) );
289 
290  //
291  // (5) Adapt the gradient directions according to the estimated transforms
292  //
293  typedef mitk::DiffusionImageCorrectionFilter CorrectionFilterType;
295 
296  mitk::Image::Pointer output = registeredWeighted;
297  corrector->SetImage( output );
298  corrector->CorrectDirections( estimated_transforms );
299 
300  //
301  // (6) Pass the corrected image to the filters output port
302  //
303  m_CurrentStep += 1;
304 
305  this->GetOutput()->Initialize( output );
306 
307  this->GetOutput()->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( static_cast<mitk::GradientDirectionsProperty*>( output->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer() ) );
308  this->GetOutput()->SetProperty( mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str(), mitk::MeasurementFrameProperty::New( static_cast<mitk::MeasurementFrameProperty*>( output->GetProperty(mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str()).GetPointer() )->GetMeasurementFrame() ) );
309  this->GetOutput()->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( static_cast<mitk::FloatProperty*>(output->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ) );
310 
311  mitk::DiffusionPropertyHelper propertyHelper( this->GetOutput() );
312  propertyHelper.InitializeImage();
313  this->GetOutput()->Modified();
314 }
315 
316 #endif // MITKDWIHEADMOTIONCORRECTIONFILTER_CPP
static const std::string REFERENCEBVALUEPROPERTYNAME
itk::SmartPointer< Self > Pointer
#define MITK_INFO
Definition: mitkLogMacros.h:22
virtual BaseGeometry::Pointer GetGeometryForTimeStep(TimeStepType timeStep) const override
Returns the geometry which corresponds to the given time step.
This class has an advanced functionality to the B0ImageExtractionImageFilter, however the b0 images a...
virtual void GenerateData() override
A version of GenerateData() specific for image processing filters.
Helper class for mitk::Images containing diffusion weighted data.
The PyramidImageRegistration class implements a multi-scale registration method.
void * GetData()
Gives full data access.
static const std::string MEASUREMENTFRAMEPROPERTYNAME
GradientDirectionsProperty::GradientDirectionsContainerType GradientDirectionsContainerType
void InitializeImage()
Make certain the owned image is up to date with all necessary properties.
virtual void Expand(TimeStepType size) override
Expands the time geometry to the given number of time steps.
virtual void SetTimeStepGeometry(BaseGeometry *geometry, TimeStepType timeStep) override
Sets the geometry for the given time step.
Splits a DW-Image passed in as input into a 3D-t image where each volume coresponds to a gradient ima...
#define mitkThrow()
Image class for storing images.
Definition: mitkImage.h:76
GradientDirectionsProperty::GradientDirectionType GradientDirectionType
static Pointer New()
static Pointer New()
#define AccessFixedDimensionByItk_1(mitkImage, itkImageTypeFunction, dimension, arg1)
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
ImageWriteAccessor class to get locked write-access for a particular image part.
mitk::BaseProperty::Pointer GetProperty(const char *propertyKey) const
Get the property (instance of BaseProperty) with key propertyKey from the PropertyList, and set it to this, respectively;.
static const std::string GRADIENTCONTAINERPROPERTYNAME
static const std::string BVALUEMAPPROPERTYNAME
static Pointer New()
static void ItkAccumulateFilter(const itk::Image< TPixel, VDimensions > *image, mitk::Image::Pointer &output)
Averages an 3d+t image along the time axis.
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.