Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
QballReconstruction.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 #include <mitkCoreObjectFactory.h>
18 #include "mitkImage.h"
19 #include "itkAnalyticalDiffusionQballReconstructionImageFilter.h"
20 #include <boost/lexical_cast.hpp>
21 #include "mitkCommandLineParser.h"
22 #include <mitkIOUtil.h>
23 #include <itksys/SystemTools.hxx>
24 #include <mitkDiffusionPropertyHelper.h>
25 #include <mitkITKImageImport.h>
26 #include <mitkImageCast.h>
27 #include <mitkProperties.h>
28 #include <mitkIOUtil.h>
29 
30 using namespace mitk;
31 using namespace std;
32 
36 int main(int argc, char* argv[])
37 {
38  mitkCommandLineParser parser;
39  parser.setArgumentPrefix("--", "-");
40  parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input file", "input raw dwi (.dwi or .fsl/.fslgz)", us::Any(), false);
41  parser.addArgument("outFile", "o", mitkCommandLineParser::OutputFile, "Output file", "output file", us::Any(), false);
42  parser.addArgument("shOrder", "sh", mitkCommandLineParser::Int, "Spherical harmonics order", "spherical harmonics order", 4, true);
43  parser.addArgument("b0Threshold", "t", mitkCommandLineParser::Int, "b0 threshold", "baseline image intensity threshold", 0, true);
44  parser.addArgument("lambda", "r", mitkCommandLineParser::Float, "Lambda", "ragularization factor lambda", 0.006, true);
45  parser.addArgument("csa", "csa", mitkCommandLineParser::Bool, "Constant solid angle consideration", "use constant solid angle consideration");
46  parser.addArgument("outputCoeffs", "shc", mitkCommandLineParser::Bool, "Output coefficients", "output file containing the SH coefficients");
47  parser.addArgument("mrtrix", "mb", mitkCommandLineParser::Bool, "MRtrix", "use MRtrix compatible spherical harmonics definition");
48 
49  parser.setCategory("Preprocessing Tools");
50  parser.setTitle("Qball Reconstruction");
51  parser.setDescription("");
52  parser.setContributor("MBI");
53 
54  map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
55  if (parsedArgs.size()==0)
56  return EXIT_FAILURE;
57 
58  std::string inFileName = us::any_cast<string>(parsedArgs["input"]);
59  std::string outfilename = us::any_cast<string>(parsedArgs["outFile"]);
60  outfilename = itksys::SystemTools::GetFilenamePath(outfilename)+"/"+itksys::SystemTools::GetFilenameWithoutExtension(outfilename);
61 
62  int threshold = 0;
63  if (parsedArgs.count("b0Threshold"))
64  threshold = us::any_cast<int>(parsedArgs["b0Threshold"]);
65 
66  int shOrder = 4;
67  if (parsedArgs.count("shOrder"))
68  shOrder = us::any_cast<int>(parsedArgs["shOrder"]);
69 
70  float lambda = 0.006;
71  if (parsedArgs.count("lambda"))
72  lambda = us::any_cast<float>(parsedArgs["lambda"]);
73 
74  int normalization = 0;
75  if (parsedArgs.count("csa") && us::any_cast<bool>(parsedArgs["csa"]))
76  normalization = 6;
77 
78  bool outCoeffs = false;
79  if (parsedArgs.count("outputCoeffs"))
80  outCoeffs = us::any_cast<bool>(parsedArgs["outputCoeffs"]);
81 
82  bool mrTrix = false;
83  if (parsedArgs.count("mrtrix"))
84  mrTrix = us::any_cast<bool>(parsedArgs["mrtrix"]);
85 
86  try
87  {
88  std::vector<BaseData::Pointer> infile = mitk::IOUtil::Load(inFileName);
89  Image::Pointer dwi = dynamic_cast<Image*>(infile.at(0).GetPointer());
90  mitk::DiffusionPropertyHelper propertyHelper(dwi);
91  propertyHelper.AverageRedundantGradients(0.001);
92  propertyHelper.InitializeImage();
93 
94  mitk::QBallImage::Pointer image = mitk::QBallImage::New();
95  mitk::Image::Pointer coeffsImage = mitk::Image::New();
96 
97  std::cout << "SH order: " << shOrder;
98  std::cout << "lambda: " << lambda;
99  std::cout << "B0 threshold: " << threshold;
100  switch ( shOrder )
101  {
102  case 4:
103  {
104  typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,4,QBALL_ODFSIZE> FilterType;
105  mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
106  mitk::CastToItkImage(dwi, itkVectorImagePointer);
107 
108  FilterType::Pointer filter = FilterType::New();
109  filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer );
110  filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi));
111  filter->SetThreshold( threshold );
112  filter->SetLambda(lambda);
113  filter->SetUseMrtrixBasis(mrTrix);
114  if (normalization==0)
115  filter->SetNormalizationMethod(FilterType::QBAR_STANDARD);
116  else
117  filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE);
118  filter->Update();
119  image->InitializeByItk( filter->GetOutput() );
120  image->SetVolume( filter->GetOutput()->GetBufferPointer() );
121  coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() );
122  coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() );
123  break;
124  }
125  case 6:
126  {
127  typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,6,QBALL_ODFSIZE> FilterType;
128  mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
129  mitk::CastToItkImage(dwi, itkVectorImagePointer);
130 
131  FilterType::Pointer filter = FilterType::New();
132  filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer );
133  filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi));
134  filter->SetThreshold( threshold );
135  filter->SetLambda(lambda);
136  filter->SetUseMrtrixBasis(mrTrix);
137  if (normalization==0)
138  filter->SetNormalizationMethod(FilterType::QBAR_STANDARD);
139  else
140  filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE);
141  filter->Update();
142  image->InitializeByItk( filter->GetOutput() );
143  image->SetVolume( filter->GetOutput()->GetBufferPointer() );
144  coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() );
145  coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() );
146  break;
147  }
148  case 8:
149  {
150  typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,8,QBALL_ODFSIZE> FilterType;
151  mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
152  mitk::CastToItkImage(dwi, itkVectorImagePointer);
153 
154  FilterType::Pointer filter = FilterType::New();
155  filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer );
156  filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi));
157  filter->SetThreshold( threshold );
158  filter->SetLambda(lambda);
159  filter->SetUseMrtrixBasis(mrTrix);
160  if (normalization==0)
161  filter->SetNormalizationMethod(FilterType::QBAR_STANDARD);
162  else
163  filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE);
164  filter->Update();
165  image->InitializeByItk( filter->GetOutput() );
166  image->SetVolume( filter->GetOutput()->GetBufferPointer() );
167  coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() );
168  coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() );
169  break;
170  }
171  case 10:
172  {
173  typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,10,QBALL_ODFSIZE> FilterType;
174  mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
175  mitk::CastToItkImage(dwi, itkVectorImagePointer);
176 
177  FilterType::Pointer filter = FilterType::New();
178  filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer );
179  filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi));
180  filter->SetThreshold( threshold );
181  filter->SetLambda(lambda);
182  filter->SetUseMrtrixBasis(mrTrix);
183  if (normalization==0)
184  filter->SetNormalizationMethod(FilterType::QBAR_STANDARD);
185  else
186  filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE);
187  filter->Update();
188  image->InitializeByItk( filter->GetOutput() );
189  image->SetVolume( filter->GetOutput()->GetBufferPointer() );
190  coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() );
191  coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() );
192  break;
193  }
194  case 12:
195  {
196  typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,12,QBALL_ODFSIZE> FilterType;
197  mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
198  mitk::CastToItkImage(dwi, itkVectorImagePointer);
199 
200  FilterType::Pointer filter = FilterType::New();
201  filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer );
202  filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi));
203  filter->SetThreshold( threshold );
204  filter->SetLambda(lambda);
205  if (normalization==0)
206  filter->SetNormalizationMethod(FilterType::QBAR_STANDARD);
207  else
208  filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE);
209  filter->Update();
210  image->InitializeByItk( filter->GetOutput() );
211  image->SetVolume( filter->GetOutput()->GetBufferPointer() );
212  coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() );
213  coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() );
214  break;
215  }
216  default:
217  {
218  std::cout << "Supplied SH order not supported. Using default order of 4.";
219  typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,4,QBALL_ODFSIZE> FilterType;
220  mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
221  mitk::CastToItkImage(dwi, itkVectorImagePointer);
222 
223  FilterType::Pointer filter = FilterType::New();
224  filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer );
225  filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi));
226  filter->SetThreshold( threshold );
227  filter->SetLambda(lambda);
228  filter->SetUseMrtrixBasis(mrTrix);
229  if (normalization==0)
230  filter->SetNormalizationMethod(FilterType::QBAR_STANDARD);
231  else
232  filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE);
233  filter->Update();
234  image->InitializeByItk( filter->GetOutput() );
235  image->SetVolume( filter->GetOutput()->GetBufferPointer() );
236  coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() );
237  coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() );
238  }
239  }
240 
241  std::string coeffout = outfilename;
242  coeffout += "_shcoeffs.nrrd";
243 
244  outfilename += ".qbi";
245  mitk::IOUtil::SaveBaseData(image, outfilename);
246 
247  if (outCoeffs)
248  mitk::IOUtil::SaveImage(coeffsImage, coeffout);
249  }
250  catch ( itk::ExceptionObject &err)
251  {
252  std::cout << "Exception: " << err;
253  }
254  catch ( std::exception err)
255  {
256  std::cout << "Exception: " << err.what();
257  }
258  catch ( ... )
259  {
260  std::cout << "Exception!";
261  }
262  return EXIT_SUCCESS;
263 }
mitk::DiffusionPropertyHelper::AverageRedundantGradients
void AverageRedundantGradients(double precision)
Definition: mitkDiffusionPropertyHelper.cpp:84
mitk::Pointer
itk::SmartPointer< Self > Pointer
Definition: mitkRenderingManager.h:389
mitkCommandLineParser::setContributor
void setContributor(std::string contributor)
Definition: mitkCommandLineParser.cpp:938
std
STL namespace.
mitk::DiffusionPropertyHelper
Helper class for mitk::Images containing diffusion weighted data.
Definition: mitkDiffusionPropertyHelper.h:37
mitk
DataCollection - Class to facilitate loading/accessing structured data.
Definition: GeometryOverview.dox:1
us::any_cast
ValueType * any_cast(Any *operand)
Definition: usAny.h:377
main
int main(int argc, char *argv[])
Definition: QballReconstruction.cpp:36
mitkCommandLineParser::parseArguments
std::map< std::string, us::Any > parseArguments(const StringContainerType &arguments, bool *ok=nullptr)
Definition: mitkCommandLineParser.cpp:355
mitk::QBallImage::New
static Pointer New()
mitk::DiffusionPropertyHelper::InitializeImage
void InitializeImage()
Make certain the owned image is up to date with all necessary properties.
Definition: mitkDiffusionPropertyHelper.cpp:285
mitk::IOUtil::SaveImage
static bool SaveImage(mitk::Image::Pointer image, const std::string &path)
SaveImage Convenience method to save an arbitrary mitkImage.
Definition: mitkIOUtil.cpp:870
mitk::IOUtil::SaveBaseData
static bool SaveBaseData(mitk::BaseData *data, const std::string &path)
SaveBaseData Convenience method to save arbitrary baseData.
Definition: mitkIOUtil.cpp:888
itk::AnalyticalDiffusionQballReconstructionImageFilter
This class takes as input one or more reference image (acquired in the absence of diffusion sensitizi...
Definition: itkAnalyticalDiffusionQballReconstructionImageFilter.h:91
mitkCommandLineParser::addArgument
void addArgument(const std::string &longarg, const std::string &shortarg, Type type, const std::string &argLabel, const std::string &argHelp=std::string(), const us::Any &defaultValue=us::Any(), bool optional=true, bool ignoreRest=false, bool deprecated=false)
Definition: mitkCommandLineParser.cpp:664
mitk::Image
Image class for storing images.
Definition: mitkImage.h:76
us::Any
Definition: usAny.h:163
mitkCommandLineParser::setCategory
void setCategory(std::string category)
Definition: mitkCommandLineParser.cpp:943
mitk::Image::New
static Pointer New()
mitkCommandLineParser::setArgumentPrefix
void setArgumentPrefix(const std::string &longPrefix, const std::string &shortPrefix)
Definition: mitkCommandLineParser.cpp:823
mitk::CastToItkImage
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
mitk::DiffusionPropertyHelper::GetGradientContainer
GradientDirectionsContainerType::Pointer GetGradientContainer() const
Definition: mitkDiffusionPropertyHelper.cpp:418
mitk::IOUtil::Load
static DataStorage::SetOfObjects::Pointer Load(const std::string &path, DataStorage &storage)
Load a file into the given DataStorage.
Definition: mitkIOUtil.cpp:483
mitkCommandLineParser::setTitle
void setTitle(std::string title)
Definition: mitkCommandLineParser.cpp:934
mitkCommandLineParser::setDescription
void setDescription(std::string description)
Definition: mitkCommandLineParser.cpp:948
mitk::New
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.