Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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.