@@ -160,6 +160,126 @@ def _gen_outfilename(self):
160160 _ , name , _ = split_filename (self .inputs .scheme_file )
161161 return name + '_recondata.Bdouble'
162162
163+ class MESDInputSpec (StdOutCommandLineInputSpec ):
164+ in_file = File (exists = True , argstr = '-inputfile %s' , mandatory = True , position = 1 ,
165+ desc = 'voxel-order data filename' )
166+ inverter = traits .Enum ('SPIKE' , 'PAS' , argstr = '-filter %s' , position = 2 , mandatory = True ,
167+ desc = ('The inversion index specifies the type of inversion to perform on the data.'
168+ 'The currently available choices are:'
169+ 'Inverter name | Inverter parameters'
170+ '---------------|------------------'
171+ 'SPIKE | bd (b-value x diffusivity along the fibre.)'
172+ 'PAS | r' ))
173+ inverter_param = traits .Float (argstr = '%f' , units = 'NA' , position = 3 , mandatory = True ,
174+ desc = ('Parameter associated with the inverter. Cf. inverter description for'
175+ 'more information.' ))
176+ fastmesd = traits .Bool (argstr = '-fastmesd' , requires = ['mepointset' ],
177+ desc = ('Turns off numerical integration checks and fixes the integration point set size at that of'
178+ 'the index specified by -basepointset..' ))
179+ mepointset = traits .Int (argstr = '-mepointset %d' , units = 'NA' ,
180+ desc = ('Use a set of directions other than those in the scheme file for the deconvolution kernel.'
181+ 'The number refers to the number of directions on the unit sphere. For example, '
182+ '"-mepointset 54" uses the directions in "camino/PointSets/Elec054.txt".' ))
183+ scheme_file = File (exists = True , argstr = '-schemefile %s' , mandatory = True ,
184+ desc = 'Specifies the scheme file for the diffusion MRI data' )
185+ bgmask = File (exists = True , argstr = '-bgmask %s' , desc = 'background mask' )
186+ inputdatatype = traits .Enum ('float' , 'char' , 'short' , 'int' , 'long' , 'double' , argstr = '-inputdatatype %s' ,
187+ desc = ('Specifies the data type of the input file: "char", "short", "int", "long",'
188+ '"float" or "double". The input file must have BIG-ENDIAN ordering.'
189+ 'By default, the input type is "float".' ))
190+
191+ class MESDOutputSpec (TraitedSpec ):
192+ mesd_data = File (exists = True , desc = 'MESD data' )
193+
194+ class MESD (StdOutCommandLine ):
195+ """
196+ MESD is a general program for maximum entropy spherical deconvolution.
197+ It also runs PASMRI, which is a special case of spherical deconvolution.
198+ The input data must be in voxel order.
199+
200+ The format of the output in each voxel is:
201+ { exitcode, ln(A^star(0)), lambda_0, lambda_1, ..., lambda_N }
202+
203+ The exitcode contains the results of three tests. The first test thresholds
204+ the maximum relative error between the numerical integrals computed at con-
205+ vergence and those computed using a larger test point set; if the error is
206+ greater than a threshold the exitcode is increased from zero to one as a
207+ warning; if it is greater than a larger threshold the exitcode is increased to
208+ two to suggest failure. The second test thresholds the predicted error in
209+ numerical integrals computed using the test point set; if the predicted error
210+ is greater than a threshold the exitcode is increased by 10. The third test
211+ thresholds the RMS error between the measurements and their predictions from
212+ the fitted deconvolution; if the errors are greater than a threshold, the exit
213+ code is increased by 100. An exitcode of 112 means that all three tests were
214+ failed and the result is likely to be unreliable. If all is well the exitcode
215+ is zero. Results are often still reliable even if one or two of the tests are
216+ failed.
217+
218+ Other possible exitcodes are:
219+ 5 - The optimization failed to converge
220+ -1 - Background
221+ -100 - Something wrong in the MRI data, e.g. negative or zero measurements,
222+ so that the optimization could not run.
223+
224+ The standard MESD implementation is computationally demanding, particularly
225+ as the number of measurements increases (computation is approximately O(N^2),
226+ where N is the number of measurements). There are two ways to obtain significant
227+ computational speed-up:
228+
229+ i) Turn off error checks and use a small point set for computing numerical
230+ integrals in the algorithm by adding the flag -fastmesd. Sakaie CDMRI 2008
231+ shows that using the smallest point set (-basepointset 0) with no
232+ error checks usually has only a minor effect on the output of the algorithm,
233+ but provides a major reduction in computation time. You can increase the point
234+ set size using -basepointset with an argument higher than 0, which may produce
235+ better results in some voxels, but will increase computation time, which
236+ approximately doubles every time the point set index increases by 1.
237+
238+ ii) Reduce the complexity of the maximum entropy encoding using -mepointset <X>.
239+ By default <X> = N, the number of measurements, and is the number of parameters
240+ in the max. ent. representation of the output function, ie the number of
241+ lambda parameters, as described in Jansons and Alexander Inverse Problems 2003.
242+ However, we can represent the function using less components and <X> here
243+ specifies the number of lambda parameters. To obtain speed-up, set <X>
244+ < N; complexity become O(<X>^2) rather than O(N^2). Note that <X> must be chosen
245+ so that the camino/PointSets directory contains a point set with that number
246+ of elements. When -mepointset decreases, the numerical integration checks
247+ make less and less of a difference and smaller point sets for numerical
248+ integration (see -basepointset) become adequate. So when <X> is low -fastmesd is
249+ worth using to get even more speed-up.
250+
251+ The choice of <X> is a parameter of the technique. Too low and you lose angular
252+ resoloution; too high and you see no computational benefit and may even suffer
253+ from overfitting. Empirically, we have found that <X>=16 often gives good
254+ results and good speed up, but it is worth trying a few values a comparing
255+ performance. The reduced encoding is described in the following ISMRM abstract:
256+ Sweet and Alexander "Reduced Encoding Persistent Angular Structure" 572 ISMRM 2010.
257+
258+ Example
259+ ---------
260+ Run MESD on every voxel of the data file SubjectA.Bfloat using the PASMRI kernel.
261+
262+ >>> import nipype.interfaces.camino as cam
263+ >>> mesd = cam.MESD()
264+ >>> mesd.inputs.in_file = 'SubjectA.Bfloat'
265+ >>> mesd.inputs.scheme_file = 'A.scheme'
266+ >>> mesd.inputs.inverter = 'PAS'
267+ >>> mesd.inputs.inverter_param = 1.4
268+ >>> mesd.run() # doctest: +SKIP
269+ """
270+ _cmd = 'mesd'
271+ input_spec = MESDInputSpec
272+ output_spec = MESDOutputSpec
273+
274+ def _list_outputs (self ):
275+ outputs = self .output_spec ().get ()
276+ outputs ['mesd_data' ] = os .path .abspath (self ._gen_outfilename ())
277+ return outputs
278+
279+ def _gen_outfilename (self ):
280+ _ , name , _ = split_filename (self .inputs .scheme_file )
281+ return name + '_MESD.Bdouble'
282+
163283class SFPeaksInputSpec (StdOutCommandLineInputSpec ):
164284 in_file = File (exists = True , argstr = '-inputfile %s' , mandatory = True ,
165285 desc = 'Voxel-order data of spherical functions' )
0 commit comments