An open framework for parameter mapping and image processing

1. An open framework for parameter mapping and image processing Apr 8, 2013 Jason Su

2. Motivation • Mapping methods share many common tasks and features • Fitting of a tissue model to data with a known signal equation • Multi-subject processing • Residuals and statistics analysis to evaluate fit • Pause and resume, recovery from crash • We should promote code sharing in the community, so groups around the world can experiment with promising new mapping methods

3. Goals • There are 2 audiences, coders and users • For developers of novel methods: • Provide a framework that fosters readability and code re-use • Easy to convert a bare-bones implementation into production-level functionality, avoid code re-writes • Provide easy ways to speed up code: embarrassing parallelization (across slices or subjects) and caching • What else do developers need? • Data/error visualization • ?

4. Goals • For clinicians and users: • A consistent workflow and usage pattern for producing a multitude of mapping methods • Support for all platforms • Scalability to large-subject studies and diverse computing resources • Amazon EC2 has built-in support for Python • A community • Open source repository to share and improve

5. The Framework: Protocol class • Like when we sit down at the scanner, we start with a protocol • Input is readable JSON, not an inscrutable stream of command-line switches • Instantiates a bunch of Series and manages access to their image data • I’ve assumed voxel-wise fits { "Mask": { "filename": “brainMask_mask.nii.gz" }, "SPGR": [ { "filename”: "spgr_fa3000.nii.gz", "flip_angle": 3.0, "repetition_time": 6.74, "weight": 2 }, ... { "filename": "spgr_fa8000.nii.gz", "flip_angle": 8.0, "repetition_time": 6.74 } ], "SSFP": [ { "filename": "ssfp_fa11000_ph0.nii.gz", "flip_angle": 11.0, "repetition_time": 3.6, "rf_phase_cycle": 0.0, "weight": 2 }, ... ] }

6. Series class • Defines important scan parameters, can pull from DICOM header • Loads image from disk • SPGR and SSFP are implemented • Flexible for any combination of phase cycles, TRs, flip angles • B0, B1/Kappa, ROI?

7. Utilities module • Memoization/caching of expensive results • pandas data structures inspired by R • Vectors and table data frames with named access • Notification of other functions when data changes >> @memo def stats(x): return (mean(x), std(x)) >> x = DataFrame(d, index=['a', 'b', 'c', 'd']) P001 P002 Types 'RRMS''SPMS' EDSS 3.5 6.0 >> x['Bob'] >> y = notify_container_factory( DataFrame)(stats, data=x) >> y.ix['EDSS'] = [4.0 7.0] -> stats() function is called and recomputes mean and std. of EDSS

8. Implementing a signal equation • Define a TissueModel • RelaxationTissueModel • Take parameters from TissueModel and compute expensive intermediate parts of the signal equation: • SPGRSystemArray and SSFPSystemArray • Finally, incorporate sequence parameters from a Series and compute the intensity at a voxel • SPGRSignal and SSFPSignal • Series are able to discover the modules that model their signal behavior • Automatically looks for the module based on its name • SPGRSeries is tied to the sequence.spgr module

9. TissueModel class • A bare-bones base class that provides a flexible dictionary hierarchy and notification > model.parameters['myelin']['fraction'] = 0.2 • Methods: • add_notifiee • del_notifiee • desync – do something when the model is changed • sync – call notifiees, propagate the changes • Notifiees are called in the order they are added, duplicates are ignored • Acts like a state machine • Classes group/organize functions and the data they manipulate under common umbrellas

10. RelaxationTissueModel • Converts parameters to pandas DataFrame for vectorized operations • N-component relaxation model • mcDESPOT is usually 2-component • Pure tissue, does not include terms like RF phase cycle • desync enforces model assumptions • Chemical equilibrium and unit total volume component_list = [ { 'name': 'intraextra', 'T1': 1000, 'T2': 70, 'off-resonance': 0.1, 'fraction': 0.75, 'exchange_rate': {'myelin': .1, 'csf': 0 } }, { 'name': 'myelin', 'T1': 800.0, 'T2': 30.0, 'off-resonance': 0.2, 'fraction': 0.2, 'exchange_rate': {'intraextra': .1, 'csf': 0 } }, { 'name': 'csf', 'T1': 5000.0, 'T2': 300.0, 'off-resonance': 0.3, 'fraction': 0.05, 'exchange_rate': {'intraextra': 0, 'myelin': 0 } } ] model = RelaxationTissueModel(component_list)

11. Simplifying Assumptions • 2 component model • Only need to find fF (the fast volume fraction) • Chemical equilibrium • Allows us to eliminate finding kSF=1/τS • Both components are on the same resonance

12. SPGR/SSFPSystemArray • Subclassing promotes code re-use • SPGR is the simpler form, so SPGR extends it • Computes shared intermediate quantities • A shared cache between different objects • Notified by TissueModel

13. SPGR/SSFPSignal • With phase cycle separated out of the tissue model matrix, matrix multiply instead of expm • Adding more phase cycles has minimal impact on processing time • Current C implementation calculates all quantities for every image

14. Solvers • Input a TissueModel, constraints, and initial guess to fit • Each of these are like TissueModel • Stochastic Region Contraction, in progress • Store summary about each contraction step for later viewing • Idea: chain a solver like Newton’s or gradient descent once SRC has narrowed in

15. What is the advantage of all this structure? • Thinking about code re-use and class structure helps to expose similarities in signal equations • Can lead to speed-ups as with SPGR/SSFP • Improved debugging via unit testing and exceptions • e.g. Series raises exceptions giving valuable feedback about what sequence parameters are missing in the input JSON • Dynamic typing can allow for unexpected clever use by others

16. Open Source: Python vs. Matlab • Matlab Central is a nice central repository for interesting code • But for the most part it is single contributor • Limits the possibility of large scale improvements to MATLAB functionality • Python Package Index (PyPI) • Has single line install for most everything • Just in the category of nearly drop-in acceleration of Python: • Theano – GPU and generated C • Numba(Pro) – LLVM • PyPy/NumPyPy – total rewrite of CPython core interpreter • Numexpr – optimizes evaluation of math expressions including considerations like CPU cache

17. Disadvantages • Need to learn a new language • Loss of existing tools in MATLAB, including: • read_pfile • Extended phase graph analysis • Development environment is different • Spyder tries to imitate MATLAB IDE • Too many packages to install? • There are bundles that have all of the best packages for Python with a double-click

18. Possible Future Features • Access to super-large files with HDF5 and PyTables, important for high-res 32ch data • Extending to handle more exotic cases • QSM • ? • P-files and raw data? • Reconstruction does share some similar desired functions: • Multi-subject processing, pause and resume • Auto-suggest of possible fits given a Protocol • GUI for users?