DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
SphericalPolarFourierImaging.cxx
Go to the documentation of this file.
1 
18 #include "SphericalPolarFourierImagingCLP.h"
19 
21 
23 #include "itkDWIReader.h"
30 #include "itkSamplingScheme3D.h"
31 
33 
34 #include "utl.h"
35 
39 int
40 main (int argc, char const* argv[])
41 {
42  // GenerateCLP
43  PARSE_ARGS;
44 
45  utl::LogLevel = _Verbose;
46  bool _Debug = _Verbose>=LOG_DEBUG;
47 
48  utlGlobalException(!_OutputEAPProfileFileArg.isSet() && !_OutputODFFileArg.isSet() && !_OutputSPFFileArg.isSet(), "no output");
49 
50  typedef double PrecisionType;
51  typedef itk::VectorImage<PrecisionType, 3> VectorImageType;
52  typedef itk::Image<PrecisionType, 3> ImageType;
53 
54  typedef itk::DWIReader<PrecisionType> DWIReaderType;
55 
56  ImageType::Pointer mdImage=NULL, maskImage=NULL;
57  if (_MDImageFileArg.isSet())
58  itk::ReadImage<ImageType>(_MDImageFile, mdImage);
59  if (_MaskFileArg.isSet())
60  itk::ReadImage<ImageType>(_MaskFile, maskImage);
61 
62  // read DWI, mdImage, mask
63  DWIReaderType::Pointer reader = DWIReaderType::New();
64  reader->GetSamplingSchemeQSpace()->SetTau(_Tau);
65  reader->SetConfigurationFile(_InputFile);
66  if (_MaskFileArg.isSet())
67  reader->SetMaskImage(maskImage);
68  reader->Update();
69 
70 
71 
72  // SPFI
73  VectorImageType::Pointer spf;
75  SPFIFilterBaseType::Pointer spfiFilter=NULL;
77  spfiFilter = SPFIFilterType::New();
78  std::cout << "Use SPF basis" << std::endl << std::flush;
79 
80  utl::InitializeThreadedLibraries(_NumberOfThreads);
81  if (_NumberOfThreads>0)
82  spfiFilter->SetNumberOfThreads(_NumberOfThreads);
83  if (_MaskFileArg.isSet())
84  spfiFilter->SetMaskImage(maskImage);
85  spfiFilter->SetInput(reader->GetOutput());
86  spfiFilter->SetSamplingSchemeQSpace(reader->GetSamplingSchemeQSpace());
87  spfiFilter->SetSHRank(_SHRank);
88  spfiFilter->SetRadialRank(_RadialRank);
89  spfiFilter->SetMD0(_MD0);
90  spfiFilter->SetBasisScale(_Scale);
91 
92  spfiFilter->SetIsAnalyticalB0(true);
93 
94  if (_EstimationType=="LS")
95  spfiFilter->SetEstimationType(SPFIFilterBaseType::LS);
96  else if (_EstimationType=="L1_2")
97  {
98  utlGlobalException(_LambdaSpherical<=0 && _LambdaRadial<=0, "need to set _LambdaSpherical and _LambdaRadial when _EstimationType=\"L1_2\".");
99  spfiFilter->SetEstimationType(SPFIFilterBaseType::L1_2);
100  }
101  else if (_EstimationType=="L1_DL")
102  {
103  utlGlobalException(_LambdaL1<0, "need to set _LambdaL1 when _EstimationType=\"L1_DL\".");
104  spfiFilter->SetEstimationType(SPFIFilterBaseType::L1_DL);
105  }
106  spfiFilter->SetLambdaSpherical(_LambdaSpherical);
107  spfiFilter->SetLambdaRadial(_LambdaRadial);
108  spfiFilter->SetLambdaL1(_LambdaL1);
109 
110  if (_SolverType=="FISTA_LS")
111  {
113  L1SolverType::Pointer l1Solver = L1SolverType::New();
114  l1Solver->SetUseL2SolverForInitialization(_SolverType=="FISTA_LS");
115  l1Solver->SetMaxNumberOfIterations(_MaxIter);
116  l1Solver->SetMinRelativeChangeOfCostFunction(_MinChange);
117  l1Solver->SetMinRelativeChangeOfPrimalResidual(_MinChange);
118  spfiFilter->SetL1FISTASolver(l1Solver);
119  spfiFilter->SetL1SolverType(SPFIFilterBaseType::FISTA_LS);
120  }
121  if (_SolverType=="SPAMS")
122  {
123  typedef itk::SpamsWeightedLassoSolver<double> L1SolverType;
124  L1SolverType::Pointer l1Solver = L1SolverType::New();
125  spfiFilter->SetL1SpamsSolver(l1Solver);
126  spfiFilter->SetL1SolverType(SPFIFilterBaseType::SPAMS);
127  }
128 
129 
130  typedef SPFIFilterBaseType::SamplingSchemeQSpaceType SamplingSchemeQSpaceType;
131  typedef SPFIFilterBaseType::SamplingSchemeRSpaceType SamplingSchemeRSpaceType;
132  SamplingSchemeQSpaceType::Pointer samplingQSpace = SamplingSchemeQSpaceType::New();
133  SamplingSchemeRSpaceType::Pointer samplingRSpace = SamplingSchemeRSpaceType::New();
134 
135  if (_MDImageFileArg.isSet())
136  spfiFilter->SetMDImage(mdImage);
137 
138  spfiFilter->SetBasisEnergyPowerDL(1.0);
139 
141  if (_ShowProgressArg.isSet())
142  spfiFilter->AddObserver( itk::ProgressEvent(), observer );
143  spfiFilter->SetDebug(_Verbose>=LOG_DEBUG);
144  spfiFilter->SetLogLevel(utl::LogLevel);
145 
146  std::cout << "SPF estimation starts" << std::endl << std::flush;
147  spfiFilter->Update();
148  std::cout << "SPF estimation ends" << std::endl << std::flush;
149  spf = spfiFilter->GetOutput();
150  // save SPF coefficients if needed
151  if (_OutputSPFFileArg.isSet())
152  {
153  itk::SaveImage<VectorImageType>(spf, _OutputSPFFile);
154  }
155 
156  // output features
158  // FeaturesFromSPFFilterType::Pointer featureFromSPFFilter=0;
159 
160  // output EAP
161  if (_OutputEAPProfileFileArg.isSet())
162  {
164  ProfileFromSPFFilterType::Pointer featureFromSPFFilter = ProfileFromSPFFilterType::New();
165  if (_MaskFileArg.isSet())
166  featureFromSPFFilter->SetMaskImage(maskImage);
167  if (_NumberOfThreads>0)
168  featureFromSPFFilter->SetNumberOfThreads(_NumberOfThreads);
169  featureFromSPFFilter->SetSHRank(_SHRank);
170  featureFromSPFFilter->SetRadialRank(_RadialRank);
171  featureFromSPFFilter->SetMD0(spfiFilter->GetMD0());
172  featureFromSPFFilter->SetTau(spfiFilter->GetSamplingSchemeQSpace()->GetTau());
173  featureFromSPFFilter->SetBasisScale(spfiFilter->GetBasisScale());
174  featureFromSPFFilter->SetBasisType(FeaturesFromSPFFilterType::SPF);
175  if (_MDImageFileArg.isSet())
176  {
177  featureFromSPFFilter->SetScaleImage(spfiFilter->GetScaleImage());
178  }
179  if (_Debug)
180  featureFromSPFFilter->DebugOn();
181  featureFromSPFFilter->SetInput(spf);
182  featureFromSPFFilter->SetRadius(_Radius);
183  featureFromSPFFilter->SetIsFourier(true);
184  featureFromSPFFilter->SetIsInQSpace(true);
185  std::cout << "EAP profile estimation starts" << std::endl << std::flush;
186  featureFromSPFFilter->Update();
187  std::cout << "EAP profile estimation ends" << std::endl << std::flush;
188  VectorImageType::Pointer eap = featureFromSPFFilter->GetOutput();
189  itk::SaveImage<VectorImageType>(eap, _OutputEAPProfileFile);
190  }
191 
192  // output ODF
193  if (_OutputODFFileArg.isSet())
194  {
196  ODFFromSPFFilterType::Pointer featureFromSPFFilter = ODFFromSPFFilterType::New();
197  if (_MaskFileArg.isSet())
198  featureFromSPFFilter->SetMaskImage(maskImage);
199  if (_NumberOfThreads>0)
200  featureFromSPFFilter->SetNumberOfThreads(_NumberOfThreads);
201  featureFromSPFFilter->SetSHRank(_SHRank);
202  featureFromSPFFilter->SetRadialRank(_RadialRank);
203  featureFromSPFFilter->SetMD0(spfiFilter->GetMD0());
204  featureFromSPFFilter->SetTau(spfiFilter->GetSamplingSchemeQSpace()->GetTau());
205  featureFromSPFFilter->SetBasisScale(spfiFilter->GetBasisScale());
206  featureFromSPFFilter->SetBasisType(FeaturesFromSPFFilterType::SPF);
207  if (_MDImageFileArg.isSet())
208  {
209  featureFromSPFFilter->SetScaleImage(spfiFilter->GetScaleImage());
210  }
211  if (_Debug)
212  featureFromSPFFilter->DebugOn();
213  featureFromSPFFilter->SetInput(spf);
214  featureFromSPFFilter->SetODFOrder(_ODFOrder);
215  featureFromSPFFilter->SetBMax(-1.0);
216  featureFromSPFFilter->SetIsFourier(true);
217  featureFromSPFFilter->SetIsInQSpace(true);
218  std::cout << "ODF estimation starts" << std::endl << std::flush;
219  featureFromSPFFilter->Update();
220  std::cout << "ODF estimation ends" << std::endl << std::flush;
221  VectorImageType::Pointer odf = featureFromSPFFilter->GetOutput();
222  itk::SaveImage<VectorImageType>(odf, _OutputODFFile);
223  }
224 
225  return 0;
226 }
227 
helper functions specifically used in dmritool
Calculate DWI/EAP profile from SPF coefficients.
solve least square problem with L1 regularization using FISTA
used for debug information. this->GetDebug()
Definition: utlCoreMacro.h:193
calculate ODFs from SPF coefficients
#define utlGlobalException(cond, expout)
Definition: utlCoreMacro.h:372
static int LogLevel
Definition: utlCoreMacro.h:203
Load gradient file, b values and DWI files (with optional index file)
Definition: itkDWIReader.h:59
static void InitializeThreadedLibraries(const int numThreads)
Definition: utl.h:327
estimate the coeffcients of generalized Spherical Polar Fourier basis which can be separated into dif...
solve weighted LASSO using spams
itk::VectorImage< ScalarType, 3 > VectorImageType
Definition: 4DImageMath.cxx:30
Compute some features (DWI/EAP profile, ODFs, scalar indices) from SPF coefficients.
int main(int argc, char const *argv[])
general SPFI using SPF basis