diff --git a/.circleci/config.yml b/.circleci/config.yml index bf445fd8d..670f6e1e3 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -339,6 +339,8 @@ jobs: root: /src/coverage paths: - .coverage.examples_pasl_multipld + # The resource_class feature allows configuring CPU and RAM resources for each job. Different resource classes are available for different executors. https://circleci.com/docs/2.0/configuration-reference/#resourceclass + resource_class: large qtab: <<: *dockersetup diff --git a/LICENSE.md b/LICENSE.md index c6a850dc5..f68c25e2d 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -25,30 +25,3 @@ SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - - -## External code - -### aslprep.utils.misc.estimate_att_pcasl - -MIT License - -Copyright (c) 2023 willtack - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. diff --git a/aslprep/data/aslprep_bids_config.json b/aslprep/data/aslprep_bids_config.json index 755620537..f0c88121b 100755 --- a/aslprep/data/aslprep_bids_config.json +++ b/aslprep/data/aslprep_bids_config.json @@ -37,9 +37,9 @@ "sub-{subject}[/ses-{session}]/[{datatype|func}/]sub-{subject}[_ses-{session}]_task-{task}[_acq-{acquisition}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_recording-{recording}]_{suffix}{extension<.tsv.gz|.json>|.tsv.gz}", "sub-{subject}[/ses-{session}]/{datatype|func}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_space-{space}][_atlas-{atlas}][_cohort-{cohort}][_desc-{desc}]_{suffix|timeseries}{extension<.json|.tsv|.csv|>|.tsv}", "sub-{subject}[/ses-{session}]/{datatype|perf}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_desc-{desc}]_{suffix}{extension<.tsv|.json>|.tsv}", - "sub-{subject}[/ses-{session}]/{datatype|perf}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_space-{space}][_res-{res}][_atlas-{atlas}][_cohort-{cohort}][_desc-{desc}]_{suffix}{extension<.nii|.nii.gz|.json|.csv|.tsv>|.nii.gz}", - "sub-{subject}[/ses-{session}]/{datatype|perf}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_space-{space}][_res-{res}][_den-{den}][_cohort-{cohort}][_desc-{desc}]_{suffix}{extension<.dtseries.nii|.dscalar.nii|.json|.dscalar.json|.dtseries.json|.func.gii|.func.json>|.dtseries.nii}", - "sub-{subject}[/ses-{session}]/{datatype|figures}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_space-{space}][_res-{res}][_den-{den}][_cohort-{cohort}][_desc-{desc}]_{suffix}{extension<.svg|.png|.html>|.svg}", + "sub-{subject}[/ses-{session}]/{datatype|perf}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_space-{space}][_res-{res}][_atlas-{atlas}][_cohort-{cohort}][_desc-{desc}]_{suffix}{extension<.nii|.nii.gz|.json|.csv|.tsv>|.nii.gz}", + "sub-{subject}[/ses-{session}]/{datatype|perf}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_space-{space}][_res-{res}][_den-{den}][_cohort-{cohort}][_desc-{desc}]_{suffix}{extension<.dtseries.nii|.dscalar.nii|.json|.dscalar.json|.dtseries.json|.func.gii|.func.json>|.dtseries.nii}", + "sub-{subject}[/ses-{session}]/{datatype|figures}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_space-{space}][_res-{res}][_den-{den}][_cohort-{cohort}][_desc-{desc}]_{suffix}{extension<.svg|.png|.html>|.svg}", "sub-{subject}/{datatype}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_run-{run}][_space-{space}][_atlas-{atlas}][_cohort-{cohort}][_desc-{desc}]_{suffix}{extension<.html|.svg|.png>}", "sub-{subject}/{datatype}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_run-{run}][_space-{space}][_atlas-{atlas}][_cohort-{cohort}][_desc-{desc}]_{suffix}{extension<.html|.svg|.png>}", "sub-{subject}[/ses-{session}]/{datatype|anat}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_run-{run}]_from-{from}_to-{to}_mode-{mode|image}_{suffix|xfm}{extension<.txt|.h5|.json>}", diff --git a/aslprep/data/boilerplate.bib b/aslprep/data/boilerplate.bib index 1acdf6503..a162a1d00 100644 --- a/aslprep/data/boilerplate.bib +++ b/aslprep/data/boilerplate.bib @@ -557,27 +557,6 @@ @article{groves2009combined publisher={Elsevier} } -@article{dai2012reduced, - title={Reduced resolution transit delay prescan for quantitative continuous arterial spin labeling perfusion imaging}, - author={Dai, Weiying and Robson, Philip M and Shankaranarayanan, Ajit and Alsop, David C}, - journal={Magnetic resonance in medicine}, - volume={67}, - number={5}, - pages={1252--1265}, - year={2012}, - publisher={Wiley Online Library} -} - -@article{wang2013multi, - title={Multi-delay multi-parametric arterial spin-labeled perfusion MRI in acute ischemic stroke—comparison with dynamic susceptibility contrast enhanced perfusion imaging}, - author={Wang, Danny JJ and Alger, Jeffry R and Qiao, Joe X and Gunther, Matthias and Pope, Whitney B and Saver, Jeffrey L and Salamon, Noriko and Liebeskind, David S and UCLA Stroke Investigators and others}, - journal={NeuroImage: Clinical}, - volume={3}, - pages={1--7}, - year={2013}, - publisher={Elsevier} -} - @article{clement2022asl, title={ASL-BIDS, the brain imaging data structure extension for arterial spin labeling}, author={Clement, Patricia and Castellaro, Marco and Okell, Thomas W and Thomas, David L and Vandemaele, Pieter and Elgayar, Sara and Oliver-Taylor, Aaron and Kirk, Thomas and Woods, Joseph G and Vos, Sjoerd B and others}, @@ -615,31 +594,6 @@ @article{wang2005amplitude doi={10.1148/radiol.2351031663} } -@article{fan2017long, - title={Long-delay arterial spin labeling provides more accurate cerebral blood flow measurements in moyamoya patients: a simultaneous positron emission tomography/MRI study}, - author={Fan, Audrey P and Guo, Jia and Khalighi, Mohammad M and Gulaka, Praveen K and Shen, Bin and Park, Jun Hyung and Gandhi, Harsh and Holley, Dawn and Rutledge, Omar and Singh, Prachi and others}, - journal={Stroke}, - volume={48}, - number={9}, - pages={2441--2449}, - year={2017}, - publisher={Am Heart Assoc}, - url={https://doi.org/10.1161/STROKEAHA.117.017773}, - doi={10.1161/STROKEAHA.117.017773} -} - -@article{juttukonda2021characterizing, - title={Characterizing cerebral hemodynamics across the adult lifespan with arterial spin labeling MRI data from the Human Connectome Project-Aging}, - author={Juttukonda, Meher R and Li, Binyin and Almaktoum, Randa and Stephens, Kimberly A and Yochim, Kathryn M and Yacoub, Essa and Buckner, Randy L and Salat, David H}, - journal={Neuroimage}, - volume={230}, - pages={117807}, - year={2021}, - publisher={Elsevier}, - url={https://doi.org/10.1016/j.neuroimage.2021.117807}, - doi={10.1016/j.neuroimage.2021.117807} -} - @article{noguchi2015technical, title={A technical perspective for understanding quantitative arterial spin-labeling MR imaging using Q2TIPS}, author={Noguchi, Tomoyuki and Nishihara, Masashi and Hara, Yukiko and Hirai, Tetsuyoshi and Egashira, Yoshiaki and Azama, Shinya and Irie, Hiroyuki}, @@ -852,3 +806,13 @@ @article{topup doi = {10.1016/S1053-8119(03)00336-7}, url = {https://www.sciencedirect.com/science/article/pii/S1053811903003367} } + +@article{woods2023recommendations, + title={Recommendations for quantitative cerebral perfusion MRI using multi-timepoint arterial spin labeling: Acquisition, quantification, and clinical applications}, + author={Woods, Joseph G and Achten, Eric and Asllani, Iris and Bolar, Divya S and Dai, Weiying and Detre, John A and Fan, Audrey P and Fern{\'a}ndez-Seara, Maria and Golay, Xavier and G{\"u}nther, Matthias and Guo, Jia and Hernandez-Garcia, Luis and Ho, Mai-Lan and Juttukonda, Meher R. and Lu, Hanzhang and MacIntosh, Bradley J. and Madhuranthakam, Ananth J. and Mutsaerts, Henk J.M.M. and Okell, Thomas W. and Parkes, Laura M. and Pinter, Nandor and Pinto, Joana and Qin, Qin and Smits, Marion and Suzuki, Yuriko and Thomas, David L. and van Osch, Matthias J.P. and Wang, Danny J.J. and Warnert, Esther A.H. and Zaharchuk, Greg and Zelaya, Fernando and Zhao, Moss and Chappell, Michael A.}, + year={2023}, + journal={OSF Preprints}, + publisher={OSF Preprints}, + doi={10.31219/osf.io/4tskr}, + url={https://doi.org/10.31219/osf.io/4tskr}, +} diff --git a/aslprep/data/reports-spec.yml b/aslprep/data/reports-spec.yml index a2491e873..7c4ed2e87 100644 --- a/aslprep/data/reports-spec.yml +++ b/aslprep/data/reports-spec.yml @@ -11,30 +11,31 @@ sections: desc: conform extension: [.html] suffix: T1w - - bids: {datatype: figures, suffix: dseg} + - subtitle: Brain mask and brain tissue segmentation of the T1w + bids: {datatype: figures, suffix: dseg} caption: | This panel shows the template T1-weighted image (if several T1w images were found), with contours delineating the detected brain mask and brain tissue segmentations. - subtitle: Brain mask and brain tissue segmentation of the T1w - - bids: {datatype: figures, space: .*, suffix: T1w, regex_search: True} + - subtitle: Spatial normalization of the anatomical T1w reference + bids: {datatype: figures, space: .*, suffix: T1w, regex_search: True} caption: Spatial normalization of the T1w image to the {space} template. description: | Results of nonlinear alignment of the T1w reference one or more template space(s). Hover on the panels with the mouse pointer to transition between both spaces. static: false - subtitle: Spatial normalization of the anatomical T1w reference - - bids: {datatype: figures, desc: reconall, suffix: T1w} + - subtitle: Surface reconstruction + bids: {datatype: figures, desc: reconall, suffix: T1w} caption: | Surfaces (white and pial) reconstructed with FreeSurfer (recon-all) overlaid on the participant's T1w template. - subtitle: Surface reconstruction - name: B0 field mapping ordering: session,acquisition,run,fmapid reportlets: - - bids: {datatype: figures, desc: mapped, suffix: fieldmap} + - subtitle: "Preprocessed B0 mapping acquisition" + bids: {datatype: figures, desc: mapped, suffix: fieldmap} caption: | Inhomogeneities of the B0 field introduce (oftentimes severe) spatial distortions along the phase-encoding direction of the image. Some scanners produce a B0 @@ -45,8 +46,8 @@ sections: Hover over the panels with the mouse pointer to also visualize the intensity of the field inhomogeneity in Hertz. static: false - subtitle: "Preprocessed B0 mapping acquisition" - - bids: {datatype: figures, desc: phasediff, suffix: fieldmap} + - subtitle: "Preprocessed mapping of phase-difference acquisition" + bids: {datatype: figures, desc: phasediff, suffix: fieldmap} caption: | Inhomogeneities of the B0 field introduce (oftentimes severe) spatial distortions along the phase-encoding direction of the image. A Gradient-Recalled Echo (GRE) scheme was included for the @@ -57,8 +58,8 @@ sections: Hover over the panels with the mouse pointer to also visualize the intensity of the field inhomogeneity in Hertz. static: false - subtitle: "Preprocessed mapping of phase-difference acquisition" - - bids: {datatype: figures, desc: pepolar, suffix: fieldmap} + - subtitle: "Preprocessed estimation with varying Phase-Encoding (PE) blips" + bids: {datatype: figures, desc: pepolar, suffix: fieldmap} caption: | Inhomogeneities of the B0 field introduce (oftentimes severe) spatial distortions along the phase-encoding direction of the image. Utilizing two or more images with different @@ -69,8 +70,8 @@ sections: Hover on the panels with the mouse pointer to also visualize the intensity of the inhomogeneity of the field in Hertz. static: false - subtitle: "Preprocessed estimation with varying Phase-Encoding (PE) blips" - - bids: {datatype: figures, desc: anat, suffix: fieldmap} + - subtitle: "Preprocessed estimation by nonlinear registration to an anatomical scan (“fieldmap-less”)" + bids: {datatype: figures, desc: anat, suffix: fieldmap} caption: | Inhomogeneities of the B0 field introduce (oftentimes severe) spatial distortions along the phase-encoding direction of the image. Utilizing an anatomically-correct acquisition @@ -81,14 +82,14 @@ sections: Hover on the panels with the mouse pointer to also visualize the intensity of the inhomogeneity of the field in Hertz. static: false - subtitle: "Preprocessed estimation by nonlinear registration to an anatomical scan (“fieldmap-less”)" - name: Arterial Spin Labeling ordering: session,task,acquisition,ceagent,reconstruction,direction,run,echo reportlets: - bids: {datatype: figures, desc: summary, suffix: asl} - bids: {datatype: figures, desc: validation, suffix: asl} - - bids: {datatype: figures, desc: fmapCoreg, suffix: asl} + - subtitle: Alignment between the anatomical reference of the fieldmap and the target EPI (debug mode) + bids: {datatype: figures, desc: fmapCoreg, suffix: asl} caption: | The estimated fieldmap was aligned to the corresponding EPI reference with a rigid-registration process of the anatomical reference of the fieldmap, @@ -96,8 +97,8 @@ sections: Overlaid on top of the co-registration results, the final BOLD mask is represented with a red contour for reference. static: false - subtitle: Alignment between the anatomical reference of the fieldmap and the target EPI (debug mode) - - bids: {datatype: figures, desc: fieldmap, suffix: asl} + - subtitle: "Reconstructed B0 map in the corresponding run's space (debug mode)" + bids: {datatype: figures, desc: fieldmap, suffix: asl} caption: | Estimated fieldmap, as reconstructed on the target BOLD run space to allow the assessment of its alignment with the distorted data. @@ -109,20 +110,20 @@ sections: Therefore, the fieldmap should be positioned relative to the anatomical reference exactly as it is positioned in the reportlet above. static: false - subtitle: "Reconstructed B0 map in the corresponding run's space (debug mode)" - - bids: {datatype: figures, desc: sdc, suffix: asl} + - subtitle: Susceptibility distortion correction + bids: {datatype: figures, desc: sdc, suffix: asl} caption: | Results of performing susceptibility distortion correction (SDC) on the EPI static: false - subtitle: Susceptibility distortion correction - - bids: {datatype: figures, desc: forcedsyn, suffix: asl} + - subtitle: Experimental fieldmap-less susceptibility distortion correction + bids: {datatype: figures, desc: forcedsyn, suffix: asl} caption: | The dataset contained some fieldmap information, but the argument --force-syn was used. The higher-priority SDC method was used. Here, we show the results of performing SyN-based SDC on the EPI for comparison. static: false - subtitle: Experimental fieldmap-less susceptibility distortion correction - - bids: {datatype: figures, desc: flirtnobbr, suffix: asl} + - subtitle: Alignment of functional and anatomical MRI data (volume based) + bids: {datatype: figures, desc: flirtnobbr, suffix: asl} caption: | mri_coreg (FreeSurfer) was used to generate transformations from EPI space to T1 Space - BBR refinement using FSL flirt rejected. @@ -130,8 +131,8 @@ sections: highlight potential spin-history and other artifacts, whereas final images are resampled using Lanczos interpolation. static: false - subtitle: Alignment of functional and anatomical MRI data (volume based) - - bids: {datatype: figures, desc: coreg, suffix: asl} + - subtitle: Alignment of functional and anatomical MRI data (volume based) + bids: {datatype: figures, desc: coreg, suffix: asl} caption: | mri_coreg (FreeSurfer) was used to generate transformations from EPI space to T1 Space - bbregister refinement rejected. Note @@ -139,8 +140,8 @@ sections: potential spin-history and other artifacts, whereas final images are resampled using Lanczos interpolation. static: false - subtitle: Alignment of functional and anatomical MRI data (volume based) - - bids: {datatype: figures, desc: flirtbbr, suffix: asl} + - subtitle: Alignment of functional and anatomical MRI data (surface driven) + bids: {datatype: figures, desc: flirtbbr, suffix: asl} caption: | FSL flirt was used to generate transformations from EPI-space to T1w-space - The white matter mask calculated with FSL fast (brain @@ -148,16 +149,16 @@ sections: is used in the reportlets in order to highlight potential spin-history and other artifacts, whereas final images are resampled using Lanczos interpolation. static: false - subtitle: Alignment of functional and anatomical MRI data (surface driven) - - bids: {datatype: figures, desc: bbregister, suffix: asl} + - subtitle: Alignment of functional and anatomical MRI data (surface driven) + bids: {datatype: figures, desc: bbregister, suffix: asl} caption: | bbregister was used to generate transformations from EPI-space to T1w-space. Note that Nearest Neighbor interpolation is used in the reportlets in order to highlight potential spin-history and other artifacts, whereas final images are resampled using Lanczos interpolation. static: false - subtitle: Alignment of functional and anatomical MRI data (surface driven) - - bids: {datatype: figures, desc: rois, suffix: asl} + - subtitle: Brain mask and (anatomical/temporal) CompCor ROIs + bids: {datatype: figures, desc: rois, suffix: asl} caption: | Brain mask calculated on the BOLD signal (red contour), along with the regions of interest (ROIs) used for the estimation of physiological and movement @@ -169,82 +170,120 @@ sections: variable voxels within the brain mask.
The brain edge (or crown) ROI (green contour) picks signals outside but close to the brain, which are decomposed into 24 principal components. - subtitle: Brain mask and (anatomical/temporal) CompCor ROIs - - bids: - datatype: figures - desc: '[at]compcor' - extension: [.html] - suffix: asl - - bids: {datatype: figures, desc: carpetplot, suffix: asl} + - subtitle: ASL Summary + bids: {datatype: figures, desc: carpetplot, suffix: asl} caption: | Summary statistics are plotted, which may reveal trends or artifacts in the ASL data. DVARS and FD show the standardized DVARS and framewise-displacement measures for each time point. A carpet plot shows the time series for all voxels within the brain mask. Voxels are grouped into cortical (blue), and subcortical (orange) gray matter, cerebellum (green) and white matter and CSF (red), indicated by the color map on the left-hand side. - subtitle: ASL Summary - - bids: {datatype: figures, desc: carpetplot, suffix: cbf} + - subtitle: CBF Summary + bids: {datatype: figures, desc: carpetplot, suffix: cbf} caption: | This carpet plot shows the time series for all voxels within the brain mask for CBF. Voxels are grouped into cortical (blue), and subcortical (orange) gray matter, cerebellum (green), white matter and CSF (red), indicated by the color map on the left-hand side. The SCORE index with values greater than zero indicates which volume(s) are removed by SCORE. - subtitle: CBF Summary - - bids: {datatype: figures, desc: cbf, suffix: cbf} + - subtitle: Cerebral Blood Flow + bids: {datatype: figures, desc: cbf, suffix: cbf} caption: | - The maps plot cerebral blood flow (CBF) for basic CBF. + The map plots cerebral blood flow (CBF) estimated by the basic model. The unit is mL/100 g/min. - subtitle: CBF - - bids: {datatype: figures, desc: cbfByTissueType, suffix: cbf} + - subtitle: Mean CBF by Tissue Type + bids: {datatype: figures, desc: cbfByTissueType, suffix: cbf} caption: | - The maps plot the distribution of cerebral blood flow (CBF) values for the basic output, + The map plots the distribution of cerebral blood flow (CBF) values for the basic output, separated by tissue type. The unit is mL/100 g/min. - subtitle: Mean CBF by Tissue Type - - bids: {datatype: figures, desc: score, suffix: cbf} + - subtitle: Arterial Transit Time + bids: {datatype: figures, desc: att, suffix: cbf} + caption: | + The map plots arterial transit time (ATT) estimated by the basic model. + The unit is seconds. + - subtitle: ATT by Tissue Type + bids: {datatype: figures, desc: attByTissueType, suffix: att} + caption: | + The map plots the distribution of arterial transit time values for the basic output, + separated by tissue type. + The unit is seconds. + - subtitle: Arterial Bolus Arrival Time + bids: {datatype: figures, desc: abat, suffix: cbf} + caption: | + The map plots arterial bolus arrival time (aBAT) estimated by the basic model. + The unit is seconds. + - subtitle: aBAT by Tissue Type + bids: {datatype: figures, desc: abatByTissueType, suffix: abat} + caption: | + The map plots the distribution of arterial bolus arrival time (aBAT) values for the basic output, + separated by tissue type. + The unit is seconds. + - subtitle: Arterial Blood Volume + bids: {datatype: figures, desc: abv, suffix: cbf} + caption: | + The map plots arterial blood volume (aBV) estimated by the basic model. + The unit is fraction. + - subtitle: aBV by Tissue Type + bids: {datatype: figures, desc: abvByTissueType, suffix: abv} + caption: | + The map plots the distribution of arterial blood volume (aBV) values for the basic output, + separated by tissue type. + The unit is fraction. + - subtitle: SCORE CBF + bids: {datatype: figures, desc: score, suffix: cbf} caption: | - The maps plot cerebral blood flow (CBF) for SCORE-corrected CBF. + The map plots cerebral blood flow (CBF) for SCORE-corrected CBF. The unit is mL/100 g/min. - subtitle: SCORE CBF - - bids: {datatype: figures, desc: scoreByTissueType, suffix: cbf} + - subtitle: SCORE CBF by Tissue Type + bids: {datatype: figures, desc: scoreByTissueType, suffix: cbf} caption: | - The maps plot the distribution of cerebral blood flow (CBF) values for the SCORE output, + The map plots the distribution of cerebral blood flow (CBF) values for the SCORE output, separated by tissue type. The unit is mL/100 g/min. - subtitle: SCORE CBF by Tissue Type - - bids: {datatype: figures, desc: scrub, suffix: cbf} + - subtitle: SCRUB CBF + bids: {datatype: figures, desc: scrub, suffix: cbf} caption: | - The maps plot cerebral blood flow (CBF) for SCRUB-corrected CBF. + The map plots cerebral blood flow (CBF) for SCRUB-corrected CBF. The unit is mL/100 g/min. - subtitle: SCRUB CBF - - bids: {datatype: figures, desc: scrubByTissueType, suffix: cbf} + - subtitle: SCRUB CBF by Tissue Type + bids: {datatype: figures, desc: scrubByTissueType, suffix: cbf} caption: | - The maps plot the distribution of cerebral blood flow (CBF) values for the SCRUB output, + The map plots the distribution of cerebral blood flow (CBF) values for the SCRUB output, separated by tissue type. The unit is mL/100 g/min. - subtitle: SCRUB CBF by Tissue Type - - bids: {datatype: figures, desc: basil, suffix: cbf} + - subtitle: BASIL CBF + bids: {datatype: figures, desc: basil, suffix: cbf} caption: | - The maps plot cerebral blood flow (CBF) for BASIL-estimated CBF. + The map plots cerebral blood flow (CBF) for BASIL-estimated CBF. The unit is mL/100 g/min. - subtitle: BASIL CBF - - bids: {datatype: figures, desc: basilByTissueType, suffix: cbf} + - subtitle: BASIL CBF by Tissue Type + bids: {datatype: figures, desc: basilByTissueType, suffix: cbf} caption: | - The maps plot the distribution of cerebral blood flow (CBF) values for the BASIL output, + The map plots the distribution of cerebral blood flow (CBF) values for the BASIL output, separated by tissue type. The unit is mL/100 g/min. - subtitle: BASIL CBF by Tissue Type - - bids: {datatype: figures, desc: basilGM, suffix: cbf} + - subtitle: BASIL GM-PVC CBF + bids: {datatype: figures, desc: basilGM, suffix: cbf} caption: | - The maps plot cerebral blood flow (CBF) for gray matter partial volume-corrected CBF. + The map plots cerebral blood flow (CBF) for gray matter partial volume-corrected CBF. The unit is mL/100 g/min. - subtitle: BASIL GM-PVC CBF - - bids: {datatype: figures, desc: basilGMByTissueType, suffix: cbf} + - subtitle: BASIL GM-PVC CBF by Tissue Type + bids: {datatype: figures, desc: basilGMByTissueType, suffix: cbf} caption: | - The maps plot the distribution of cerebral blood flow (CBF) values for the + The map plots the distribution of cerebral blood flow (CBF) values for the BASIL gray matter partial volume-corrected output, separated by tissue type. The unit is mL/100 g/min. - subtitle: BASIL GM-PVC CBF by Tissue Type + - subtitle: BASIL WM-PVC CBF + bids: {datatype: figures, desc: basilWM, suffix: cbf} + caption: | + The map plots cerebral blood flow (CBF) for white matter partial volume-corrected CBF. + The unit is mL/100 g/min. + - subtitle: BASIL WM-PVC CBF by Tissue Type + bids: {datatype: figures, desc: basilWMByTissueType, suffix: cbf} + caption: | + The map plots the distribution of cerebral blood flow (CBF) values for the + BASIL white matter partial volume-corrected output, separated by tissue type. + The unit is mL/100 g/min. - name: About reportlets: diff --git a/aslprep/interfaces/cbf.py b/aslprep/interfaces/cbf.py index 071b81ceb..955fd3353 100644 --- a/aslprep/interfaces/cbf.py +++ b/aslprep/interfaces/cbf.py @@ -25,12 +25,7 @@ estimate_labeling_efficiency, pcasl_or_pasl, ) -from aslprep.utils.cbf import ( - _getcbfscore, - _scrubcbf, - estimate_cbf_pcasl_multipld, - estimate_t1, -) +from aslprep.utils.cbf import _getcbfscore, _scrubcbf, estimate_t1, fit_deltam_multipld class _RefineMaskInputSpec(BaseInterfaceInputSpec): @@ -343,6 +338,16 @@ class _ComputeCBFOutputSpec(TraitedSpec): None, desc="Arterial transit time map, in seconds. Only generated for multi-delay data.", ) + abat = traits.Either( + File(exists=True), + None, + desc="Arterial bolus arrival time map, in seconds. Only generated for multi-delay data.", + ) + abv = traits.Either( + File(exists=True), + None, + desc="Arterial blood volume map. Only generated for multi-delay data.", + ) plds = traits.Either( File(exists=True), None, @@ -361,13 +366,8 @@ class ComputeCBF(SimpleInterface): Single-delay CBF, for both (P)CASL and QUIPSSII PASL is calculated according to :footcite:t:`alsop_recommended_2015`. - Multi-delay CBF is handled using a weighted average, - based on :footcite:t:`dai2012reduced,wang2013multi`. - Multi-delay CBF is calculated according to :footcite:t:`fan2017long`, - although CBF is averaged across PLDs according to the method in - :footcite:t:`juttukonda2021characterizing`. - Arterial transit time is estimated according to :footcite:t:`dai2012reduced`. + Multi-delay CBF is calculated according to :footcite:t:`woods2023recommendations`. If slice timing information is detected, then PLDs will be shifted by the slice times. @@ -377,7 +377,9 @@ class ComputeCBF(SimpleInterface): :func:`~aslprep.utils.asl.determine_multi_pld` :func:`~aslprep.utils.cbf.estimate_t1` :func:`~aslprep.utils.asl.estimate_labeling_efficiency` - :func:`~aslprep.utils.cbf.estimate_cbf_pcasl_multipld` + :func:`~aslprep.utils.cbf.calculate_deltam_pasl` + :func:`~aslprep.utils.cbf.calculate_deltam_pcasl` + :func:`~aslprep.utils.cbf.fit_deltam_multipld` References ---------- @@ -501,34 +503,52 @@ def _run_interface(self, runtime): self._results["plds"] = pld_file elif is_multi_pld: - # Broadcast PLDs to voxels by PLDs + # 3D acquisition multi-PLD + # Broadcast PLDs to voxels by PLDs, even though there's no slice timing to account for. plds = np.dot(plds[:, None], np.ones((1, deltam_arr.shape[0]))).T - if is_casl: - tau = np.array(metadata["LabelingDuration"]) - + # Now estimate CBF and any other metrics if is_multi_pld: - if is_casl: - att, mean_cbf = estimate_cbf_pcasl_multipld( - deltam_arr, - scaled_m0data, - plds, - tau, - labeleff, - t1blood=t1blood, - t1tissue=t1tissue, - unit_conversion=UNIT_CONV, - partition_coefficient=PARTITION_COEF, - ) + ti1 = None + tau = None + if is_casl: # (P)CASL needs tau, but not ti1 + tau = metadata["LabelingDuration"] + + else: # PASL needs ti1, but not tau + if metadata["BolusCutOffTechnique"] == "QUIPSSII": + # PASL + QUIPSSII + # Only one BolusCutOffDelayTime allowed. + assert isinstance(metadata["BolusCutOffDelayTime"], Number) + ti1 = metadata["BolusCutOffDelayTime"] + + elif metadata["BolusCutOffTechnique"] == "Q2TIPS": + # PASL + Q2TIPS + # Q2TIPS should have two BolusCutOffDelayTimes. + assert len(metadata["BolusCutOffDelayTime"]) == 2 + ti1 = metadata["BolusCutOffDelayTime"][0] - else: - # Dai's approach can't be used on PASL data, so we'll need another method. - raise ValueError( - "Multi-delay data are not supported for PASL sequences at the moment." - ) + else: + raise ValueError( + f"Unsupported BolusCutOffTechnique ({metadata['BolusCutOffTechnique']}) " + "for multi-PLD data." + ) + + cbf, att, abat, abv = fit_deltam_multipld( + deltam_arr=deltam_arr, + scaled_m0data=scaled_m0data, + plds=plds, + labeleff=labeleff, + t1blood=t1blood, + partition_coefficient=PARTITION_COEF, + is_casl=is_casl, + tau=tau, # defined for (P)CASL + ti1=ti1, # defined for PASL + ) - mean_cbf_img = masker.inverse_transform(mean_cbf) + mean_cbf_img = masker.inverse_transform(cbf) att_img = masker.inverse_transform(att) + abat_img = masker.inverse_transform(abat) + abv_img = masker.inverse_transform(abv) # Multi-delay data won't produce a CBF time series self._results["cbf_ts"] = None @@ -538,9 +558,22 @@ def _run_interface(self, runtime): newpath=runtime.cwd, ) att_img.to_filename(self._results["att"]) + self._results["abat"] = fname_presuffix( + self.inputs.deltam, + suffix="_abat", + newpath=runtime.cwd, + ) + abat_img.to_filename(self._results["abat"]) + self._results["abv"] = fname_presuffix( + self.inputs.deltam, + suffix="_abv", + newpath=runtime.cwd, + ) + abv_img.to_filename(self._results["abv"]) else: # Single-delay if is_casl: + tau = metadata["LabelingDuration"] denom_factor = t1blood * (1 - np.exp(-(tau / t1blood))) elif not metadata["BolusCutOffFlag"]: diff --git a/aslprep/interfaces/plotting.py b/aslprep/interfaces/plotting.py index 22c084785..ca0cde4e0 100644 --- a/aslprep/interfaces/plotting.py +++ b/aslprep/interfaces/plotting.py @@ -129,10 +129,11 @@ def _run_interface(self, runtime): class _CBFSummaryPlotInputSpec(BaseInterfaceInputSpec): - cbf = File(exists=True, mandatory=True, desc="") - label = traits.Str(exists=True, mandatory=True, desc="label") - vmax = traits.Int(exists=True, default_value=90, mandatory=True, desc="max value of asl") - ref_vol = File(exists=True, mandatory=True, desc="") + cbf = File(exists=True, mandatory=True, desc="Image to plot") + label = traits.Str(mandatory=True, desc="label") + vmin = traits.Float(mandatory=True, desc="Minimum value for figure") + vmax = traits.Float(mandatory=True, desc="Maximum value for figure") + ref_vol = File(exists=True, mandatory=True, desc="Reference volume to plot") class _CBFSummaryPlotOutputSpec(TraitedSpec): @@ -159,6 +160,7 @@ def _run_interface(self, runtime): cbf=self.inputs.cbf, label=self.inputs.label, ref_vol=self.inputs.ref_vol, + vmin=self.inputs.vmin, vmax=self.inputs.vmax, outfile=self._results["out_file"], ).plot() @@ -166,8 +168,16 @@ def _run_interface(self, runtime): class _CBFByTissueTypePlotInputSpec(BaseInterfaceInputSpec): - cbf = File(exists=True, mandatory=True, desc="") + in_file = File(exists=True, mandatory=True, desc="CBF, ATT, aBAT, or ABV file.") seg_file = File(exists=True, mandatory=True, desc="Segmentation file") + img_type = traits.Enum( + "cbf", + "att", + "abat", + "abv", + mandatory=True, + desc="Image type", + ) class _CBFByTissueTypePlotOutputSpec(TraitedSpec): @@ -186,12 +196,20 @@ def _run_interface(self, runtime): from nilearn import image, masking self._results["out_file"] = fname_presuffix( - self.inputs.cbf, - suffix="_cbfplot.svg", + self.inputs.in_file, + suffix=f"_{self.inputs.img_type}plot.svg", use_ext=False, newpath=runtime.cwd, ) + column_names = { + "cbf": "Cerebral Blood Flow\n(mL/100 g/min)", + "att": "Arterial Transit Time\n(s)", + "abat": "Arterial Bolus Arrival Time\n(s)", + "abv": "Arterial Blood Volume\n(fraction)", + } + unit_str = column_names[self.inputs.img_type] + dfs = [] for i_tissue_type, tissue_type in enumerate(["GM", "WM", "CSF"]): tissue_type_val = i_tissue_type + 1 @@ -199,9 +217,9 @@ def _run_interface(self, runtime): f"(img == {tissue_type_val}).astype(int)", img=self.inputs.seg_file, ) - tissue_type_vals = masking.apply_mask(self.inputs.cbf, mask_img) + tissue_type_vals = masking.apply_mask(self.inputs.in_file, mask_img) df = pd.DataFrame( - columns=["CBF\n(mL/100 g/min)", "Tissue Type"], + columns=[unit_str, "Tissue Type"], data=list( map(list, zip(*[tissue_type_vals, [tissue_type] * tissue_type_vals.size])) ), @@ -215,13 +233,13 @@ def _run_interface(self, runtime): fig, ax = plt.subplots(figsize=(16, 8)) sns.despine(ax=ax, bottom=True, left=True) sns.boxenplot( - y="CBF\n(mL/100 g/min)", + y=unit_str, data=df, width=0.6, showfliers=True, palette={"GM": "#1b60a5", "WM": "#2da467", "CSF": "#9d8f25"}, hue="Tissue Type", - legend=False, + legend=True, ax=ax, ) fig.tight_layout() diff --git a/aslprep/tests/data/expected_outputs_examples_pasl_multipld.txt b/aslprep/tests/data/expected_outputs_examples_pasl_multipld.txt index 4ede4c3be..f437e8f73 100644 --- a/aslprep/tests/data/expected_outputs_examples_pasl_multipld.txt +++ b/aslprep/tests/data/expected_outputs_examples_pasl_multipld.txt @@ -77,143 +77,118 @@ logs/CITATION.md logs/CITATION.tex sub-01 sub-01.html -sub-01/ses-1 -sub-01/ses-1/perf -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S1056Parcels_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S1056Parcels_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S1056Parcels_desc-basilGM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S1056Parcels_desc-basilWM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S1056Parcels_desc-basil_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S1056Parcels_desc-basil_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S1056Parcels_desc-score_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S1056Parcels_desc-scrub_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S156Parcels_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S156Parcels_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S156Parcels_desc-basilGM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S156Parcels_desc-basilWM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S156Parcels_desc-basil_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S156Parcels_desc-basil_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S156Parcels_desc-score_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S156Parcels_desc-scrub_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S256Parcels_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S256Parcels_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S256Parcels_desc-basilGM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S256Parcels_desc-basilWM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S256Parcels_desc-basil_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S256Parcels_desc-basil_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S256Parcels_desc-score_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S256Parcels_desc-scrub_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S356Parcels_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S356Parcels_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S356Parcels_desc-basilGM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S356Parcels_desc-basilWM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S356Parcels_desc-basil_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S356Parcels_desc-basil_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S356Parcels_desc-score_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S356Parcels_desc-scrub_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S456Parcels_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S456Parcels_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S456Parcels_desc-basilGM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S456Parcels_desc-basilWM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S456Parcels_desc-basil_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S456Parcels_desc-basil_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S456Parcels_desc-score_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S456Parcels_desc-scrub_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S556Parcels_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S556Parcels_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S556Parcels_desc-basilGM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S556Parcels_desc-basilWM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S556Parcels_desc-basil_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S556Parcels_desc-basil_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S556Parcels_desc-score_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S556Parcels_desc-scrub_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S656Parcels_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S656Parcels_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S656Parcels_desc-basilGM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S656Parcels_desc-basilWM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S656Parcels_desc-basil_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S656Parcels_desc-basil_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S656Parcels_desc-score_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S656Parcels_desc-scrub_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S756Parcels_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S756Parcels_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S756Parcels_desc-basilGM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S756Parcels_desc-basilWM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S756Parcels_desc-basil_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S756Parcels_desc-basil_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S756Parcels_desc-score_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S756Parcels_desc-scrub_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S856Parcels_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S856Parcels_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S856Parcels_desc-basilGM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S856Parcels_desc-basilWM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S856Parcels_desc-basil_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S856Parcels_desc-basil_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S856Parcels_desc-score_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S856Parcels_desc-scrub_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S956Parcels_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S956Parcels_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S956Parcels_desc-basilGM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S956Parcels_desc-basilWM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S956Parcels_desc-basil_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S956Parcels_desc-basil_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S956Parcels_desc-score_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-4S956Parcels_desc-scrub_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Glasser_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Glasser_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Glasser_desc-basilGM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Glasser_desc-basilWM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Glasser_desc-basil_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Glasser_desc-basil_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Glasser_desc-score_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Glasser_desc-scrub_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Gordon_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Gordon_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Gordon_desc-basilGM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Gordon_desc-basilWM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Gordon_desc-basil_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Gordon_desc-basil_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Gordon_desc-score_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Gordon_desc-scrub_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-HCP_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-HCP_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-HCP_desc-basilGM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-HCP_desc-basilWM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-HCP_desc-basil_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-HCP_desc-basil_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-HCP_desc-score_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-HCP_desc-scrub_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Tian_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Tian_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Tian_desc-basilGM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Tian_desc-basilWM_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Tian_desc-basil_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Tian_desc-basil_coverage.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Tian_desc-score_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_atlas-Tian_desc-scrub_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_att.json -sub-01/ses-1/perf/sub-01_ses-1_att.nii.gz -sub-01/ses-1/perf/sub-01_ses-1_cbf.json -sub-01/ses-1/perf/sub-01_ses-1_cbf.nii.gz -sub-01/ses-1/perf/sub-01_ses-1_desc-basilGM_cbf.json -sub-01/ses-1/perf/sub-01_ses-1_desc-basilGM_cbf.nii.gz -sub-01/ses-1/perf/sub-01_ses-1_desc-basilWM_cbf.json -sub-01/ses-1/perf/sub-01_ses-1_desc-basilWM_cbf.nii.gz -sub-01/ses-1/perf/sub-01_ses-1_desc-basil_att.json -sub-01/ses-1/perf/sub-01_ses-1_desc-basil_att.nii.gz -sub-01/ses-1/perf/sub-01_ses-1_desc-basil_cbf.json -sub-01/ses-1/perf/sub-01_ses-1_desc-basil_cbf.nii.gz -sub-01/ses-1/perf/sub-01_ses-1_desc-brain_mask.json -sub-01/ses-1/perf/sub-01_ses-1_desc-brain_mask.nii.gz -sub-01/ses-1/perf/sub-01_ses-1_desc-confounds_timeseries.tsv -sub-01/ses-1/perf/sub-01_ses-1_desc-coreg_aslref.json -sub-01/ses-1/perf/sub-01_ses-1_desc-coreg_aslref.nii.gz -sub-01/ses-1/perf/sub-01_ses-1_desc-hmc_aslref.json -sub-01/ses-1/perf/sub-01_ses-1_desc-hmc_aslref.nii.gz -sub-01/ses-1/perf/sub-01_ses-1_desc-preproc_asl.json -sub-01/ses-1/perf/sub-01_ses-1_desc-preproc_asl.nii.gz -sub-01/ses-1/perf/sub-01_ses-1_desc-qualitycontrol_cbf.tsv -sub-01/ses-1/perf/sub-01_ses-1_from-aslref_to-T1w_mode-image_xfm.json -sub-01/ses-1/perf/sub-01_ses-1_from-aslref_to-T1w_mode-image_xfm.txt -sub-01/ses-1/perf/sub-01_ses-1_from-orig_to-aslref_mode-image_xfm.json -sub-01/ses-1/perf/sub-01_ses-1_from-orig_to-aslref_mode-image_xfm.txt +sub-01/perf +sub-01/perf/sub-01_abat.json +sub-01/perf/sub-01_abat.nii.gz +sub-01/perf/sub-01_abv.json +sub-01/perf/sub-01_abv.nii.gz +sub-01/perf/sub-01_atlas-4S1056Parcels_cbf.tsv +sub-01/perf/sub-01_atlas-4S1056Parcels_coverage.tsv +sub-01/perf/sub-01_atlas-4S1056Parcels_desc-basilGM_cbf.tsv +sub-01/perf/sub-01_atlas-4S1056Parcels_desc-basilWM_cbf.tsv +sub-01/perf/sub-01_atlas-4S1056Parcels_desc-basil_cbf.tsv +sub-01/perf/sub-01_atlas-4S1056Parcels_desc-basil_coverage.tsv +sub-01/perf/sub-01_atlas-4S156Parcels_cbf.tsv +sub-01/perf/sub-01_atlas-4S156Parcels_coverage.tsv +sub-01/perf/sub-01_atlas-4S156Parcels_desc-basilGM_cbf.tsv +sub-01/perf/sub-01_atlas-4S156Parcels_desc-basilWM_cbf.tsv +sub-01/perf/sub-01_atlas-4S156Parcels_desc-basil_cbf.tsv +sub-01/perf/sub-01_atlas-4S156Parcels_desc-basil_coverage.tsv +sub-01/perf/sub-01_atlas-4S256Parcels_cbf.tsv +sub-01/perf/sub-01_atlas-4S256Parcels_coverage.tsv +sub-01/perf/sub-01_atlas-4S256Parcels_desc-basilGM_cbf.tsv +sub-01/perf/sub-01_atlas-4S256Parcels_desc-basilWM_cbf.tsv +sub-01/perf/sub-01_atlas-4S256Parcels_desc-basil_cbf.tsv +sub-01/perf/sub-01_atlas-4S256Parcels_desc-basil_coverage.tsv +sub-01/perf/sub-01_atlas-4S356Parcels_cbf.tsv +sub-01/perf/sub-01_atlas-4S356Parcels_coverage.tsv +sub-01/perf/sub-01_atlas-4S356Parcels_desc-basilGM_cbf.tsv +sub-01/perf/sub-01_atlas-4S356Parcels_desc-basilWM_cbf.tsv +sub-01/perf/sub-01_atlas-4S356Parcels_desc-basil_cbf.tsv +sub-01/perf/sub-01_atlas-4S356Parcels_desc-basil_coverage.tsv +sub-01/perf/sub-01_atlas-4S456Parcels_cbf.tsv +sub-01/perf/sub-01_atlas-4S456Parcels_coverage.tsv +sub-01/perf/sub-01_atlas-4S456Parcels_desc-basilGM_cbf.tsv +sub-01/perf/sub-01_atlas-4S456Parcels_desc-basilWM_cbf.tsv +sub-01/perf/sub-01_atlas-4S456Parcels_desc-basil_cbf.tsv +sub-01/perf/sub-01_atlas-4S456Parcels_desc-basil_coverage.tsv +sub-01/perf/sub-01_atlas-4S556Parcels_cbf.tsv +sub-01/perf/sub-01_atlas-4S556Parcels_coverage.tsv +sub-01/perf/sub-01_atlas-4S556Parcels_desc-basilGM_cbf.tsv +sub-01/perf/sub-01_atlas-4S556Parcels_desc-basilWM_cbf.tsv +sub-01/perf/sub-01_atlas-4S556Parcels_desc-basil_cbf.tsv +sub-01/perf/sub-01_atlas-4S556Parcels_desc-basil_coverage.tsv +sub-01/perf/sub-01_atlas-4S656Parcels_cbf.tsv +sub-01/perf/sub-01_atlas-4S656Parcels_coverage.tsv +sub-01/perf/sub-01_atlas-4S656Parcels_desc-basilGM_cbf.tsv +sub-01/perf/sub-01_atlas-4S656Parcels_desc-basilWM_cbf.tsv +sub-01/perf/sub-01_atlas-4S656Parcels_desc-basil_cbf.tsv +sub-01/perf/sub-01_atlas-4S656Parcels_desc-basil_coverage.tsv +sub-01/perf/sub-01_atlas-4S756Parcels_cbf.tsv +sub-01/perf/sub-01_atlas-4S756Parcels_coverage.tsv +sub-01/perf/sub-01_atlas-4S756Parcels_desc-basilGM_cbf.tsv +sub-01/perf/sub-01_atlas-4S756Parcels_desc-basilWM_cbf.tsv +sub-01/perf/sub-01_atlas-4S756Parcels_desc-basil_cbf.tsv +sub-01/perf/sub-01_atlas-4S756Parcels_desc-basil_coverage.tsv +sub-01/perf/sub-01_atlas-4S856Parcels_cbf.tsv +sub-01/perf/sub-01_atlas-4S856Parcels_coverage.tsv +sub-01/perf/sub-01_atlas-4S856Parcels_desc-basilGM_cbf.tsv +sub-01/perf/sub-01_atlas-4S856Parcels_desc-basilWM_cbf.tsv +sub-01/perf/sub-01_atlas-4S856Parcels_desc-basil_cbf.tsv +sub-01/perf/sub-01_atlas-4S856Parcels_desc-basil_coverage.tsv +sub-01/perf/sub-01_atlas-4S956Parcels_cbf.tsv +sub-01/perf/sub-01_atlas-4S956Parcels_coverage.tsv +sub-01/perf/sub-01_atlas-4S956Parcels_desc-basilGM_cbf.tsv +sub-01/perf/sub-01_atlas-4S956Parcels_desc-basilWM_cbf.tsv +sub-01/perf/sub-01_atlas-4S956Parcels_desc-basil_cbf.tsv +sub-01/perf/sub-01_atlas-4S956Parcels_desc-basil_coverage.tsv +sub-01/perf/sub-01_atlas-Glasser_cbf.tsv +sub-01/perf/sub-01_atlas-Glasser_coverage.tsv +sub-01/perf/sub-01_atlas-Glasser_desc-basilGM_cbf.tsv +sub-01/perf/sub-01_atlas-Glasser_desc-basilWM_cbf.tsv +sub-01/perf/sub-01_atlas-Glasser_desc-basil_cbf.tsv +sub-01/perf/sub-01_atlas-Glasser_desc-basil_coverage.tsv +sub-01/perf/sub-01_atlas-Gordon_cbf.tsv +sub-01/perf/sub-01_atlas-Gordon_coverage.tsv +sub-01/perf/sub-01_atlas-Gordon_desc-basilGM_cbf.tsv +sub-01/perf/sub-01_atlas-Gordon_desc-basilWM_cbf.tsv +sub-01/perf/sub-01_atlas-Gordon_desc-basil_cbf.tsv +sub-01/perf/sub-01_atlas-Gordon_desc-basil_coverage.tsv +sub-01/perf/sub-01_atlas-HCP_cbf.tsv +sub-01/perf/sub-01_atlas-HCP_coverage.tsv +sub-01/perf/sub-01_atlas-HCP_desc-basilGM_cbf.tsv +sub-01/perf/sub-01_atlas-HCP_desc-basilWM_cbf.tsv +sub-01/perf/sub-01_atlas-HCP_desc-basil_cbf.tsv +sub-01/perf/sub-01_atlas-HCP_desc-basil_coverage.tsv +sub-01/perf/sub-01_atlas-Tian_cbf.tsv +sub-01/perf/sub-01_atlas-Tian_coverage.tsv +sub-01/perf/sub-01_atlas-Tian_desc-basilGM_cbf.tsv +sub-01/perf/sub-01_atlas-Tian_desc-basilWM_cbf.tsv +sub-01/perf/sub-01_atlas-Tian_desc-basil_cbf.tsv +sub-01/perf/sub-01_atlas-Tian_desc-basil_coverage.tsv +sub-01/perf/sub-01_att.json +sub-01/perf/sub-01_att.nii.gz +sub-01/perf/sub-01_cbf.json +sub-01/perf/sub-01_cbf.nii.gz +sub-01/perf/sub-01_desc-basilGM_cbf.json +sub-01/perf/sub-01_desc-basilGM_cbf.nii.gz +sub-01/perf/sub-01_desc-basilWM_cbf.json +sub-01/perf/sub-01_desc-basilWM_cbf.nii.gz +sub-01/perf/sub-01_desc-basil_att.json +sub-01/perf/sub-01_desc-basil_att.nii.gz +sub-01/perf/sub-01_desc-basil_cbf.json +sub-01/perf/sub-01_desc-basil_cbf.nii.gz +sub-01/perf/sub-01_desc-brain_mask.json +sub-01/perf/sub-01_desc-brain_mask.nii.gz +sub-01/perf/sub-01_desc-confounds_timeseries.tsv +sub-01/perf/sub-01_desc-coreg_aslref.json +sub-01/perf/sub-01_desc-coreg_aslref.nii.gz +sub-01/perf/sub-01_desc-hmc_aslref.json +sub-01/perf/sub-01_desc-hmc_aslref.nii.gz +sub-01/perf/sub-01_desc-preproc_asl.json +sub-01/perf/sub-01_desc-preproc_asl.nii.gz +sub-01/perf/sub-01_desc-qualitycontrol_cbf.tsv +sub-01/perf/sub-01_from-aslref_to-T1w_mode-image_xfm.json +sub-01/perf/sub-01_from-aslref_to-T1w_mode-image_xfm.txt +sub-01/perf/sub-01_from-orig_to-aslref_mode-image_xfm.json +sub-01/perf/sub-01_from-orig_to-aslref_mode-image_xfm.txt diff --git a/aslprep/tests/data/expected_outputs_examples_pcasl_multipld.txt b/aslprep/tests/data/expected_outputs_examples_pcasl_multipld.txt index b76e98ae7..12487f0bf 100644 --- a/aslprep/tests/data/expected_outputs_examples_pcasl_multipld.txt +++ b/aslprep/tests/data/expected_outputs_examples_pcasl_multipld.txt @@ -164,6 +164,10 @@ sub-01/perf/sub-01_atlas-Tian_desc-basil_cbf.tsv sub-01/perf/sub-01_atlas-Tian_desc-basil_coverage.tsv sub-01/perf/sub-01_att.json sub-01/perf/sub-01_att.nii.gz +sub-01/perf/sub-01_abat.json +sub-01/perf/sub-01_abat.nii.gz +sub-01/perf/sub-01_abv.json +sub-01/perf/sub-01_abv.nii.gz sub-01/perf/sub-01_cbf.json sub-01/perf/sub-01_cbf.nii.gz sub-01/perf/sub-01_desc-basilGM_cbf.json diff --git a/aslprep/tests/test_cli.py b/aslprep/tests/test_cli.py index 9755f4a98..fcb495774 100644 --- a/aslprep/tests/test_cli.py +++ b/aslprep/tests/test_cli.py @@ -47,7 +47,6 @@ def test_examples_pasl_multipld(data_dir, output_dir, working_dir): "--output-spaces=asl", "--scorescrub", "--basil", - "--use-syn-sdc", "--m0_scale=10", "--fs-no-reconall", f"--fs-subjects-dir={os.path.join(data_dir, 'anatomical/freesurfer')}", @@ -55,7 +54,7 @@ def test_examples_pasl_multipld(data_dir, output_dir, working_dir): f"{os.path.join(data_dir, 'anatomical/smriprep')}", ] - _run_and_fail(parameters) + _run_and_generate(TEST_NAME, PARTICIPANT_LABEL, parameters, out_dir) @pytest.mark.examples_pcasl_multipld diff --git a/aslprep/tests/test_interfaces_cbf.py b/aslprep/tests/test_interfaces_cbf.py index 6216bc4d6..2bd132488 100644 --- a/aslprep/tests/test_interfaces_cbf.py +++ b/aslprep/tests/test_interfaces_cbf.py @@ -223,7 +223,7 @@ def test_computecbf_pasl(datasets, tmp_path_factory): m0_file=m0_file, mask=mask_file, ) - with pytest.raises(ValueError, match="Multi-delay data are not supported"): + with pytest.raises(ValueError, match="Unsupported BolusCutOffTechnique"): results = interface.run(cwd=tmpdir) # Scenario 4: QUIPSS PASL with one PostLabelingDelay for each deltam volume (good) @@ -243,7 +243,7 @@ def test_computecbf_pasl(datasets, tmp_path_factory): m0_file=m0_file, mask=mask_file, ) - with pytest.raises(ValueError, match="Multi-delay data are not supported"): + with pytest.raises(ValueError, match="Unsupported BolusCutOffTechnique"): results = interface.run(cwd=tmpdir) # Scenario 5: QUIPSSII PASL with one PostLabelingDelay @@ -290,8 +290,7 @@ def test_computecbf_pasl(datasets, tmp_path_factory): m0_file=m0_file, mask=mask_file, ) - with pytest.raises(ValueError, match="Multi-delay data are not supported"): - results = interface.run(cwd=tmpdir) + results = interface.run(cwd=tmpdir) # Scenario 7: Q2TIPS PASL with one PostLabelingDelay # This should produce CBF time series and mean CBF, but no ATT @@ -337,8 +336,7 @@ def test_computecbf_pasl(datasets, tmp_path_factory): m0_file=m0_file, mask=mask_file, ) - with pytest.raises(ValueError, match="Multi-delay data are not supported"): - results = interface.run(cwd=tmpdir) + results = interface.run(cwd=tmpdir) def test_compare_slicetiming(datasets, tmp_path_factory): diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index b24e33e8b..df3320625 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -3,7 +3,6 @@ from __future__ import annotations import numpy as np -from scipy.interpolate import interp1d from scipy.stats import median_abs_deviation from aslprep import config @@ -486,337 +485,437 @@ def _scrubcbf(cbf_ts, gm, wm, csf, mask, wfun="huber", thresh=0.7): return newcbf -def estimate_att_pcasl(deltam_arr, plds, lds, t1blood, t1tissue): - """Estimate arterial transit time using the weighted average method. +def estimate_t1(metadata): + """Estimate the relaxation rates of blood and gray matter based on magnetic field strength. + + t1blood is set based on the scanner's field strength, + according to :footcite:t:`zhang2013vivo,alsop_recommended_2015`. + If recommended values from these publications cannot be used + (i.e., if the field strength isn't 1.5T, 3T, 7T), + then the formula from :footcite:t:`zhang2013vivo` will be applied. - The weighted average method comes from :footcite:t:`dai2012reduced`. + t1tissue is set based on the scanner's field strength as well, + according to :footcite:t:`wright2008water`. + At the moment, field strengths other than 1.5T, 3T, and 7T are not supported and will + raise an exception. Parameters ---------- - deltam_arr : :obj:`numpy.ndarray` of shape (S, D) - Delta-M array, averaged by PLD. - plds : :obj:`numpy.ndarray` of shape (S, D) - Post-labeling delays. w in Dai 2012. - In case of a 2D acquisition, PLDs may vary by slice, and thus the plds array will vary - in the spatial dimension. - For 3D acquisitions, or 2D acquisitions without slice timing info, plds will only vary - along the second dimension. - lds : :obj:`numpy.ndarray` - Labeling durations. tau in Dai 2012. + metadata : :obj:`dict` + Dictionary of metadata from the ASL file. + + Returns + ------- t1blood : :obj:`float` - T1 relaxation rate for blood. + Estimated relaxation rate of blood based on magnetic field strength. t1tissue : :obj:`float` - T1 relaxation rate for tissue. + Estimated relaxation rate of gray matter based on magnetic field strength. + + References + ---------- + .. footbibliography:: + """ + T1BLOOD_DICT = { + 1.5: 1.35, + 3: 1.65, + 7: 2.087, + } + t1blood = T1BLOOD_DICT.get(metadata["MagneticFieldStrength"]) + if not t1blood: + config.loggers.interface.warning( + f"T1blood cannot be inferred for {metadata['MagneticFieldStrength']}T data. " + "Defaulting to formula from Zhang et al. (2013)." + ) + t1blood = (110 * metadata["MagneticFieldStrength"] + 1316) / 1000 + + # TODO: Supplement with formula for other field strengths + T1TISSUE_DICT = { + 1.5: 1.197, + 3: 1.607, + 7: 1.939, + } + t1tissue = T1TISSUE_DICT.get(metadata["MagneticFieldStrength"]) + if not t1tissue: + raise ValueError( + f"T1tissue cannot be inferred for {metadata['MagneticFieldStrength']}T data." + ) + + return t1blood, t1tissue + + +def calculate_deltam_pcasl(X, cbf, att, abat, abv): + """Specify a model for use with scipy curve fitting. + + The model is fit separately for each voxel. + + This model is used to estimate ``cbf``, ``att``, ``abat``, and ``abv`` for multi-PLD PCASL + data. + + Parameters + ---------- + cbf : :obj:`float` + Cerebral blood flow. + Used for deltam_tiss calculation. + Estimated by the model. + att : :obj:`float` + Arterial transit time. + Used for deltam_tiss calculation. + Estimated by the model. + abat : :obj:`float` + The arterial bolus arrival time, describing the first arrival of the labeled blood + water in the voxel within the arterial vessel. Used for deltam_art calculation. + Estimated by the model. + abv : :obj:`float` + Arterial blood volume. + The fraction of the voxel occupied by the arterial vessel. + Used for deltam_art calculation. + Estimated by the model. + X : :obj:`numpy.ndarray` of shape (n_plds, 6) + Dependent variables. Returns ------- - att_arr : :obj:`numpy.ndarray` - Arterial transit time array. + deltam : :obj:`numpy.ndarray` + Delta-M values for the voxel. Notes ----- - This function was originally written in MATLAB by Jianxun Qu and William Tackett. - It was translated to Python by Taylor Salo. - Taylor Salo modified the code to loop over voxels, in order to account for slice timing-shifted - post-labeling delays. + X breaks down into six variables: - Please see https://shorturl.at/wCO56 and https://shorturl.at/aKQU3 for the original MATLAB - code. + labeleff : :obj:`float` + Labeling efficiency. + If background suppression was performed, then this is already adjusted and combines + alpha_bs as well. + Used for both deltam_tiss and deltam_art calculation. + t1blood : :obj:`float` + Longitudinal relaxation time of arterial blood in seconds. + Used for both deltam_tiss and deltam_art calculation. + m0_a : :obj:`float` + Equilibrium magnetization of tissue, calculated as M0a = SIpd / lambda, + where SIpd is a proton density weighted image and lambda is the tissue/blood partition + coefficient. + Used for deltam_tiss calculation. + m0_b : :obj:`float` + Equilibrium magnetization of arterial blood. + Used for deltam_art calculation. + plds : :obj:`numpy.ndarray` + Post-labeling delays. + Used for both deltam_tiss and deltam_art calculation. + taus : :obj:`numpy.ndarray` + Labeling durations. Also known as LD. + Used for both deltam_tiss and deltam_art calculation. + """ + # float parameters + labeleff = X[0, 0] + t1blood = X[0, 1] + m0_a = X[0, 2] + m0_b = X[0, 3] + # array parameters + plds = X[:, 4] + taus = X[:, 5] + + deltam = np.zeros(plds.size) + + for i_pld, pld in enumerate(plds): + tau = taus[i_pld] + + # Equation 2: the tissue component + if 0 < (tau + pld) < att: + # 0 < LD + PLD < ATT + deltam_tiss = 0 + elif att <= (tau + pld) < (att + tau): + # ATT < LD + PLD < ATT + LD + deltam_tiss = ( + 2 + * labeleff + * t1blood + * m0_a + * cbf + * np.exp(-att / t1blood) + * (1 - np.exp(-(tau + pld - att) / t1blood)) + ) / 6000 + else: + # 0 < ATT < PLD + deltam_tiss = ( + 2 + * labeleff + * t1blood + * m0_a + * cbf + * np.exp(-pld / t1blood) + * (1 - np.exp(-tau / t1blood)) + ) / 6000 + + # Equation 4: the intravascular component + if 0 < (tau + pld) < abat: + # 0 < LD + PLD < aBAT + deltam_art = 0 + elif abat <= (tau + pld) < (abat + tau): + # aBAT < LD + PLD < aBAT + LD + deltam_art = 2 * labeleff * m0_b * abv * np.exp(-abat / t1blood) + else: + # aBAT < PLD + deltam_art = 0 - The code could probably be improved by operating on arrays, rather than looping over voxels. - It is also overkill for 3D acquisitions, where PLD doesn't vary by voxel. + deltam[i_pld] = deltam_tiss + deltam_art - License - ------- - MIT License + return deltam - Copyright (c) 2023 willtack - Permission is hereby granted, free of charge, to any person obtaining a copy - of this software and associated documentation files (the "Software"), to deal - in the Software without restriction, including without limitation the rights - to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - copies of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: +def calculate_deltam_pasl(X, cbf, att, abat, abv): + """Specify a model for use with scipy curve fitting. - The above copyright notice and this permission notice shall be included in all - copies or substantial portions of the Software. + The model is fit separately for each voxel. - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - SOFTWARE. + This model is used to estimate ``cbf``, ``att``, ``abat``, and ``abv`` for multi-PLD PCASL + data. - References + Parameters ---------- - .. footbibliography:: - """ - n_voxels, n_plds = plds.shape - assert deltam_arr.shape == plds.shape, f"{deltam_arr.shape} != {plds.shape}" + cbf : :obj:`float` + Cerebral blood flow. + Used for deltam_tiss calculation. + Estimated by the model. + att : :obj:`float` + Arterial transit time. + Used for deltam_tiss calculation. + Estimated by the model. + abat : :obj:`float` + The arterial bolus arrival time, describing the first arrival of the labeled blood + water in the voxel within the arterial vessel. Used for deltam_art calculation. + Estimated by the model. + abv : :obj:`float` + Arterial blood volume. + The fraction of the voxel occupied by the arterial vessel. + Used for deltam_art calculation. + Estimated by the model. + X : :obj:`numpy.ndarray` of shape (n_plds/tis, 6) + Dependent variables. - # Beginning of auxil_asl_gen_wsum - assert lds.size == n_plds, f"{lds.size} != {n_plds}" - - att_arr = np.empty(n_voxels) - for i_voxel in range(n_voxels): - deltam_by_voxel = deltam_arr[i_voxel, :] - plds_by_voxel = plds[i_voxel, :] - - # Define the possible transit times to evaluate for this voxel - transit_times = np.arange( - np.round(np.min(plds_by_voxel), 3), - np.round(np.max(plds_by_voxel), 3) + 0.001, - 0.001, - ) - - sig_sum = np.zeros((transit_times.size)) - sig_pld_sum = np.zeros((transit_times.size)) - - for j_pld in range(n_plds): - pld = plds_by_voxel[j_pld] - ld = lds[j_pld] - - # e ^ (-delta / T1a) - exp_tt_T1b = np.exp(-transit_times / t1blood) - # e ^ (-max(w - delta, 0) / T1t) - exp_pld_tt_T1t = np.exp(-np.maximum(0, pld - transit_times) / t1tissue) - # e ^ (-max(tau + w - delta, 0) / T1t) - exp_pld_ld_tt_T1t = np.exp(-np.maximum(0, pld + ld - transit_times) / t1tissue) - - # The combination of exponent terms in Equation 1 - sig = exp_tt_T1b * (exp_pld_tt_T1t - exp_pld_ld_tt_T1t) - - # The elements in Equation 1 that aren't touched here are: - # 2M0t (2 * equilibrium magnetization of brain tissue), - # Beta (a term to compensate for static tissue signal loss caused by vessel supp. - # pulses), - # alpha (labeling efficiency), - # T1t (t1tissue), - # f (CBF; perfusion rate) - # lambda (partition coefficient) - # It seems like they are cancelled out though, so it's probably fine. - - # Numerator in Equation 4, for a range of transit times - sig_pld_sum += sig * pld - # Denominator in Equation 4, for a range of transit times - sig_sum += sig - - # Predicted weighted delay values for a range of transit times - weighted_delay_predicted = sig_pld_sum / sig_sum # TT - # End of auxil_asl_gen_wsum - - # Calculate the observed weighted delay for each voxel - weighted_delay_denom = np.sum(deltam_by_voxel) - weighted_delay_num = np.sum(deltam_by_voxel * plds_by_voxel) - weighted_delay_observed = weighted_delay_num / ( - np.abs(weighted_delay_denom) + np.finfo(float).eps - ) + Returns + ------- + deltam : :obj:`float` + Delta-M value for the voxel. - # Truncate extreme transit time value to the PLD limits - weighted_delay_min = min(weighted_delay_predicted) - weighted_delay_max = max(weighted_delay_predicted) - weighted_delay_observed = np.maximum(weighted_delay_observed, weighted_delay_min) - weighted_delay_observed = np.minimum(weighted_delay_observed, weighted_delay_max) + Notes + ----- + X breaks down into six variables: - # Use linear interpolation to get the ATT for each weighted delay value (i.e., each voxel), - # using the predicted weighted delay and associated transit time arrays. - interp_func = interp1d(weighted_delay_predicted, transit_times) + labeleff : :obj:`float` + Labeling efficiency. + In the case of background suppression, this combines alpha and alpha_bs, + so we don't need a separate term for alpha_bs. + Used for both deltam_tiss and deltam_art calculation. + t1blood : :obj:`float` + Longitudinal relaxation time of arterial blood in seconds. + Used for both deltam_tiss and deltam_art calculation. + m0_a : :obj:`float` + Equilibrium magnetization of tissue, calculated as M0a = SIpd / lambda, + where SIpd is a proton density weighted image and lambda is the tissue/blood partition + coefficient. + SIpd is also already scaled by any user-provided scaling factor. + Used for deltam_tiss calculation. + m0_b : :obj:`float` + Equilibrium magnetization of arterial blood. + Used for deltam_art calculation. + tis : :obj:`numpy.ndarray` + Inversion times in seconds. + ti1 : :obj:`float` + Bolus cutoff delay time. For Q2TIPS, this is the first delay time. + For QUIPSSII, it's the only one. + """ + # float parameters + labeleff = X[0, 0] + t1blood = X[0, 1] + m0_a = X[0, 2] + m0_b = X[0, 3] + tis = X[:, 4] # array parameter + ti1 = X[0, 5] + + deltam = np.zeros(tis.size) + + for i_ti, ti in enumerate(tis): + # Equation 3: the tissue component + if 0 < ti < att: + # 0 < TI < ATT + deltam_tiss = 0 + elif att < ti < (att + ti1): + # ATT < TI < ATT + TI1 + deltam_tiss = (2 * labeleff * m0_a * cbf * np.exp(-ti / t1blood) * (ti - att)) / 6000 + else: + # ATT + TI1 < TI + deltam_tiss = (2 * labeleff * m0_a * cbf * np.exp(-ti / t1blood) * ti1) / 6000 + + # Equation 5: the intravascular component + if 0 < ti < abat: + # 0 < TI < aBAT + deltam_art = 0 + elif abat < ti < (abat + ti1): + # aBAT < TI < aBAT + TI1 + deltam_art = 2 * labeleff * m0_b * abv * np.exp(-ti / t1blood) + else: + # aBAT + TI1 < TI + deltam_art = 0 - att_arr[i_voxel] = interp_func(weighted_delay_observed) + deltam[i_ti] = deltam_tiss + deltam_art - return att_arr + return deltam -def estimate_cbf_pcasl_multipld( +def fit_deltam_multipld( deltam_arr, scaled_m0data, plds, - tau, labeleff, t1blood, - t1tissue, - unit_conversion, partition_coefficient, + is_casl, + tau=None, + ti1=None, ): - """Estimate CBF and ATT for multi-delay PCASL data. + """Estimate CBF, ATT, aBAT, and aBV from multi-PLD PCASL data. + + This function uses the general kinetic model specified by + :footcite:t:`woods2023recommendations`. Parameters ---------- - deltam_arr : :obj:`numpy.ndarray` of shape (S, P) - Control-label values for each voxel and PLDs. - S = sample (i.e., voxel). - P = Post-labeling delay (i.e., volume). - scaled_m0data : :obj:`numpy.ndarray` of shape (S,) - The M0 volume, after scaling based on the M0-scale value. - plds : :obj:`numpy.ndarray` of shape (S, P) - Post-labeling delays. One value for each volume in ``deltam_arr``. - tau : :obj:`numpy.ndarray` of shape (P,) or (0,) - Label duration. May be a single value or may vary across volumes/PLDs. + deltam_arr : :obj:`numpy.ndarray` of shape (n_voxels, n_volumes) + Subtracted (control - label) signal. + scaled_m0data : :obj:`numpy.ndarray` of shape (n_voxels,) + M0 data, after any scaling requested by the user. + plds : :obj:`numpy.ndarray` of shape (n_voxels, n_volumes) + Post-labeling delay values. This uses the TI values for PASL data. labeleff : :obj:`float` - Estimated labeling efficiency. + Labeling efficiency, also referred to as alpha. + In the case of background suppression, this combines alpha and alpha_bs, + so we don't need a separate term for alpha_bs. t1blood : :obj:`float` - The longitudinal relaxation time of blood in seconds. - t1tissue : :obj:`float` - The longitudinal relaxation time of tissue in seconds. - unit_conversion : :obj:`float` - The factor to convert CBF units from mL/g/s to mL/ (100 g)/min. 6000. + Longitudinal relaxation time of arterial blood in seconds. + Used for both deltam_tiss and deltam_art calculation. partition_coefficient : :obj:`float` - The brain/blood partition coefficient in mL/g. Called lambda in the literature. + Brain partition coefficient (lambda in Alsop 2015). Generally 0.9. + is_casl : :obj:`bool` + True if is (P)CASL sequence. + tau : :obj:`float` or :obj:`numpy.ndarray` or None, optional + Label duration. May be a single value or may vary across volumes/PLDs. + Only defined for (P)CASL data. + ti1 : :obj:`float` or None, optional + TI1. Only defined for PASL data. Returns ------- - att_arr : :obj:`numpy.ndarray` of shape (S,) - Arterial transit time map. - cbf : :obj:`numpy.ndarray` of shape (S,) - Estimated cerebrospinal fluid map, after estimating for each PLD and averaging across - PLDs. + cbf + att + abat + abv Notes ----- - Delta-M values are first averaged over time for each unique post-labeling delay value. - - Arterial transit time is estimated on a voxel-wise basis according to - :footcite:t:`dai2012reduced`. - - CBF is then calculated for each unique PLD value using the mean delta-M values and the - estimated ATT. - - CBF is then averaged over PLDs according to :footcite:t:`juttukonda2021characterizing`, - in which an unweighted average is calculated for each voxel across all PLDs in which - PLD + tau > ATT. + This model does not include a separate term for T1,tissue, but we could expand it to do so + in the future. References ---------- .. footbibliography:: """ - n_voxels, n_volumes = deltam_arr.shape - n_voxels_pld, n_plds = plds.shape - assert n_voxels_pld == n_voxels - first_voxel_plds = plds[0, :] + import scipy - if n_plds != n_volumes: + n_voxels, n_volumes = deltam_arr.shape + if deltam_arr.shape != plds.shape: raise ValueError( - f"Number of PostLabelingDelays ({n_plds}) does not match " - f"number of delta-M volumes ({n_volumes})." + f"Number of PostLabelingDelays ({plds.shape}) does not match " + f"number of delta-M volumes ({deltam_arr.shape})." ) - # Formula from Fan 2017 (equation 2) - # Determine unique original post-labeling delays (ignoring slice timing shifts) - unique_first_voxel_plds, unique_pld_idx = np.unique(first_voxel_plds, return_index=True) - unique_plds = plds[:, unique_pld_idx] # S x unique PLDs - n_unique_plds = unique_pld_idx.size - - # tau should be a 1D array, with one volume per unique PLD - if tau.size > 1: - if tau.size != n_plds: - raise ValueError( - f"Number of LabelingDurations ({tau.size}) != " - f"number of PostLabelingDelays ({n_plds})" - ) - - tau = tau[unique_pld_idx] - else: - tau = np.full(n_unique_plds, tau) - - mean_deltam_by_pld = np.zeros((n_voxels, n_unique_plds)) - for i_pld, first_voxel_pld in enumerate(unique_first_voxel_plds): - pld_idx = first_voxel_plds == first_voxel_pld - mean_deltam_by_pld[:, i_pld] = np.mean(deltam_arr[:, pld_idx], axis=1) - - # Estimate ATT for each voxel - att_arr = estimate_att_pcasl( - deltam_arr=mean_deltam_by_pld, - plds=unique_plds, - lds=tau, - t1blood=t1blood, - t1tissue=t1tissue, - ) - - # Start calculating CBF - num_factor = unit_conversion * partition_coefficient - denom_factor = 2 * labeleff * scaled_m0data * t1blood - - # Loop over PLDs and calculate CBF for each, accounting for ATT. - cbf_by_pld = np.zeros((n_voxels, n_unique_plds)) - for i_pld in range(n_unique_plds): - pld_by_voxel = unique_plds[:, i_pld] - tau_for_pld = tau[i_pld] - - pld_num_factor = num_factor * mean_deltam_by_pld[:, i_pld] * np.exp(att_arr / t1blood) - pld_denom_factor = denom_factor * ( - np.exp(-np.maximum(pld_by_voxel - att_arr, 0)) - - np.exp(-np.maximum(tau_for_pld + pld_by_voxel - att_arr, 0)) + if scaled_m0data.size != n_voxels: + raise ValueError( + f"Number of voxels in M0 data ({scaled_m0data.size}) does not match " + f"number of delta-M voxels ({n_voxels})." ) - cbf_by_pld[:, i_pld] = pld_num_factor / pld_denom_factor - - # Average CBF across PLDs, but only include PLDs where PLD + tau > ATT for that voxel, - # per Juttukonda 2021 (section 2.6). - cbf = np.zeros(n_voxels) # mean CBF - for i_voxel in range(n_voxels): - plds_voxel = unique_plds[i_voxel, :] - cbf_by_pld_voxel = cbf_by_pld[i_voxel, :] - att_voxel = att_arr[i_voxel] - cbf[i_voxel] = np.mean(cbf_by_pld_voxel[(plds_voxel + tau) > att_voxel]) - return att_arr, cbf - - -def estimate_t1(metadata): - """Estimate the relaxation rates of blood and gray matter based on magnetic field strength. - - t1blood is set based on the scanner's field strength, - according to :footcite:t:`zhang2013vivo,alsop_recommended_2015`. - If recommended values from these publications cannot be used - (i.e., if the field strength isn't 1.5T, 3T, 7T), - then the formula from :footcite:t:`zhang2013vivo` will be applied. - - t1tissue is set based on the scanner's field strength as well, - according to :footcite:t:`wright2008water`. - At the moment, field strengths other than 1.5T, 3T, and 7T are not supported and will - raise an exception. - - Parameters - ---------- - metadata : :obj:`dict` - Dictionary of metadata from the ASL file. + if is_casl: + assert ti1 is None + assert tau is not None + if isinstance(tau, float): + tau = np.full((n_volumes,), tau) + else: + tau = np.ndarray(tau) - Returns - ------- - t1blood : :obj:`float` - Estimated relaxation rate of blood based on magnetic field strength. - t1tissue : :obj:`float` - Estimated relaxation rate of gray matter based on magnetic field strength. + if tau.size != n_volumes: + raise ValueError( + f"Number of values in tau ({tau.shape}) does not match number of " + f"delta-M volumes ({n_volumes})." + ) + else: + assert ti1 is not None + assert tau is None - References - ---------- - .. footbibliography:: - """ - T1BLOOD_DICT = { - 1.5: 1.35, - 3: 1.65, - 7: 2.087, - } - t1blood = T1BLOOD_DICT.get(metadata["MagneticFieldStrength"]) - if not t1blood: - config.loggers.interface.warning( - f"T1blood cannot be inferred for {metadata['MagneticFieldStrength']}T data. " - "Defaulting to formula from Zhang et al. (2013)." - ) - t1blood = (110 * metadata["MagneticFieldStrength"] + 1316) / 1000 + m0_a = scaled_m0data / partition_coefficient - # TODO: Supplement with formula for other field strengths - T1TISSUE_DICT = { - 1.5: 1.197, - 3: 1.607, - 7: 1.939, - } - t1tissue = T1TISSUE_DICT.get(metadata["MagneticFieldStrength"]) - if not t1tissue: - raise ValueError( - f"T1tissue cannot be inferred for {metadata['MagneticFieldStrength']}T data." + cbf = np.zeros(n_voxels) + att = np.zeros(n_voxels) + abat = np.zeros(n_voxels) + abv = np.zeros(n_voxels) + fail_count = 0 + for i_voxel in range(n_voxels): + deltam_voxel = deltam_arr[i_voxel, :] + plds_voxel = plds[i_voxel, :] + m0_a_voxel = m0_a[i_voxel] + + # TODO: Figure out what m0_b should be. + # The independent variables used to estimate cbf, etc. are either floats or arrays, + # but curve_fit needs them all to be the same size/shape. + xdata = np.zeros((n_volumes, 6)) + xdata[0, 0] = labeleff + xdata[0, 1] = t1blood + xdata[0, 2] = m0_a_voxel + xdata[0, 3] = m0_a_voxel + xdata[:, 4] = plds_voxel + + if is_casl: + xdata[:, 5] = tau + func = calculate_deltam_pcasl + else: + xdata[0, 5] = ti1 + func = calculate_deltam_pasl + + try: + popt = scipy.optimize.curve_fit( + func, + xdata=xdata, + ydata=deltam_voxel, + # Initial estimates for DVs (CBF, ATT, aBAT, aBV) + # Values provided by Manuel Taso + p0=(60, 1.2, 1, 0.02), + # Lower and upper bounds for DVs + # Upper bounds provided by Manuel Taso + bounds=( + (0, 0, 0, 0), + # Manuel Taso recommended 5, 5, 0.2 for ATT, aBAT, and aBV, + # but our test data maxed out around 12 when left unbounded. + (300, 5, 5, 0.1), + ), + )[0] + cbf[i_voxel] = popt[0] + att[i_voxel] = popt[1] + abat[i_voxel] = popt[2] + abv[i_voxel] = popt[3] + + except (RuntimeError, ValueError): + # If curve fit fails to converge, set voxel's value to NaN + cbf[i_voxel] = np.nan + att[i_voxel] = np.nan + abat[i_voxel] = np.nan + abv[i_voxel] = np.nan + + fail_count += 1 + + if fail_count: + fail_percent = 100 * fail_count / n_voxels + config.loggers.workflow.warning( + f"General kinetic model fit failed on " + f"{fail_count}/{n_voxels} ({fail_percent:.2f}%) voxel(s)." ) - return t1blood, t1tissue + return cbf, att, abat, abv diff --git a/aslprep/utils/plotting.py b/aslprep/utils/plotting.py index b24a07a29..643f15512 100644 --- a/aslprep/utils/plotting.py +++ b/aslprep/utils/plotting.py @@ -25,13 +25,14 @@ class CBFPlot(object): This plot restricts CBF values to -20 (if there are negative values) or 0 (if not) to 100. """ - __slots__ = ["cbf", "ref_vol", "label", "outfile", "vmax"] + __slots__ = ["cbf", "ref_vol", "label", "outfile", "vmin", "vmax"] - def __init__(self, cbf, ref_vol, label, outfile, vmax): + def __init__(self, cbf, ref_vol, label, outfile, vmin, vmax): self.cbf = cbf self.ref_vol = ref_vol self.label = label self.outfile = outfile + self.vmin = vmin self.vmax = vmax def plot(self): @@ -41,8 +42,8 @@ def plot(self): """ cbf_img = nb.load(self.cbf) cbf_data = cbf_img.get_fdata() - cbf_data[cbf_data < -20] = -20 - cbf_data[cbf_data > 100] = 100 + cbf_data[cbf_data < self.vmin] = self.vmin + cbf_data[cbf_data > self.vmax] = self.vmax cbf_img = nb.Nifti1Image(cbf_data, affine=cbf_img.affine, header=cbf_img.header) statfile = plot_stat_map( cbf=cbf_img, @@ -89,7 +90,7 @@ def plot_stat_map( display_mode=mode, cut_coords=cuts[mode], vmax=vmax, - threshold=0.02, + threshold=0.00001, draw_cross=False, colorbar=True, symmetric_cbar=False, diff --git a/aslprep/workflows/asl/base.py b/aslprep/workflows/asl/base.py index 432a3a9d7..e3926faaf 100644 --- a/aslprep/workflows/asl/base.py +++ b/aslprep/workflows/asl/base.py @@ -6,6 +6,7 @@ import numpy as np from fmriprep.workflows.bold.apply import init_bold_volumetric_resample_wf from nipype.interfaces import utility as niu +from nipype.interfaces.afni.utils import Resample from nipype.pipeline import engine as pe from aslprep import config @@ -232,7 +233,11 @@ def init_asl_wf( cbf_4d_derivs = [] if is_multi_pld: - att_derivs += ["att"] + att_derivs += [ + "att", + "abat", + "abv", + ] else: cbf_4d_derivs += ["cbf_ts"] @@ -317,10 +322,25 @@ def init_asl_wf( inputnode.inputs.aslcontext = run_data["aslcontext"] inputnode.inputs.m0scan_metadata = run_data["m0scan_metadata"] + m0scan_buffer = pe.Node( + niu.IdentityInterface(fields=["m0scan"]), + name="m0scan_buffer", + ) + # Resample M0 scan to ASL FOV and resolution + if metadata["M0Type"] == "Separate": + resample = pe.Node( + Resample( + in_file=run_data["m0scan"], + master=asl_file, + outputtype="NIFTI", + ), + name="resample_m0scan_to_asl", + ) + workflow.connect([(resample, m0scan_buffer, [("out_file", "m0scan")])]) + # Perform minimal preprocessing of the ASL data, including HMC and SDC asl_fit_wf = init_asl_fit_wf( asl_file=asl_file, - m0scan=run_data["m0scan"], use_ge=use_ge, precomputed=precomputed, fieldmap_id=fieldmap_id, @@ -344,6 +364,7 @@ def init_asl_wf( ("fmap_id", "inputnode.fmap_id"), ("sdc_method", "inputnode.sdc_method"), ]), + (m0scan_buffer, asl_fit_wf, [("m0scan", "inputnode.m0scan")]), ]) # fmt:skip # Resample to aslref space. @@ -351,7 +372,7 @@ def init_asl_wf( # because ASLPrep's main output is CBF and we need aslref-space data to calculate CBF. asl_native_wf = init_asl_native_wf( asl_file=asl_file, - m0scan=run_data["m0scan"], + m0type=metadata["M0Type"], fieldmap_id=fieldmap_id, omp_nthreads=omp_nthreads, name="asl_native_wf", @@ -364,6 +385,7 @@ def init_asl_wf( ("fmap_coeff", "inputnode.fmap_coeff"), ("fmap_id", "inputnode.fmap_id"), ]), + (m0scan_buffer, asl_native_wf, [("m0scan", "inputnode.m0scan")]), (asl_fit_wf, asl_native_wf, [ ("outputnode.coreg_aslref", "inputnode.aslref"), ("outputnode.asl_mask", "inputnode.asl_mask"), @@ -488,6 +510,7 @@ def init_asl_wf( plot_timeseries=not (is_multi_pld or use_ge or (config.workflow.level == "resampling")), scorescrub=scorescrub, basil=basil, + is_multi_pld=is_multi_pld, name="cbf_reporting_wf", ) workflow.connect([ diff --git a/aslprep/workflows/asl/cbf.py b/aslprep/workflows/asl/cbf.py index 648ddcaaf..e916848ac 100644 --- a/aslprep/workflows/asl/cbf.py +++ b/aslprep/workflows/asl/cbf.py @@ -109,7 +109,7 @@ def init_cbf_wf( workflow = Workflow(name=name) workflow.__desc__ = """ -### Cerebral blood flow computation and denoising +Cerebral blood flow computation and denoising """ @@ -147,70 +147,48 @@ def init_cbf_wf( if processing_target == "cbf": workflow.__desc__ += """\ -*ASLPrep* loaded pre-calculated cerebral blood flow (CBF) data from the ASL file. +: *ASLPrep* loaded pre-calculated cerebral blood flow (CBF) data from the ASL file. """ - elif is_casl: - if is_multi_pld: - workflow.__desc__ += f"""\ -*ASLPrep* calculated cerebral blood flow (CBF) from the multi-delay -{metadata['ArterialSpinLabelingType']} data using the following method. + elif is_multi_pld: + workflow.__desc__ += f"""\ +: *ASLPrep* calculated cerebral blood flow (CBF) from the multi-delay +{metadata['ArterialSpinLabelingType']} data using a two-compartment general kinetic model (GKM) +[@buxton1998general], as recommended and extended in @woods2023recommendations. +This model contains separate terms for tissue and macrovascular delta-M signals. -First, delta-M values were averaged over time for each post-labeling delay (PLD). {m0_str} +The voxel-wise M0 values were used as both M0a and M0b in the GKM. -Next, arterial transit time (ATT) was estimated on a voxel-wise basis according to -@dai2012reduced. - -CBF was then calculated for each delay using the mean delta-M values and the estimated ATT, -according to the formula from @fan2017long. - -CBF was then averaged over delays according to @juttukonda2021characterizing, -in which an unweighted average is calculated for each voxel across all delays in which -PLD + labeling duration > ATT. +CBF, arterial transit time (ATT), arterial bolus arrival time (aBAT), +and arterial blood volume (aBV) were estimated using a nonlinear model fit with SciPy's +``curve_fit``. """ - else: - # Single-delay (P)CASL data - workflow.__desc__ += f"""\ -*ASLPrep* calculated cerebral blood flow (CBF) from the single-delay + + elif is_casl: + # Single-delay (P)CASL data + workflow.__desc__ += f"""\ +: *ASLPrep* calculated cerebral blood flow (CBF) from the single-delay {metadata['ArterialSpinLabelingType']} using a single-compartment general kinetic model [@buxton1998general]. {m0_str} """ else: - bcut = metadata.get("BolusCutOffTechnique") - - if is_multi_pld: - raise ValueError( - "Multi-delay data are not supported for PASL sequences at the moment." - ) - # Single-delay PASL data, with different bolus cut-off techniques - if bcut == "QUIPSS": - workflow.__desc__ += f"""\ -*ASLPrep* calculated cerebral blood flow (CBF) from the single-delay PASL -using a single-compartment general kinetic model [@buxton1998general] -using the QUIPSS modification, as described in @wong1998quantitative. -{m0_str} -""" - elif bcut == "QUIPSSII": - workflow.__desc__ += f"""\ -*ASLPrep* calculated cerebral blood flow (CBF) from the single-delay PASL -using a single-compartment general kinetic model [@buxton1998general] -using the QUIPSS II modification, as described in @alsop_recommended_2015. -{m0_str} -""" - elif bcut == "Q2TIPS": - workflow.__desc__ += f"""\ -*ASLPrep* calculated cerebral blood flow (CBF) from the single-delay PASL + bcut = metadata.get("BolusCutOffTechnique") + singlepld_pasl_strs = { + "QUIPSS": "@wong1998quantitative", + "QUIPSSII": "@alsop_recommended_2015", + "Q2TIPS": "@noguchi2015technical", + } + + workflow.__desc__ += f"""\ +: *ASLPrep* calculated cerebral blood flow (CBF) from the single-delay PASL using a single-compartment general kinetic model [@buxton1998general] -using the Q2TIPS modification, as described in @noguchi2015technical. +using the {bcut} modification, as described in {singlepld_pasl_strs[bcut]}. {m0_str} """ - else: - # No bolus cutoff delay technique - raise ValueError("PASL without a bolus cut-off technique is not supported in ASLPrep.") if "SliceTiming" in metadata: workflow.__desc__ += ( @@ -240,6 +218,8 @@ def init_cbf_wf( "mean_cbf", "cbf_ts", # Only calculated for single-delay data "att", # Only calculated for multi-delay data + "abat", # Only calculated for multi-delay data + "abv", # Only calculated for multi-delay data "plds", # SCORE/SCRUB outputs "cbf_ts_score", @@ -425,6 +405,8 @@ def _getfiledir(file): ("mean_cbf", "mean_cbf"), ("att", "att"), ("plds", "plds"), + ("abat", "abat"), + ("abv", "abv"), ]), ]) # fmt:skip diff --git a/aslprep/workflows/asl/confounds.py b/aslprep/workflows/asl/confounds.py index 45f907922..d98a32681 100644 --- a/aslprep/workflows/asl/confounds.py +++ b/aslprep/workflows/asl/confounds.py @@ -105,7 +105,7 @@ def init_asl_confounds_wf( workflow.__desc__ = """\ Several confounding timeseries were calculated, including both framewise displacement (FD) and DVARS. FD and DVARS are calculated using the implementations in in *Nipype* -(following the definition by [@power_fd_dvars]) for each ASL run. ASLPrep summarizes +(following the definition by [@power_fd_dvars]) for each ASL run. *ASLPrep* summarizes in-scanner motion as the mean framewise displacement and relative root-mean square displacement. """ inputnode = pe.Node( diff --git a/aslprep/workflows/asl/fit.py b/aslprep/workflows/asl/fit.py index 1f891ef2e..4832f6d22 100644 --- a/aslprep/workflows/asl/fit.py +++ b/aslprep/workflows/asl/fit.py @@ -86,7 +86,6 @@ def get_sbrefs( def init_asl_fit_wf( *, asl_file: str, - m0scan: ty.Union[str, None], use_ge: bool, precomputed: dict = {}, fieldmap_id: ty.Optional[str] = None, @@ -111,7 +110,6 @@ def init_asl_fit_wf( ) wf = init_asl_fit_wf( asl_file=str(asl_file), - m0scan=None, use_ge=False, ) @@ -119,8 +117,6 @@ def init_asl_fit_wf( ---------- asl_file Path to ASL NIfTI file. - m0scan - Path to M0 NIfTI file. use_ge If True, the M0 scan (when available) will be prioritized as the volume type for the reference image, as GE deltam volumes exhibit extreme background noise. @@ -134,6 +130,7 @@ def init_asl_fit_wf( ------ asl_file ASL series NIfTI file + m0scan t1w_preproc Bias-corrected structural template image t1w_mask @@ -243,6 +240,7 @@ def init_asl_fit_wf( fields=[ "asl_file", "aslcontext", + "m0scan", # only defined if metadata['M0Type'] is Separate # Fieldmap registration "fmap", "fmap_ref", @@ -365,10 +363,14 @@ def init_asl_fit_wf( m0scan=(metadata["M0Type"] == "Separate"), use_ge=use_ge, ) - hmc_aslref_wf.inputs.inputnode.m0scan = m0scan hmc_aslref_wf.inputs.inputnode.dummy_scans = config.workflow.dummy_scans - workflow.connect([(inputnode, hmc_aslref_wf, [("aslcontext", "inputnode.aslcontext")])]) + workflow.connect([ + (inputnode, hmc_aslref_wf, [ + ("aslcontext", "inputnode.aslcontext"), + ("m0scan", "inputnode.m0scan"), + ]), + ]) # fmt:skip ds_hmc_aslref_wf = init_ds_aslref_wf( bids_root=layout.root, @@ -632,7 +634,7 @@ def init_asl_fit_wf( def init_asl_native_wf( *, asl_file: str, - m0scan: ty.Optional[str] = None, + m0type: str, fieldmap_id: ty.Optional[str] = None, omp_nthreads: int = 1, name: str = "asl_native_wf", @@ -662,7 +664,8 @@ def init_asl_native_wf( ---------- asl_file List of paths to NIfTI files. - m0scan + m0type + Type of M0 scan. fieldmap_id ID of the fieldmap to use to correct this ASL series. If :obj:`None`, no correction will be applied. @@ -815,6 +818,10 @@ def init_asl_native_wf( ]), ]) # fmt:skip + workflow.__desc__ = """\ +The ASL time series was resampled to ASL-reference space. +""" + # Resample ASL to aslref aslref_asl = pe.Node( ResampleSeries(jacobian="fmap-jacobian" not in config.workflow.ignore), @@ -855,9 +862,14 @@ def init_asl_native_wf( (aslref_asl, outputnode, [("out_file", "asl_native")]), ]) # fmt:skip - if m0scan: + if m0type == "Separate": from niworkflows import data as nw_data + workflow.__desc__ += """\ +The separate calibration volume (M0) was resampled to the same resolution and orientation as +the ASL time series. +""" + # Resample separate M0 file to aslref # No HMC identity_xfm = nw_data.load("itkIdentityTransform.txt") @@ -872,7 +884,10 @@ def init_asl_native_wf( ) workflow.connect([ - (inputnode, aslref_m0scan, [("aslref", "ref_file")]), + (inputnode, aslref_m0scan, [ + ("aslref", "ref_file"), + ("m0scan", "in_file"), + ]), (aslbuffer, aslref_m0scan, [ ("ro_time", "ro_time"), ("pe_dir", "pe_dir"), diff --git a/aslprep/workflows/asl/hmc.py b/aslprep/workflows/asl/hmc.py index 3f2390d40..a2a4bd695 100644 --- a/aslprep/workflows/asl/hmc.py +++ b/aslprep/workflows/asl/hmc.py @@ -90,7 +90,7 @@ def init_asl_hmc_wf( in order to account for intensity differences between different contrasts, which, when motion corrected together, can conflate intensity differences with head motions [@wang2008empirical]. -Next, ASLPrep concatenated the motion parameters across volume types and +Next, *ASLPrep* concatenated the motion parameters across volume types and re-calculated relative root mean-squared deviation. """ diff --git a/aslprep/workflows/asl/outputs.py b/aslprep/workflows/asl/outputs.py index 0be25749d..cd5978171 100644 --- a/aslprep/workflows/asl/outputs.py +++ b/aslprep/workflows/asl/outputs.py @@ -34,6 +34,15 @@ }, "att": { "suffix": "att", + "Units": "s", + }, + "abat": { + "suffix": "abat", + "Units": "s", + }, + "abv": { + "suffix": "abv", + "Units": "fraction", }, # SCORE/SCRUB outputs "cbf_ts_score": { @@ -64,6 +73,7 @@ "att_basil": { "desc": "basil", "suffix": "att", + "Units": "s", }, } @@ -566,10 +576,6 @@ def init_ds_asl_native_wf( workflow.connect([(inputnode, ds_cbf, [(cbf_name, "in_file")])]) for att_name in att: - # TODO: Add EstimationReference and EstimationAlgorithm - att_meta = { - "Units": "s", - } fields = BASE_INPUT_FIELDS[att_name] ds_att = pe.Node( @@ -578,7 +584,6 @@ def init_ds_asl_native_wf( compress=True, dismiss_entities=("echo",), **fields, - **att_meta, ), name=f"ds_{att_name}", run_without_submitting=True, @@ -759,9 +764,6 @@ def init_ds_volumes_wf( for att_name in att: # TODO: Add EstimationReference and EstimationAlgorithm - att_meta = { - "Units": "s", - } fields = BASE_INPUT_FIELDS[att_name] resample_att = pe.Node( @@ -782,7 +784,6 @@ def init_ds_volumes_wf( compress=True, dismiss_entities=("echo",), **fields, - **att_meta, ), name=f"ds_{att_name}", run_without_submitting=True, diff --git a/aslprep/workflows/asl/plotting.py b/aslprep/workflows/asl/plotting.py index c622edf59..e41824344 100644 --- a/aslprep/workflows/asl/plotting.py +++ b/aslprep/workflows/asl/plotting.py @@ -20,6 +20,7 @@ def init_cbf_reporting_wf( plot_timeseries=True, scorescrub=False, basil=False, + is_multi_pld=False, name="cbf_reporting_wf", ): """Generate CBF reports. @@ -36,6 +37,7 @@ def init_cbf_reporting_wf( "RepetitionTime": 4, "RepetitionTimePreparation": 4, }, + is_multi_pld=True, ) """ from niworkflows.interfaces.images import SignalExtraction @@ -63,6 +65,8 @@ def init_cbf_reporting_wf( "cifti_cbf_ts", # Multi-delay outputs "att", + "abat", + "abv", # SCORE/SCRUB outputs "cbf_ts_score", # unused "mean_cbf_score", @@ -176,7 +180,15 @@ def init_cbf_reporting_wf( (cbf_confounds, carpetplot_wf, [("confounds_file", "inputnode.confounds_file")]), ]) # fmt:skip - cbf_summary = pe.Node(CBFSummaryPlot(label="cbf", vmax=100), name="cbf_summary", mem_gb=1) + cbf_summary = pe.Node( + CBFSummaryPlot( + label="cbf", + vmin=-20, + vmax=100, + ), + name="cbf_summary", + mem_gb=1, + ) workflow.connect([ (inputnode, cbf_summary, [ ("mean_cbf", "cbf"), @@ -193,11 +205,11 @@ def init_cbf_reporting_wf( workflow.connect([(cbf_summary, ds_report_cbf, [("out_file", "in_file")])]) cbf_by_tt_plot = pe.Node( - CBFByTissueTypePlot(), + CBFByTissueTypePlot(img_type="cbf"), name="cbf_by_tt_plot", ) workflow.connect([ - (inputnode, cbf_by_tt_plot, [("mean_cbf", "cbf")]), + (inputnode, cbf_by_tt_plot, [("mean_cbf", "in_file")]), (warp_t1w_dseg_to_aslref, cbf_by_tt_plot, [("output_image", "seg_file")]), ]) # fmt:skip @@ -214,9 +226,69 @@ def init_cbf_reporting_wf( ) workflow.connect([(cbf_by_tt_plot, ds_report_cbf_by_tt, [("out_file", "in_file")])]) + if is_multi_pld: + # Limits for the different figures. + # Make sure these match the hardcoded limits in the model-fitting function. + lims = { + "att": (0, 5), + "abat": (0, 5), + "abv": (0, 0.1), + } + for img_type in ["att", "abat", "abv"]: + img_summary = pe.Node( + CBFSummaryPlot( + label=img_type, + vmin=lims[img_type][0], + vmax=lims[img_type][1], + ), + name=f"{img_type}_summary", + mem_gb=1, + ) + workflow.connect([ + (inputnode, img_summary, [ + (img_type, "cbf"), + ("aslref", "ref_vol"), + ]), + ]) # fmt:skip + + ds_report_img = pe.Node( + DerivativesDataSink( + datatype="figures", + desc=img_type, + suffix="cbf", + keep_dtype=True, + ), + name=f"ds_report_{img_type}", + run_without_submitting=True, + mem_gb=config.DEFAULT_MEMORY_MIN_GB, + ) + workflow.connect([(img_summary, ds_report_img, [("out_file", "in_file")])]) + + img_by_tt_plot = pe.Node( + CBFByTissueTypePlot(img_type=img_type), + name=f"{img_type}_by_tt_plot", + ) + workflow.connect([ + (inputnode, img_by_tt_plot, [(img_type, "in_file")]), + (warp_t1w_dseg_to_aslref, img_by_tt_plot, [("output_image", "seg_file")]), + ]) # fmt:skip + + ds_report_img_by_tt = pe.Node( + DerivativesDataSink( + datatype="figures", + desc=f"{img_type}ByTissueType", + suffix=img_type, + keep_dtype=True, + ), + name=f"ds_report_{img_type}_by_tt", + run_without_submitting=True, + mem_gb=config.DEFAULT_MEMORY_MIN_GB, + ) + workflow.connect([(img_by_tt_plot, ds_report_img_by_tt, [("out_file", "in_file")])]) + if scorescrub: score_summary = pe.Node( - CBFSummaryPlot(label="score", vmax=100), + CBFSummaryPlot(label="score", vmin=-20, vmax=100), name="score_summary", mem_gb=1, ) @@ -236,11 +308,11 @@ def init_cbf_reporting_wf( workflow.connect([(score_summary, ds_report_score, [("out_file", "in_file")])]) score_by_tt_plot = pe.Node( - CBFByTissueTypePlot(), + CBFByTissueTypePlot(img_type="cbf"), name="score_by_tt_plot", ) workflow.connect([ - (inputnode, score_by_tt_plot, [("mean_cbf_score", "cbf")]), + (inputnode, score_by_tt_plot, [("mean_cbf_score", "in_file")]), (warp_t1w_dseg_to_aslref, score_by_tt_plot, [("output_image", "seg_file")]), ]) # fmt:skip @@ -258,7 +330,7 @@ def init_cbf_reporting_wf( workflow.connect([(score_by_tt_plot, ds_report_score_by_tt, [("out_file", "in_file")])]) scrub_summary = pe.Node( - CBFSummaryPlot(label="scrub", vmax=100), + CBFSummaryPlot(label="scrub", vmin=-20, vmax=100), name="scrub_summary", mem_gb=1, ) @@ -278,11 +350,11 @@ def init_cbf_reporting_wf( workflow.connect([(scrub_summary, ds_report_scrub, [("out_file", "in_file")])]) scrub_by_tt_plot = pe.Node( - CBFByTissueTypePlot(), + CBFByTissueTypePlot(img_type="cbf"), name="scrub_by_tt_plot", ) workflow.connect([ - (inputnode, scrub_by_tt_plot, [("mean_cbf_scrub", "cbf")]), + (inputnode, scrub_by_tt_plot, [("mean_cbf_scrub", "in_file")]), (warp_t1w_dseg_to_aslref, scrub_by_tt_plot, [("output_image", "seg_file")]), ]) # fmt:skip @@ -301,7 +373,7 @@ def init_cbf_reporting_wf( if basil: basil_summary = pe.Node( - CBFSummaryPlot(label="basil", vmax=100), + CBFSummaryPlot(label="basil", vmin=0, vmax=100), name="basil_summary", mem_gb=1, ) @@ -321,11 +393,11 @@ def init_cbf_reporting_wf( workflow.connect([(basil_summary, ds_report_basil, [("out_file", "in_file")])]) basil_by_tt_plot = pe.Node( - CBFByTissueTypePlot(), + CBFByTissueTypePlot(img_type="cbf"), name="basil_by_tt_plot", ) workflow.connect([ - (inputnode, basil_by_tt_plot, [("mean_cbf_basil", "cbf")]), + (inputnode, basil_by_tt_plot, [("mean_cbf_basil", "in_file")]), (warp_t1w_dseg_to_aslref, basil_by_tt_plot, [("output_image", "seg_file")]), ]) # fmt:skip @@ -343,7 +415,7 @@ def init_cbf_reporting_wf( workflow.connect([(basil_by_tt_plot, ds_report_basil_by_tt, [("out_file", "in_file")])]) pvc_summary = pe.Node( - CBFSummaryPlot(label="pvc", vmax=120), + CBFSummaryPlot(label="pvc", vmin=0, vmax=120), name="pvc_summary", mem_gb=1, ) @@ -363,11 +435,11 @@ def init_cbf_reporting_wf( workflow.connect([(pvc_summary, ds_report_pvc, [("out_file", "in_file")])]) pvc_by_tt_plot = pe.Node( - CBFByTissueTypePlot(), + CBFByTissueTypePlot(img_type="cbf"), name="pvc_by_tt_plot", ) workflow.connect([ - (inputnode, pvc_by_tt_plot, [("mean_cbf_gm_basil", "cbf")]), + (inputnode, pvc_by_tt_plot, [("mean_cbf_gm_basil", "in_file")]), (warp_t1w_dseg_to_aslref, pvc_by_tt_plot, [("output_image", "seg_file")]), ]) # fmt:skip diff --git a/aslprep/workflows/asl/resampling.py b/aslprep/workflows/asl/resampling.py index 7cbf1192d..d34ec9869 100644 --- a/aslprep/workflows/asl/resampling.py +++ b/aslprep/workflows/asl/resampling.py @@ -82,6 +82,9 @@ def init_asl_surf_wf( medial_surface_nan=False, metadata={}, output_dir=".", + cbf_3d=["mean_cbf"], + cbf_4d=[], + att=["att"], ) Parameters @@ -179,7 +182,8 @@ def select_target(subject_id, space): kwargs["dimension"] = 3 if cbf_deriv in att: - meta = {"Units": "s"} + # units provided in another variable + meta = {} else: meta = {"Units": "mL/100 g/min"} diff --git a/docs/installation.rst b/docs/installation.rst index 4481c754f..4637f0f89 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -6,8 +6,9 @@ Installation There are two ways to install *ASLPrep*: -* using `Container Technologies`_ (RECOMMENDED) -* within a `Manually Prepared Environment (Python 3.10+)`_, also known as *bare-metal installation* +* using `Container Technologies`_ (RECOMMENDED) +* within a `Manually Prepared Environment (Python 3.10)`_ + also known as *bare-metal installation* *ASLPrep* is not just a Python package; it also depends on a number of other neuroimaging libraries written in other languages. diff --git a/docs/links.rst b/docs/links.rst index c85764037..589149a7d 100644 --- a/docs/links.rst +++ b/docs/links.rst @@ -9,6 +9,8 @@ .. _`Docker installation`: https://docs.docker.com/install/ .. _FreeSurfer: https://surfer.nmr.mgh.harvard.edu/ .. _FSL: https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/ +.. _`HCP Pipelines`: https://humanconnectome.org/software/hcp-mr-pipelines/ .. _Nipype: http://nipype.readthedocs.io/en/latest/ +.. _QSIPrep: https://qsiprep.readthedocs.io .. _Singularity: https://github.com/singularityware/singularity .. _SPM: https://www.fil.ion.ucl.ac.uk/spm/software/spm12/ diff --git a/docs/outputs.rst b/docs/outputs.rst index 3b67cb767..26d45fe38 100644 --- a/docs/outputs.rst +++ b/docs/outputs.rst @@ -44,7 +44,7 @@ The outputs will be a `BIDS Derivatives`_ dataset of the form:: For each participant in the dataset, a directory of derivatives (``sub-