Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merge SHARDRecon into main MRtrix3 codebase #2994

Open
wants to merge 434 commits into
base: dev
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 250 commits
Commits
Show all changes
434 commits
Select commit Hold shift + click to select a range
38795b9
Merge branch 'master' into tvreg2
dchristiaens Nov 20, 2017
52dc716
motionfilter: revert to median filtering.
dchristiaens Nov 20, 2017
cbecd53
motionfilter: support volume-level transformations.
dchristiaens Nov 20, 2017
1e6ec97
Integrate motion filter in motion correction script. Parameters now u…
dchristiaens Nov 20, 2017
1ef3bd9
mssh2amp: fix output gradient table at 4 columns.
dchristiaens Nov 21, 2017
917be84
dwimotioncorrect: increment SVR progress bar less often.
dchristiaens Nov 22, 2017
9f248d5
motionfilter: add unused options to support mrtrix script multi-threa…
dchristiaens Nov 24, 2017
21c4ea0
Fix cubic interpolation kernel.
dchristiaens Dec 2, 2017
79f8021
dwimotioncorrect: blur recon for registration.
dchristiaens Dec 5, 2017
824c23c
motionfilter: catch numeric errors in matrix logarithm.
dchristiaens Dec 8, 2017
9343d3d
dwimotioncorrect: add optional intensity matching.
dchristiaens Dec 11, 2017
f23c947
dwirecon: small bug fix in jacobian modulation.
dchristiaens Dec 11, 2017
86ade6a
dwirecon: exploit multiband factor for efficient q-space basis evalua…
dchristiaens Dec 12, 2017
467f538
dwirecon: bug fix in motion lookup and symbol renaming
dchristiaens Dec 12, 2017
e45b6b4
motionfilter: adapt to new motion file format.
dchristiaens Dec 12, 2017
f170a81
dwirecon: bug fix in MB loop indexing.
dchristiaens Dec 12, 2017
a4ce2a8
mrfieldunwarp: support motion tables of any multiband dimension.
dchristiaens Dec 12, 2017
0b3ff55
dwirecon: print progress at info level.
dchristiaens Dec 19, 2017
27e5f55
dwirecon: code cleanup.
dchristiaens Dec 19, 2017
75e67fc
Working first prototype for Lie-algebra rigid registration.
dchristiaens Dec 19, 2017
9e9709b
mrreg_asym: Fix Jacobian scaling.
dchristiaens Dec 19, 2017
6b9a2f5
mrreg_asym: add mask.
dchristiaens Dec 20, 2017
1a36b3d
mrreg_asym: code organization.
dchristiaens Dec 20, 2017
489ab4b
mrreg_asym: interface improvements.
dchristiaens Dec 20, 2017
07ff722
mrreg_asym: add user option for starting point.
dchristiaens Dec 20, 2017
ffbee9a
Framework for new slice-to-volume registrtaion command.
dchristiaens Dec 21, 2017
8eb1506
dwislicealign: set up inrastructure for s2v registration.
dchristiaens Dec 21, 2017
7b4b4c2
dwislicealign: implement registration.
dchristiaens Dec 23, 2017
4f835de
dwislicealign: bug fixes
dchristiaens Jan 2, 2018
c33dc18
dwislicealign: add progress bar.
dchristiaens Jan 2, 2018
e7dcd62
dwimotioncorrect: integrate new registration command
dchristiaens Jan 2, 2018
af78fc8
dwislicealign: critical bug fix
dchristiaens Jan 2, 2018
17533dd
mrfieldunwarp: add optional index for field alignment.
dchristiaens Jan 3, 2018
224ef2b
dwimotioncorrect: cleanup
dchristiaens Jan 3, 2018
7b34a79
dwisliceoutliers: use MB factor.
dchristiaens Jan 3, 2018
993f854
dwimotioncorrect: reinstate multi-scale registration
dchristiaens Jan 3, 2018
c13c619
dwisliceoutliers: fix bug.
dchristiaens Jan 3, 2018
7c6fa72
Merge branch 'master' into customreg
dchristiaens Jan 3, 2018
803f99f
dwirecon: print residual error on console.
dchristiaens Jan 4, 2018
94dee12
Merge branch 'customreg' into intensitymatching
dchristiaens Jan 4, 2018
e282b36
intensity matching: use rank-reduced prediction and mask
dchristiaens Jan 9, 2018
172f2eb
intensity matching: separate outlier reweighting and bias field scali…
dchristiaens Jan 9, 2018
65ad2f0
dwislicealign: incorporate SSP in registration.
dchristiaens Jan 9, 2018
818121a
dwimotioncorrect: fix multi-level registration.
dchristiaens Jan 9, 2018
b892d2c
Merge branch 'customreg' into intensitymatching2
dchristiaens Jan 11, 2018
76c2716
dwirecon: use custom SSP.
dchristiaens Jan 15, 2018
3b289e0
Merge branch 'master' into customreg
dchristiaens Jan 15, 2018
2827ca7
dwislicealign: use exact SSP.
dchristiaens Jan 16, 2018
e324f9d
Move SVR to namespace; clean PSF code.
dchristiaens Jan 18, 2018
2d0d18f
dwimotioncorrect: switch everything to Lie representation.
dchristiaens Jan 18, 2018
dbb2399
maskslice: remove deprecated command
dchristiaens Jan 18, 2018
bb5874d
motionfiltter overhaul
dchristiaens Jan 18, 2018
bcedba2
dwislicealign: remove global scaling factor.
dchristiaens Jan 19, 2018
e08272e
dwimotioncorrect: enable motion filter only if slice order is explici…
dchristiaens Jan 20, 2018
dd83c37
dwirecon: disable minimum norm regularizer
dchristiaens Jan 20, 2018
deb4392
dwirecon: minimum norm regularizer at machine epsilon.
dchristiaens Jan 20, 2018
4c15a41
dwirecon: normalise regularization weight to no. volumes in the data.
dchristiaens Jan 22, 2018
b15b703
mrfieldunwarp: bug fix and improvement in s2v setting.
dchristiaens Jan 22, 2018
0430101
dwisliceoutliers: output 0 when mask is empty.
dchristiaens Jan 22, 2018
745a972
Merge branch 'master' into intensitymatching3
dchristiaens Feb 5, 2018
9771fe3
intensity matching: required operation order for motion filtering.
dchristiaens Feb 5, 2018
ec16911
intensity matching: normalise bias field in mask.
dchristiaens Feb 8, 2018
be0ef97
intensity matching: mask bias field.
dchristiaens Feb 8, 2018
a370b97
dwisliceoutliers: Correct loss function (derivative instead of ratio).
dchristiaens Feb 22, 2018
4ecc49a
dwimotioncorrect: avoid double interpolation in field mapping.
dchristiaens Feb 23, 2018
af44e4a
dwirecon: try new Jacobian approximation
dchristiaens Feb 26, 2018
5ebcfca
dwirecon: fix Jacobian modulation in field mapping.
dchristiaens Feb 28, 2018
75a6f1b
dwimotioncorrect: revert to field unwarping before recon.
dchristiaens Mar 1, 2018
9c59f91
slice-level intensity matching with normalisation
dchristiaens Mar 2, 2018
fdff015
Merge branch 'intensitymatching2a' into imscale2
dchristiaens Mar 19, 2018
305acc4
dwimotioncorrect: include local intensity matching.
dchristiaens Mar 19, 2018
d43c458
intensity matching: normalise in mask.
dchristiaens Mar 19, 2018
f9788ce
intensity matching: get bias in full image, normalise in amsk only.
dchristiaens Mar 19, 2018
e10415e
intensity matching: make full user option.
dchristiaens Mar 21, 2018
8092fd5
dwimotioncorrcet: increase no. inner iterations
dchristiaens Mar 23, 2018
83f8664
dwimotioncorrect: first setup for svr config list
dchristiaens Mar 23, 2018
c7fcd63
dwimotioncorrect: bug fix
dchristiaens Mar 23, 2018
343dfcf
dwirecon: switch to core library interpolators.
dchristiaens Apr 3, 2018
ba7faa8
dwirecon: clip at edges.
dchristiaens Apr 3, 2018
19c8987
dwirecon: set boundary conditions in Laplacian regularization.
dchristiaens Apr 3, 2018
383d94b
dwimotioncorrect: test single-shell recon.
dchristiaens Apr 4, 2018
17c8ac6
msshsvd: apply Laplace-Beltrami filter defore decomposition.
dchristiaens Apr 4, 2018
1d100d7
code cleanup
dchristiaens Apr 4, 2018
17e385c
Merge branch 'dev' into recon-core
dchristiaens Apr 4, 2018
981ce38
code cleanup
dchristiaens Apr 5, 2018
7f0f6f8
Merge branch 'recon-core' into dev
dchristiaens Apr 5, 2018
548ce88
dwimotioncorrect: multi-scale reconstruction using regularizer instea…
dchristiaens Apr 6, 2018
d4afe43
small bug fix
dchristiaens Apr 6, 2018
5a91ac9
dwisliceoutliers: force mask to be 3-D.
dchristiaens Apr 16, 2018
79ca22e
dwimotioncorrect: new default config with multi-scale reconstruction.
dchristiaens Apr 16, 2018
2a650f8
dwimotioncorrect: revert multi-scale registrtaion to separate filtering.
dchristiaens Apr 17, 2018
cd7e8e3
dwirecon: fix edge case for SSP width = 0.
dchristiaens Apr 17, 2018
fe4a30e
dwimotioncorrect: update for dev-RC3
dchristiaens May 2, 2018
c9d4653
dwirecon: change to LoG regularizer.
dchristiaens May 4, 2018
a159b48
dwislicealign: avoid unnecessary memory allocation in multi-threading.
dchristiaens May 4, 2018
d122138
dwisliceoutliers: let mask move alongside the image.
dchristiaens May 7, 2018
162a578
dwimotioncorrect: fix previous commit
dchristiaens May 7, 2018
0dfa5e2
dwirecon: switch to DoG regularizer
dchristiaens May 7, 2018
2a11cc6
dwirecon: code cleanup
dchristiaens May 8, 2018
cec1867
dwimotioncorrect: disable basis update in outer loop.
dchristiaens May 8, 2018
7052089
dwirecon: further code cleanup
dchristiaens May 8, 2018
c9b44dd
dwisliceoutliers: let mask move alongside the image (for real this ti…
dchristiaens May 9, 2018
5a6cd37
dwisliceoutliers: remove risk of race condition.
dchristiaens May 9, 2018
9d99648
dwislicealign: align mask with recon, using initial motion estimate.
dchristiaens May 9, 2018
210a7fb
dwirecon: switch to isotropic Laplacian regularizer at grid resolution.
dchristiaens May 11, 2018
a9b075e
dwimotioncorrect: fewer iterations in early stages
dchristiaens May 11, 2018
cb212fc
motionstats: new script for analysing motion trajectories and outliers.
dchristiaens May 13, 2018
6c405fa
motionstats: bug fixes and improvements
dchristiaens May 13, 2018
2e4f802
motionstats: measure intra-volume gradient dispersion.
dchristiaens May 13, 2018
676f9a8
dwisliceoutliers: use RMSE with scaling from voxel-wise MAD.
dchristiaens May 15, 2018
fc36416
dwisliceoutliers: fix multi-threading bug
dchristiaens May 16, 2018
f8db1d6
dwisliceoutliers: fix bug in dynamic masking.
dchristiaens May 16, 2018
e1e3169
dwimotioncorrect: adapt optimization scheme for Chi2 outlier reweight…
dchristiaens May 19, 2018
7f1dc56
dwirecon: incorporate slice-direction regularizer.
dchristiaens May 19, 2018
d60124e
Merge branch 'orchi2' into dev
dchristiaens May 19, 2018
8f54c0f
dwirecon: lower default regularization values
dchristiaens May 19, 2018
cfd46cf
dwirecon: bug fix and update motion corrcetion script
dchristiaens May 19, 2018
4cdce1d
dwimotioncorrect: take prior slice weights at input
dchristiaens Aug 9, 2018
b29fcdd
motionfilter: reinstate median filtering
dchristiaens Aug 9, 2018
82a94cf
dwisliceoutliers2: new bayesian slice outlier detection.
dchristiaens Aug 10, 2018
0d78288
dwisliceoutliers: avoid zero-division NaN.
dchristiaens Aug 10, 2018
a2514c3
dwisliceoutliers: trial new GMM-style OR.
dchristiaens Aug 28, 2018
7734cba
dwisliceoutliers: bug fix and initialisation
dchristiaens Aug 29, 2018
acd300c
dwimotioncorrect: switch outlier rejection to log-normal distributions.
dchristiaens Aug 30, 2018
8dffbac
Merge branch 'orprior' of gitlab.com:ChD/dwi-recon into orprior
dchristiaens Sep 4, 2018
94f033c
dwisliceoutliers2: fix EM implementation
dchristiaens Sep 4, 2018
dbbaff8
dwisliceoutliers2: refactor.
dchristiaens Sep 4, 2018
741e0d3
dwisliceoutliergmm: rename command
dchristiaens Sep 4, 2018
6410a59
dwisliceoutliergmm: switch to log-normal distributions.
dchristiaens Sep 4, 2018
a6118df
dwimotioncorrect: integrate new dwisliceoutliergmm command.
dchristiaens Sep 4, 2018
149ba3f
dwisliceoutliergmm: fix bugs.
dchristiaens Sep 6, 2018
343a0b7
sliceoutliergmm: avoid NaN output due to numeric overflow.
dchristiaens Sep 18, 2018
bc48474
dwimotioncorrect: add support for fixed slice weights.
dchristiaens Oct 30, 2018
b391c62
dwirecon: revert Laplacial regularizer to minimal convolution filter.
dchristiaens Nov 16, 2018
0eebd00
mrfieldunwarp: add option to disable Jacobian modulation.
dchristiaens Nov 16, 2018
18a654b
dwirecon: code cleanup
dchristiaens Nov 16, 2018
ff99b72
dwisliceoutliergmm: fix GMM initialisation for log-normal data.
dchristiaens Dec 6, 2018
b871aef
dwirecon: normalise slice outlier weights in reconstruction.
dchristiaens Dec 13, 2018
0375f5c
dwirecon: scale differential operators in regularization.
dchristiaens Dec 19, 2018
23cc3eb
dwimotioncorrect: add flag for volume-level correction
dchristiaens Apr 9, 2019
ed5fc5a
Update README.md
dchristiaens May 3, 2019
0705aa3
Merge branch 'master' of gitlab.com:ChD/dwi-recon
dchristiaens May 3, 2019
52b185e
Remove deprecated commands
dchristiaens May 3, 2019
a8af209
Update copyright notice.
dchristiaens May 3, 2019
d64a289
Fix typo in copyright notice.
dchristiaens May 29, 2019
2849313
msshsvd: simplify and speed-up
dchristiaens Aug 20, 2019
c016b16
dwirecon: add capability for voxel-level outlier weights
dchristiaens Sep 23, 2019
e19d2d3
dwimotioncorrect: unwarp voxel weights to ensure spatial correspondence.
dchristiaens Sep 23, 2019
535f9f9
mrfieldunwarp: linear interpolation
dchristiaens Sep 23, 2019
bb8d2a4
dwirecon: fix bug with local outlier weight handling.
dchristiaens Sep 24, 2019
062c74c
update gitignore
dchristiaens Oct 2, 2019
55f60d5
dwirecon: switch to LSCG solver
dchristiaens Oct 2, 2019
579fda1
dwirecon: use regularisers in LSCG
dchristiaens Oct 3, 2019
180f7a0
dwirecon: balance weights symmetrically in LSCG solver
dchristiaens Oct 3, 2019
b5f127a
dwirecon: reinstate weight normalisation
dchristiaens Oct 3, 2019
0fb8f40
dwirecon: code cleanup
dchristiaens Oct 3, 2019
cdeef45
dwislicealign: fix assert error
dchristiaens Oct 4, 2019
f209890
Merge branch 'master' into recon-refactoring
dchristiaens Oct 4, 2019
1b74cf8
dwirecon: improve run time efficiency
dchristiaens Oct 11, 2019
2bfb4ad
dwisliceoutliergmm: fix floating point rounding error
dchristiaens Oct 11, 2019
8d4372d
Merge branch 'master' into recon-refactoring
dchristiaens Oct 11, 2019
f4a2001
dwirecon: fix byte alignment issues
dchristiaens Oct 11, 2019
f95ee97
dwirecon: encapsulate q-space basis in seperate class.
dchristiaens Oct 15, 2019
1857204
dwirecon: start work on adapter pipes.
dchristiaens Oct 16, 2019
2e42977
dwirecon: fix bugs introduced in the refactoring
dchristiaens Oct 17, 2019
9281398
dwirecon: run time improvement
dchristiaens Oct 17, 2019
49e71c8
dwirecon: continue refactoring
dchristiaens Oct 17, 2019
4477273
dwirecon refactor: implement y->x mapping
dchristiaens Oct 18, 2019
153e54d
dwirecon refactoring: first working version
dchristiaens Oct 18, 2019
1ead869
dwirecon refactor: ensure thread safety
dchristiaens Oct 21, 2019
f604b0e
dwirecon refactor: fix bytesize allocation for atomic_flag
dchristiaens Oct 21, 2019
684097f
dwirecon refactor: rename Cache adapter and use memset for speed.
dchristiaens Oct 22, 2019
4e7c0f0
Merge branch 'recon-refactoring' of https://gitlab.com/ChD/dwi-recon …
dchristiaens Oct 23, 2019
83fcef7
dwirecon refactor: fix bug in memset call
dchristiaens Oct 23, 2019
6908626
dwirecon refactor: reinstate SSP handling
dchristiaens Oct 24, 2019
8e0dc41
dwirecon refactor: run time optimization
dchristiaens Oct 24, 2019
1bebda3
dwirecon refactor: improve (double) cache efficiency
dchristiaens Oct 30, 2019
90505c9
dwirecon refactor: add recon_proj command for testing
dchristiaens Oct 30, 2019
b085569
Merge branch 'recon-refactoring' into voxelweights
dchristiaens Nov 5, 2019
e5b37e8
dwirecon: enable reconstruction on arbitrary grid.
dchristiaens Nov 25, 2019
a01e31a
dwirecon: fix issue with negative strides.
dchristiaens Nov 27, 2019
da4daba
dwirecon: fix issue with multi-resolution grid header
dchristiaens Dec 6, 2019
7238592
Merge branch 'recon-refactoring' into voxelweights
dchristiaens Dec 10, 2019
416b72a
dwirecon: clean up commented legacy code
dchristiaens Mar 3, 2020
63fc542
dwimotioncorrect: refactor to new scripting syntax
dchristiaens Mar 4, 2020
db2a4c6
dwimotioncorrect: fix scripting bug due to use of 'global'
dchristiaens Mar 5, 2020
feaab19
Merge branch 'recon-refactoring' into voxelweights
dchristiaens Mar 25, 2020
3c38ebe
dwimotioncorrect: fix bug in voxel weights file name handling
dchristiaens Apr 1, 2020
6d890eb
dwimotioncorrect: add option to fix motion traces
dchristiaens Apr 27, 2020
9339a42
dwimotioncorrect: make motion filter part of default setup
dchristiaens Apr 27, 2020
ffc98d3
Merge branch 'recon-refactoring' of gitlab.com:ChD/dwi-recon into rec…
dchristiaens Apr 28, 2020
dd1de50
Update codebase for MRtrix 3.0.0 interface.
dchristiaens May 4, 2020
0fba182
dwimotioncorrect: temporarily increase no. iterations in default config
dchristiaens May 4, 2020
d279c5e
dwimotioncorrect: fix typo
dchristiaens May 4, 2020
0a82b1d
msshsvd: fix segfault when no mask is provided
dchristiaens May 4, 2020
c43ab5b
Revert "dwimotioncorrect: temporarily increase no. iterations in defa…
dchristiaens May 5, 2020
04fae54
dwimotioncorrect: set up config json file
dchristiaens May 18, 2020
637f760
dwimotioncorrect: fix stupid and less stupid bugs
dchristiaens May 18, 2020
9c1e126
dwirecon: allow dynamic SSP length
dchristiaens May 18, 2020
9b91ec8
cosmetic changes and error checking
dchristiaens May 19, 2020
6124301
dwimotioncorrect: override config with custom -reg and -zreg params
dchristiaens May 21, 2020
13f18a9
dwirecon: compatibility with mrtrix 3.0.1
dchristiaens Jul 6, 2020
3c95095
Update README.md
dchristiaens Oct 23, 2020
27cac3f
dwimotioncorrect: error updated to mrtrix api
maxpietsch Nov 2, 2020
d51a898
README: fix install instructions
dchristiaens Nov 3, 2020
576808a
Merge branch 'dwimotioncorrect_errorfix' into 'master'
dchristiaens Nov 9, 2020
bae4595
dwimotioncorrect: fix bug in config file handling
dchristiaens Dec 14, 2020
2d5076e
mrfieldunwarp: warning when dimensions don't match
dchristiaens Jan 14, 2021
9e1946b
mrfieldunwarp: improve docs
dchristiaens Jan 14, 2021
71ab3c5
Merge pull request #5 from dchristiaens/fieldmap_checks
dchristiaens Jan 14, 2021
0ce8a11
dwimotioncorrect: save command history to output
dchristiaens Jan 14, 2021
821ab7b
README.md updates
dchristiaens Jan 14, 2021
7cba1cd
dwirecon: fix segfault reported in #4
dchristiaens Jan 15, 2021
dc01d94
dwimotioncorrect: check scipy availability upfront
maxpietsch Jan 27, 2021
f878133
Merge pull request #7 from dchristiaens/save_cmdhistory
maxpietsch Jan 27, 2021
82b2e6c
Merge pull request #8 from dchristiaens/fix_segfault
maxpietsch Jan 27, 2021
55da420
dwimotioncorrect: use mssh as reference for command history
maxpietsch Jan 28, 2021
162f921
added drift motion filter option
maxpietsch Jan 28, 2021
ab5c742
Merge pull request #12 from dchristiaens/save_cmdhistory
dchristiaens Feb 18, 2021
5735923
Merge pull request #10 from dchristiaens/import_check
dchristiaens Feb 18, 2021
af087e4
dwimotioncorrect: default LB regularization 0.001
dchristiaens Jun 25, 2021
4d58d97
add support for Eigen 3.3.7 and above
dchristiaens Oct 29, 2021
219c93a
dwimotioncorrect: fix bug on NaN voxel size in 4th dimension
dchristiaens Nov 2, 2021
9bf9ba2
mssh2amp: fix bug on [ 0 0 0 ] input directions.
dchristiaens Nov 2, 2021
2beb61f
Merge branch 'master' of gitlab.com:ChD/dwi-recon
dchristiaens Nov 2, 2021
df699d7
dwirecon: support Eigen 3.4.0
dchristiaens Nov 17, 2021
d91962e
dwirecon: remove deprecated Loop({4,3}) syntax #8
dchristiaens Feb 21, 2022
d06c72b
Merge pull request #16 from dchristiaens/fix_loop
dchristiaens Feb 21, 2022
0cbcd9d
fix nthreads
treanus Feb 21, 2022
3dd671a
Merge pull request #17 from treanus/master
dchristiaens Feb 21, 2022
30ac0b1
dwimotioncorrect: smarter default for -rlmax
dchristiaens Apr 11, 2022
79afb72
dwimotioncorrect: smarter default for -lmax
dchristiaens Apr 11, 2022
dfcfaf1
Merge branch 'master' into drift_filter
dchristiaens Apr 11, 2022
239a03d
Merge pull request #13 from dchristiaens/drift_filter
dchristiaens Apr 11, 2022
0269bb9
dwimotioncorrect: bug fix
dchristiaens Apr 13, 2022
ca24e73
Merge pull request #19 from dchristiaens/rlmax_default
dchristiaens Apr 13, 2022
49a6efa
Merge remote-tracking branch 'shard/master' into merge_shard
dchristiaens Sep 11, 2024
8e4bd40
Remove command dwirecon_proj
dchristiaens Sep 12, 2024
f1f9770
Update C++ code for recent changes on dev.
dchristiaens Sep 12, 2024
3d8f36f
Update Python code for recent changes on dev.
dchristiaens Sep 12, 2024
b80bb5d
Fix runtime errors.
dchristiaens Sep 13, 2024
2cf67ed
motionstats: output mean rotation in degrees
dchristiaens Sep 13, 2024
72b736a
shard: use constexpressions
dchristiaens Sep 13, 2024
3659512
shard: use pragma once
dchristiaens Sep 13, 2024
6f45946
shard: use nested namespaces
dchristiaens Sep 13, 2024
6f011e4
dwirecon: remove deprecated -field option
dchristiaens Sep 13, 2024
ee8ffb9
shard: apply formatting rules
dchristiaens Sep 13, 2024
e7e0f61
Apply first batch of clang-tidy suggestions
dchristiaens Sep 18, 2024
cc51416
Apply second batch of clang-tidy suggestions
dchristiaens Sep 18, 2024
2172535
Replace C-style casts with C++ casts in dwirecon
daljit46 Sep 26, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
372 changes: 372 additions & 0 deletions cmd/dwirecon.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,372 @@
/* Copyright (c) 2017-2019 Daan Christiaens
*
* MRtrix and this add-on module are distributed in the hope
* that it will be useful, but WITHOUT ANY WARRANTY; without
* even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE.
*/

#include <algorithm>
#include <sstream>

#include "adapter/extract.h"
#include "command.h"
#include "dwi/gradient.h"
#include "dwi/shells.h"
#include "file/matrix.h"
#include "image.h"
#include "math/SH.h"
#include "phase_encoding.h"

#include "dwi/svr/qspacebasis.h"
#include "dwi/svr/recon.h"

constexpr int DEFAULT_LMAX = 4;
constexpr float DEFAULT_SSPW = 1.0F;
constexpr float DEFAULT_REG = 1e-3F;
constexpr float DEFAULT_ZREG = 1e-3F;
constexpr float DEFAULT_TOL = 1e-4F;
constexpr int DEFAULT_MAXITER = 10;

using namespace MR;
using namespace App;

// clang-format off
void usage ()
{
AUTHOR = "Daan Christiaens ([email protected])";

SYNOPSIS = "Reconstruct DWI signal from a series of scattered slices with associated motion parameters.";

DESCRIPTION
+ "";

ARGUMENTS
+ Argument ("DWI", "the input DWI image.").type_image_in()
+ Argument ("SH", "the output spherical harmonics coefficients image.").type_image_out();


OPTIONS
+ Option ("motion", "The motion parameters associated with input slices or volumes. "
"These are supplied as a matrix of 6 columns encoding the rigid "
"transformations w.r.t. scanner space in se(3) Lie algebra." )
+ Argument ("file").type_file_in()

+ Option ("rf", "Basis functions for the radial (multi-shell) domain, provided as matrices in which "
"rows correspond with shells and columns with SH harmonic bands.").allow_multiple()
+ Argument ("b").type_file_in()

+ Option ("lmax", "The maximum harmonic order for the output series. (default = " + str(DEFAULT_LMAX) + ")")
+ Argument ("order").type_integer(0, 30)

+ Option ("weights", "Slice weights, provided as a matrix of dimensions Nslices x Nvols.")
+ Argument ("W").type_file_in()

+ Option ("voxweights", "Voxel weights, provided as an image of same dimensions as dMRI data.")
+ Argument ("W").type_image_in()

+ Option ("ssp", "Slice sensitivity profile, either as text file or as a scalar slice thickness for a "
"Gaussian SSP, relative to the voxel size. (default = " + str(DEFAULT_SSPW) + ")")
+ Argument ("w").type_text()

+ Option ("reg", "Isotropic Laplacian regularization. (default = " + str(DEFAULT_REG) + ")")
+ Argument ("l").type_float()

+ Option ("zreg", "Regularization in the slice direction. (default = " + str(DEFAULT_ZREG) + ")")
+ Argument ("l").type_float()

+ Option ("template", "Template header to determine the reconstruction grid.")
+ Argument ("header").type_image_in()

+ DWI::GradImportOptions()

+ PhaseEncoding::ImportOptions

+ DWI::ShellsOption

+ OptionGroup ("Output options")

+ Option ("spred",
"output source prediction of all scattered slices. (useful for diagnostics)")
+ Argument ("out").type_image_out()

+ Option ("padding", "zero-padding output coefficients to given dimension.")
+ Argument ("rank").type_integer(0)

+ Option ("complete", "complete (zero-filled) source prediction.")

+ OptionGroup ("CG Optimization options")

+ Option ("tolerance", "the tolerance on the conjugate gradient solver. (default = " + str(DEFAULT_TOL) + ")")
+ Argument ("t").type_float(0.0, 1.0)

+ Option ("maxiter",
"the maximum number of iterations of the conjugate gradient solver. (default = " + str(DEFAULT_MAXITER) + ")")
+ Argument ("n").type_integer(1)

+ Option ("init",
"initial guess of the reconstruction parameters.")
+ Argument ("img").type_image_in();

}
// clang-format on

typedef float value_type;

void run() {
dchristiaens marked this conversation as resolved.
Show resolved Hide resolved
dchristiaens marked this conversation as resolved.
Show resolved Hide resolved
// Load input data
auto dwi = Image<value_type>::open(argument[0]).with_direct_io({1, 2, 3, 4});

// Read motion parameters
auto opt = get_options("motion");
Eigen::MatrixXf motion;
if (!opt.empty())
motion = File::Matrix::load_matrix<float>(opt[0][0]);
else
motion = Eigen::MatrixXf::Zero(dwi.size(3), 6);

// Check dimensions
if ((motion.size() != 0) && (motion.cols() != 6))
throw Exception("No. columns in motion parameters must equal 6.");
if ((motion.size() != 0) && (((dwi.size(3) * dwi.size(2)) % motion.rows()) != 0))
throw Exception("No. rows in motion parameters does not match image dimensions.");

// Select shells
auto grad = DWI::get_DW_scheme(dwi);
DWI::Shells shells(grad);
shells.select_shells(false, false, false);

// Read multi-shell basis
int lmax = 0;
std::vector<Eigen::MatrixXf> rf;
opt = get_options("rf");
for (size_t k = 0; k < opt.size(); k++) {
Eigen::MatrixXf const t = File::Matrix::load_matrix<float>(opt[k][0]);
if (t.rows() != shells.count())
throw Exception("No. shells does not match no. rows in basis function " + opt[k][0] + ".");
lmax = std::max(2 * (int(t.cols()) - 1), lmax);
rf.push_back(t);
}

// Read slice weights
Eigen::MatrixXf W = Eigen::MatrixXf::Ones(dwi.size(2), dwi.size(3));
opt = get_options("weights");
if (!opt.empty()) {
W = File::Matrix::load_matrix<float>(opt[0][0]);
if (W.rows() != dwi.size(2) || W.cols() != dwi.size(3))
throw Exception("Weights matrix dimensions don't match image dimensions.");
}

// Get volume indices
std::vector<size_t> idx;
if (rf.empty()) {
idx = shells.largest().get_volumes();
} else {
for (size_t k = 0; k < shells.count(); k++)
idx.insert(idx.end(), shells[k].get_volumes().begin(), shells[k].get_volumes().end());
std::sort(idx.begin(), idx.end());
}

// Select subset
auto dwisub = Adapter::make<Adapter::Extract1D>(dwi, 3, container_cast<std::vector<uint32_t>>(idx));

Eigen::MatrixXf gradsub(idx.size(), grad.cols());
for (ssize_t i = 0; i < idx.size(); i++)
gradsub.row(i) = grad.row(ssize_t(idx[i])).template cast<float>();

ssize_t const ne = motion.rows() / dwi.size(3);
Eigen::MatrixXf motionsub(ne * idx.size(), 6);
for (ssize_t i = 0; i < idx.size(); i++)
for (ssize_t j = 0; j < ne; j++)
motionsub.row(i * ne + j) = motion.row(ssize_t(idx[i]) * ne + j);

Eigen::MatrixXf Wsub(W.rows(), idx.size());
for (ssize_t i = 0; i < idx.size(); i++)
Wsub.col(i) = W.col(ssize_t(idx[i]));

// SSP
DWI::SVR::SSP<float> ssp(DEFAULT_SSPW);
opt = get_options("ssp");
if (!opt.empty()) {
std::string const t = opt[0][0];
try {
ssp = DWI::SVR::SSP<float>(std::stof(t));
} catch (std::invalid_argument &e) {
try {
Eigen::VectorXf const v = File::Matrix::load_vector<float>(t);
ssp = DWI::SVR::SSP<float>(v);
} catch (...) {
throw Exception("Invalid argument for SSP.");
}
}
}

// Read voxel weights
Eigen::VectorXf Wvox = Eigen::VectorXf::Ones(dwisub.size(0) * dwisub.size(1) * dwisub.size(2) * dwisub.size(3));
opt = get_options("voxweights");
if (!opt.empty()) {
auto voxweights = Image<value_type>::open(opt[0][0]);
check_dimensions(dwisub, voxweights, 0, 4);
ssize_t j = 0;
for (auto l = Loop("loading voxel weights data", {0, 1, 2, 3})(voxweights); l; l++, j++) {
Wvox[j] = voxweights.value();
dchristiaens marked this conversation as resolved.
Show resolved Hide resolved
}
}

// Other parameters
if (rf.empty())
lmax = get_option_value("lmax", DEFAULT_LMAX);
else
lmax = std::min(lmax, (int)get_option_value("lmax", lmax));

float const reg = get_option_value("reg", DEFAULT_REG);
float const zreg = get_option_value("zreg", DEFAULT_ZREG);

value_type const tol = get_option_value("tolerance", DEFAULT_TOL);
ssize_t const maxiter = get_option_value("maxiter", DEFAULT_MAXITER);

DWI::SVR::QSpaceBasis const qbasis{gradsub, lmax, rf, motionsub};

ssize_t const ncoefs = ssize_t(qbasis.get_ncoefs());
size_t const padding = get_option_value("padding", Math::SH::NforL(lmax));
if (padding < Math::SH::NforL(lmax))
throw Exception("user-provided padding too small.");

// Create source header - needed due to stride handling
Header srchdr(dwisub);
Stride::set(srchdr, {1, 2, 3, 4});
DWI::set_DW_scheme(srchdr, gradsub);
srchdr.datatype() = DataType::Float32;
srchdr.sanitise();

// Create recon header
Header rechdr(dwisub);
opt = get_options("template");
if (!opt.empty()) {
rechdr = Header::open(opt[0][0]);
}
rechdr.ndim() = 4;
rechdr.size(3) = ncoefs;
dchristiaens marked this conversation as resolved.
Show resolved Hide resolved
Stride::set(rechdr, {2, 3, 4, 1});
rechdr.datatype() = DataType::Float32;
rechdr.sanitise();

// Create mapping
DWI::SVR::ReconMapping const map(rechdr, srchdr, qbasis, motionsub, ssp);

// Set up scattered data matrix
INFO("initialise reconstruction matrix");
DWI::SVR::ReconMatrix R(map, reg, zreg);
R.setWeights(Wsub);

R.setVoxelWeights(Wvox);

// Read input data to vector (this enforces positive strides!)
Eigen::VectorXf y(R.rows());
y.setZero();
ssize_t j = 0;
for (auto lv = Loop("loading image data", {0, 1, 2, 3})(dwisub); lv; lv++, j++) {
float const w = Wsub(dwisub.index(2), dwisub.index(3)) * Wvox[j];

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no viable conversion from 'const CwiseBinaryOp<internal::scalar_product_op<typename internal::traits<IndexedView<Matrix<float, -1, -1, 0, -1, -1>, Index<Extract1D<Image>>, Index<Extract1D<Image>>>>::Scalar, typename internal::promote_scalar_arg<Scalar, float, (Eigen::internal::has_ReturnType<Eigen::ScalarBinaryOpTraits<Scalar, float, Eigen::internal::scalar_product_op<Scalar, float>>>::value)>::type>, const IndexedView<Matrix<float, -1, -1, 0, -1, -1>, Index<Extract1D<Image>>, Index<Extract1D<Image>>>, const typename internal::plain_constant_type<IndexedView<Matrix<float, -1, -1, 0, -1, -1>, Index<Extract1D<Image>>, Index<Extract1D<Image>>>, typename internal::promote_scalar_arg<Scalar, float, (Eigen::internal::has_ReturnType<Eigen::ScalarBinaryOpTraits<Scalar, float, Eigen::internal::scalar_product_op<Scalar, float>>>::value)>::type>::type>' (aka 'const CwiseBinaryOp<scalar_product_op<float, float>, const Eigen::IndexedView<Eigen::Matrix<float, -1, -1, 0, -1, -1>, MR::Helper::Index<MR::Adapter::Extract1D<MR::Image>>, MR::Helper::Index<MR::Adapter::Extract1D<MR::Image>>>, const CwiseNullaryOp<scalar_constant_op, const Eigen::Matrix<float, -1, -1, 0, -1, -1>>>') to 'const float' [clang-diagnostic-error]

    float const w = Wsub(dwisub.index(2), dwisub.index(3)) * Wvox[j];
                ^

y[j] = std::sqrt(w) * dwisub.value();
dchristiaens marked this conversation as resolved.
Show resolved Hide resolved
}

// Fit scattered data in basis...
INFO("initialise conjugate gradient solver");

Eigen::LeastSquaresConjugateGradient<DWI::SVR::ReconMatrix, Eigen::IdentityPreconditioner> cg;
cg.compute(R);

cg.setTolerance(tol);
cg.setMaxIterations(maxiter);
dchristiaens marked this conversation as resolved.
Show resolved Hide resolved

// Solve y = M x
Eigen::VectorXf x(R.cols());
opt = get_options("init");
if (!opt.empty()) {
// load initialisation
auto init = Image<value_type>::open(opt[0][0]).with_direct_io({3, 4, 5, 2, 1});
check_dimensions(rechdr, init, 0, 3);
if ((init.size(3) != shells.count()) || (init.size(4) < Math::SH::NforL(lmax)))
throw Exception("dimensions of init image don't match.");
// init vector
Eigen::VectorXf x0(R.cols());
// convert from mssh
Eigen::VectorXf c(shells.count() * Math::SH::NforL(lmax));
Eigen::MatrixXf x2mssh(c.size(), ncoefs);
x2mssh.setZero();
for (int k = 0; k < shells.count(); k++)
x2mssh.middleRows(ssize_t(k * Math::SH::NforL(lmax)), ssize_t(Math::SH::NforL(lmax))) = qbasis.getShellBasis(k).transpose();
auto mssh2x = x2mssh.fullPivHouseholderQr();
ssize_t j = 0;
ssize_t k = 0;
for (auto l = Loop("loading initialisation", {0, 1, 2})(init); l; l++, j += ncoefs) {
k = 0;
for (auto l2 = Loop(3)(init); l2; l2++) {
for (init.index(4) = 0; init.index(4) < Math::SH::NforL(lmax); init.index(4)++)
c[k++] = std::isfinite((float)init.value()) ? init.value() : 0.0F;
}
x0.segment(j, ncoefs) = mssh2x.solve(c);
dchristiaens marked this conversation as resolved.
Show resolved Hide resolved
}
INFO("solve from given starting point");
x = cg.solveWithGuess(y, x0);
} else {
INFO("solve from zero starting point");
x = cg.solve(y);
}

CONSOLE("CG: #iterations: " + str(cg.iterations()));
CONSOLE("CG: estimated error: " + str(cg.error()));

// Write result to output file
Header msshhdr(rechdr);
msshhdr.ndim() = 5;
msshhdr.size(3) = ssize_t(shells.count());
msshhdr.size(4) = ssize_t(padding);
Stride::set_from_command_line(msshhdr, {3, 4, 5, 2, 1});
msshhdr.datatype() = DataType::from_command_line(DataType::Float32);
PhaseEncoding::set_scheme(msshhdr, Eigen::MatrixXf());
// store b-values and counts
{
std::stringstream ss;
for (auto b : shells.get_bvalues())
ss << b << ",";
std::string const key = "shells";
std::string val = ss.str();
val.erase(val.length() - 1);
msshhdr.keyval()[key] = val;
}
{
std::stringstream ss;
for (auto c : shells.get_counts())
ss << c << ",";
std::string const key = "shellcounts";
std::string val = ss.str();
val.erase(val.length() - 1);
msshhdr.keyval()[key] = val;
}

auto out = Image<value_type>::create(argument[1], msshhdr);

j = 0;
Eigen::VectorXf c(ncoefs);
Eigen::VectorXf sh(padding);
sh.setZero();
for (auto l = Loop("writing result to image", {0, 1, 2})(out); l; l++, j += ncoefs) {
c = x.segment(j, ncoefs);
dchristiaens marked this conversation as resolved.
Show resolved Hide resolved
for (int k = 0; k < shells.count(); k++) {
out.index(3) = k;
sh.head(Math::SH::NforL(lmax)) = qbasis.getShellBasis(k).transpose() * c;
out.row(4) = sh;
}
}

// Output source prediction
bool const complete = !get_options("complete").empty();
opt = get_options("spred");
if (!opt.empty()) {
srchdr.size(3) = (complete) ? dwi.size(3) : dwisub.size(3);
auto spred = Image<value_type>::create(opt[0][0], srchdr);
auto recon = ImageView<value_type>(rechdr, x.data());
map.x2y(recon, spred);
}
}
Loading
Loading