diff --git a/docs/source/conf.py b/docs/source/conf.py index 946a89ae7..88033a488 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -90,7 +90,7 @@ # # This is also used if you do content translation via gettext catalogs. # Usually you set "language" from the command line for these cases. -language = None +language = "en" # There are two options for replacing |today|: either, you set today to some # non-false value, then it is used: diff --git a/mriqc/cli/parser.py b/mriqc/cli/parser.py index 349dcfab6..907b7e43f 100644 --- a/mriqc/cli/parser.py +++ b/mriqc/cli/parser.py @@ -538,7 +538,7 @@ def parse_args(args=None, namespace=None): config.execution.participant_label, session_id=config.execution.session_id, task=config.execution.task_id, - group_echos=False, + group_echos=True, bids_filters=config.execution.bids_filters, queries={mod: DEFAULT_BIDS_QUERIES[mod] for mod in lc_modalities} ) @@ -575,7 +575,7 @@ def parse_args(args=None, namespace=None): # Estimate the biggest file size / leave 1GB if some file does not exist (datalad) with suppress(FileNotFoundError): config.workflow.biggest_file_gb = _get_biggest_file_size_gb( - [i for sublist in config.workflow.inputs.values() for i in sublist] + config.workflow.inputs.values() ) # set specifics for alternative populations @@ -593,11 +593,14 @@ def parse_args(args=None, namespace=None): def _get_biggest_file_size_gb(files): + """Identify the largest file size (allows multi-echo groups).""" + import os - max_size = 0 + sizes = [] for file in files: - size = os.path.getsize(file) / (1024**3) - if size > max_size: - max_size = size - return max_size + if isinstance(file, (list, tuple)): + sizes.append(_get_biggest_file_size_gb(file)) + else: + sizes.append(os.path.getsize(file)) + return max(sizes) / (1024**3) diff --git a/mriqc/data/bootstrap-func.yml b/mriqc/data/bootstrap-func.yml index ae07531f3..659c3b5e8 100644 --- a/mriqc/data/bootstrap-func.yml +++ b/mriqc/data/bootstrap-func.yml @@ -26,18 +26,14 @@ ########################################################################### packagename: mriqc -title: '{filename} :: Functional MRIQC report' +title: "{filename} :: MRIQC's BOLD fMRI report" sections: - name: Summary reportlets: - bids: {datatype: figures, desc: summary, extension: [.html]} -- name: Basic visual report +- name: Basic echo-wise reports + ordering: echo reportlets: - - bids: {datatype: figures, desc: zoomed} - caption: This panel shows a mosaic of the brain. This mosaic is the most suitable to - screen head-motion intensity inhomogeneities, global/local noise, signal leakage - (for example, from the eyeballs and across the phase-encoding axis), etc. - subtitle: Voxel-wise average of BOLD time-series, zoomed-in covering just the brain - bids: {datatype: figures, desc: stdev} subtitle: Standard deviation of signal through time caption: The voxel-wise standard deviation of the signal (variability along time). @@ -45,43 +41,12 @@ sections: subtitle: Carpetplot and nuisance signals caption: The so-called «carpetplot» may assist in assessing head-motion derived artifacts and respiation effects. -- name: Extended visual report +- name: Extended reports shared across echos reportlets: - - bids: {datatype: figures, desc: background} - caption: This panel shows a mosaic enhancing the background around the head. - Artifacts usually unveil themselves in the air surrounding the head, where no signal - sources are present. - subtitle: View of the background of the voxel-wise average of the BOLD timeseries - - bids: {datatype: figures, desc: mean} - subtitle: Average signal through time - caption: The average signal calculated across the last axis (time). - - bids: {datatype: figures, desc: airmask} - caption: The hat-mask calculated internally by MRIQC. Some metrics will use this - mask, for instance, to find out artifacts and estimate the spread of gaussian noise - added to the signal. This mask leaves out the air around the face to avoid measuring - noise sourcing from the eyeballs and their movement. - subtitle: '«Hat»-mask' - - bids: {datatype: figures, desc: noisefit} - caption: The noise fit internally estimated by MRIQC to calculate the QI1 index - proposed by Mortamet et al. (2009). - subtitle: Distribution of the noise within the hat mask - style: - max-width: 450px - - bids: {datatype: figures, desc: artifacts} - caption: Mask of artifactual intensities identified within the hat-mask. - subtitle: Artifactual intensities on the background - bids: {datatype: figures, desc: brainmask} caption: Brain mask as internally extracted by MRIQC. Defects on the brainmask could indicate problematic aspects of the image quality-wise. subtitle: Brain extraction performance - - bids: {datatype: figures, desc: head} - caption: A mask of the head calculated internally by MRIQC. - subtitle: Head mask - - bids: {datatype: figures, desc: segmentation} - caption: Brain tissue segmentation, as internally extracted by MRIQC. - Defects on this segmentation, as well as noisy tissue labels could - indicate problematic aspects of the image quality-wise. - subtitle: Brain tissue segmentation - bids: {datatype: figures, desc: norm} caption: This panel shows a quick-and-dirty nonlinear registration into the MNI152NLin2009cAsym template accessed with @@ -89,6 +54,22 @@ sections: subtitle: Spatial normalization of the anatomical image static: false +- name: Extended echo-wise reports + ordering: echo + reportlets: + - bids: {datatype: figures, desc: background} + caption: This panel shows a mosaic enhancing the background around the head. + Artifacts usually unveil themselves in the air surrounding the head, where no signal + sources are present. + subtitle: View of the background of the voxel-wise average of the BOLD timeseries + - bids: {datatype: figures, desc: mean} + subtitle: Average signal through time + caption: The average signal calculated across the last axis (time). + - bids: {datatype: figures, desc: zoomed} + caption: This panel shows a mosaic of the brain. This mosaic is the most suitable to + screen head-motion intensity inhomogeneities, global/local noise, signal leakage + (for example, from the eyeballs and across the phase-encoding axis), etc. + subtitle: Voxel-wise average of BOLD time-series, zoomed-in covering just the brain - name: About nested: true diff --git a/mriqc/data/testdata/group_bold.tsv b/mriqc/data/testdata/group_bold.tsv index 779bba82a..7506e9995 100644 --- a/mriqc/data/testdata/group_bold.tsv +++ b/mriqc/data/testdata/group_bold.tsv @@ -1,10 +1,10 @@ bids_name aor aqi dummy_trs dvars_nstd dvars_std dvars_vstd efc fber fd_mean fd_num fd_perc fwhm_avg fwhm_x fwhm_y fwhm_z gcor gsr_x gsr_y size_t size_x size_y size_z snr spacing_tr spacing_x spacing_y spacing_z summary_bg_k summary_bg_mad summary_bg_mean summary_bg_median summary_bg_n summary_bg_p05 summary_bg_p95 summary_bg_stdv summary_fg_k summary_fg_mad summary_fg_mean summary_fg_median summary_fg_n summary_fg_p05 summary_fg_p95 summary_fg_stdv tsnr -sub-ds205s03_task-functionallocalizer_run-01_bold 0.006175853658536585 0.00782700987804878 0 23.982385659135808 1.1029845019753086 1.0482825533333335 0.5579308153780409 899.73681640625 0.26118389619646404 32 39.02439024390244 2.043014166666667 1.92755 2.1932675 2.008225 0.0592534 -0.00898530799895525 0.04696183279156685 82 53 53 27 4.774209217164637 2.200000047683716 4.0 4.0 4.0 26.65223575628058 7.067350153923593 48.65173053536583 25.51134490966797 55113.0 18.871950149536133 173.51107788085938 62.96822585876589 0.598545770862756 139.55821373035255 769.8158763353686 765.2284240722656 20730.0 494.1331787109375 1041.6568603515625 160.27993998869655 56.32470703125 -sub-ds205s03_task-view_run-01_bold 0.0038563043478260866 0.010398965217391304 0 21.001368679560443 1.0080583680219781 0.9704656724175822 0.555268631792101 955.6907958984375 0.16520736608810047 18 19.565217391304348 1.9890683333333332 1.86912 2.1692175 1.9288675 0.0188429 -0.008540591225028038 0.04196126013994217 92 53 53 27 4.8924206716943655 2.200000047683716 4.0 4.0 4.0 24.78328036616209 7.284480140006654 48.01004600493635 25.337934494018555 55115.0 13.708492279052734 173.4195098876953 63.24374935408871 0.6193718723673731 135.63136294470496 785.625103328459 783.30322265625 20728.0 509.40557861328125 1057.23828125 160.10158982252509 56.46196937561035 -sub-ds205s03_task-view_run-02_bold 0.000538804347826087 0.00718998347826087 0 23.167324632417582 0.9661779338461531 0.9044185629670329 0.5566515976243607 906.2960815429688 0.22854716046970905 44 47.82608695652174 2.0270425000000003 1.898665 2.1941825 1.98828 0.350186 -0.008191297762095928 0.045305926352739334 92 53 53 27 4.819771892578718 2.200000047683716 4.0 4.0 4.0 24.552996566765586 7.135003381683008 48.0985481401069 25.58615016937256 55126.0 16.843839645385742 171.75970458984375 61.92351887941012 0.654934961442704 137.18513622307717 774.2342621452095 770.2647094726562 20717.0 498.0149230957031 1046.4920654296875 159.8096624186355 44.19522476196289 -sub-ds205s07_task-functionallocalizer_run-01_bold 0.005997083333333333 0.0056276609722222225 0 21.73490091760564 1.1161214411267604 1.0305108769014089 0.5398076828636086 979.1392822265625 0.1550487344317312 15 20.833333333333332 1.7913891666666668 1.7012025 1.8870975 1.7858675 0.0356344 -0.006339103449136019 0.029900137335062027 72 53 53 27 5.051380236662852 2.200000047683716 4.0 4.0 4.0 23.998434132266535 6.359841724478726 40.97200099097582 22.96632194519043 56199.0 13.510886192321777 148.19041442871094 54.7943334817377 0.42502710232227514 134.22454681982717 722.5941279030118 718.6438293457031 19644.0 477.47076416015625 953.461181640625 142.26320406635728 68.8161506652832 -sub-ds205s07_task-view_run-01_bold 0.0037941304347826085 0.009004425434782607 0 21.771277794505497 1.0383257392307694 0.9654323329670328 0.5374363883770776 996.05810546875 0.1520993052375479 15 16.304347826086957 1.8006975 1.698995 1.8907725 1.812325 0.0408048 -0.005473143886774778 0.029460299760103226 92 53 53 27 5.072217446804501 2.200000047683716 4.0 4.0 4.0 27.23682449098333 6.229042844827164 39.84597218066951 22.68582057952881 56266.0 13.13893985748291 139.88406372070312 53.19358553170623 0.3440927682297894 134.06510192792 720.2195639126317 715.9733276367188 19577.0 476.44183349609375 949.2680053710938 141.1522768548468 57.291603088378906 -sub-ds205s07_task-view_run-02_bold 0.001657826086956522 0.008996689347826086 0 21.519903371868143 1.03221585 0.9588135309890113 0.5394046827226318 978.5082397460938 0.15451170465747024 18 19.565217391304348 1.8057741666666667 1.7069875 1.8966275 1.8137075 0.0433865 -0.006386041175574064 0.03140100836753845 92 53 53 27 5.056518385567245 2.200000047683716 4.0 4.0 4.0 29.670872739050793 6.462488049835068 40.8786782498388 22.991168975830078 56204.0 13.375846862792969 145.70977783203125 55.307788548732276 0.4335730121828525 135.24980824168193 723.5399965684157 719.189453125 19639.0 481.1661376953125 953.2984008789062 142.22654558293314 59.76730728149414 -sub-ds205s09_task-view_acq-LR_run-01_bold 0.002976986301369863 0.013025251780821916 0 21.405761824861116 1.0611945656944441 0.9949810694444441 0.513063203363028 1207.71533203125 0.3258302280827085 45 61.64383561643836 1.7912675 1.6990475 1.91971 1.755045 0.0362147 -0.005499124526977539 0.01900400035083294 73 53 53 27 5.854353348400624 2.200000047683716 4.0 4.0 4.0 25.199608221336092 5.712747254148226 39.895556611046395 23.25909996032715 58097.0 15.119007110595703 139.44947814941406 56.095377999258176 1.178096485114681 117.15050549333893 796.3205729290771 808.3048706054688 17746.0 526.5880737304688 995.4718627929688 138.0651368134319 58.34690284729004 -sub-ds205s09_task-view_acq-LR_run-02_bold 0.005996575342465754 0.015215327260273973 0 34.39419094722222 1.6192232063888892 1.6501877419444442 0.5136460861041277 1192.3255615234375 0.3345943377548133 47 64.38356164383562 1.7920308333333335 1.6996125 1.920615 1.755865 0.0313875 -0.005142268259078264 0.018914325162768364 73 53 53 27 5.851459215059082 2.200000047683716 4.0 4.0 4.0 25.07199243533304 5.957362424448292 40.199861866481136 23.448966026306152 58096.0 15.263423919677734 139.87855529785156 56.15791257667716 1.1742077271466247 117.28293886411426 797.5247600734042 809.6943969726562 17747.0 526.8685302734375 997.0792846679688 138.3708840384822 51.352577209472656 -sub-ds205s09_task-view_acq-RL_run-01_bold 0.02578232876712329 0.02322782465753425 0 30.880940159027773 1.2222362369444444 1.0929697275 0.5115397239733915 1229.69287109375 0.5699059645824405 53 72.6027397260274 1.8186083333333334 1.72236 1.9520325 1.7814325 0.0447734 -0.003605152480304241 0.017525115981698036 73 53 53 27 5.840571587502318 2.200000047683716 4.0 4.0 4.0 28.28470212917516 5.378516442052904 39.71076207043191 23.19632911682129 58165.0 17.226463317871094 134.9324951171875 56.23862202860845 1.1305372215596474 117.89153513001754 800.2769560314895 813.4251403808594 17678.0 528.4014892578125 1001.2677612304688 139.26755646464403 45.58633995056152 +sub-ds205s03_task-functionallocalizer_run-01_bold 0.0038873170731707324 0.007676269512195123 0 23.556566426419753 1.111688350987654 1.037835274938272 0.5480388108408939 946.6738891601562 0.26118389619646404 32 39.02439024390244 2.460305 2.3309525 2.749775 2.3001875 0.0680646 -0.009326859377324581 0.05734974890947342 82 53 53 27 4.452167473312103 2.200000047683716 4.0 4.0 4.0 32.7759659672557 7.720380258975168 49.68507569151857 24.975608825683594 55688.0 2.341463327407837 179.5365753173828 79.87659545341923 1.3897322764135298 138.87641036027622 767.1765602872674 768.451171875 20155.0 461.5243835449219 1042.731689453125 172.5973051795346 60.897049377672374 +sub-ds205s03_task-view_run-01_bold 0.004324782608695651 0.008262664456521738 0 21.120018382857136 1.0831709981318685 1.0370256236263742 0.5459090166346154 1002.2031860351562 0.16520736608810047 18 19.565217391304348 2.342805 2.234215 2.613075 2.181125 0.0235103 -0.0050179013051092625 0.04793250188231468 92 53 53 27 4.4829699247170405 2.200000047683716 4.0 4.0 4.0 29.71988008857921 7.139052847581331 46.66574967004908 24.7608699798584 55362.0 0.032608695328235626 164.9891357421875 71.66784333938988 1.596677545134333 133.64368600192023 779.8773281882851 783.8695678710938 20481.0 457.5434875488281 1054.61962890625 174.85070038268285 62.58223795099184 +sub-ds205s03_task-view_run-02_bold 0.0012682608695652171 0.007114478586956522 0 33.47404192791209 1.3371308202197805 0.9764774098901093 0.5494999624875053 955.73681640625 0.22854716046970905 44 47.82608695652174 2.4598116666666665 2.3444425 2.7437 2.2912925 0.359354 -0.009432817809283733 0.052887145429849625 92 53 53 27 4.419438774661141 2.200000047683716 4.0 4.0 4.0 29.03249120482654 7.7192010499810335 49.54353785439428 24.934783935546875 55415.0 0.20652174949645996 181.78260803222656 80.45599936028407 1.4009872052560288 139.07454009400536 768.1427012059421 770.8587036132812 20428.0 459.5326232910156 1046.8260498046875 174.42029972293017 45.2559957196936 +sub-ds205s07_task-functionallocalizer_run-01_bold 0.006487361111111111 0.006664996805555556 0 26.95528650943662 1.4143280252112678 1.3278985984507037 0.5345065899755163 1009.7532348632812 0.1550487344317312 15 20.833333333333332 2.2598908333333334 2.16087 2.40925 2.2095525 0.0386767 -0.005986093543469906 0.03482493385672569 72 53 53 27 4.598492679701605 2.200000047683716 4.0 4.0 4.0 35.90517987432788 6.32165175453679 43.5521938439157 22.56944465637207 56450.0 0.0 157.3333282470703 77.35003442698779 1.656221989802022 134.60795658544046 716.0395720200634 717.1805419921875 19393.0 440.0694580078125 951.6666870117188 155.95589707318132 73.35539675597101 +sub-ds205s07_task-view_run-01_bold 0.005007608695652174 0.010561140108695653 0 26.932940011208803 1.2896941338461543 1.168551503516483 0.5321868933927048 1023.4715576171875 0.1520993052375479 15 16.304347826086957 2.266245 2.184035 2.4189475 2.1957525 0.0478984 -0.006112692411988974 0.03277275711297989 92 53 53 27 4.595802249004055 2.200000047683716 4.0 4.0 4.0 33.63529478261993 6.2365988314754315 42.79862974194479 22.380435943603516 56588.0 0.3804347813129425 153.56521606445312 75.85331035278197 1.585165284213823 135.14239559316104 713.6245609062549 715.9891357421875 19255.0 437.3913269042969 947.5869750976562 155.78793525432897 63.332828449318185 +sub-ds205s07_task-view_run-02_bold 0.0027991304347826083 0.008415782065217392 0 27.396022419560428 1.345632817142857 1.2136871409890106 0.5333849667366366 1016.1077880859375 0.15451170465747024 18 19.565217391304348 2.271710833333333 2.1726275 2.433375 2.20913 0.0505144 -0.005830460228025913 0.03532366082072258 92 53 53 27 4.598664236595992 2.200000047683716 4.0 4.0 4.0 36.078683930448996 6.575017673595487 43.39332640278381 22.565217971801758 56589.0 1.22826087474823 155.56521606445312 75.49560044499137 1.6540993322199018 136.17376514793506 717.751513360383 719.2989501953125 19254.0 443.8260803222656 952.5978393554688 156.41069531029106 64.89047080953605 +sub-ds205s09_task-view_acq-LR_run-01_bold 0.003159041095890411 0.010028379999999998 0 20.396152721388898 1.090989740416667 1.0059027894444443 0.5048353890125871 1243.320556640625 0.3258302280827085 45 61.64383561643836 2.0852608333333333 2.0588275 2.2995125 1.8974425 0.0461463 -0.0019858605228364468 0.021903924643993378 73 53 53 27 5.292682271102816 2.200000047683716 4.0 4.0 4.0 33.98009379704732 5.949264906571304 39.07313696894552 22.956466674804688 58369.0 1.0936830043792725 132.37034606933594 62.99641049623399 2.7398312528413937 115.58686871143783 791.6591133225957 809.4622802734375 17474.0 483.6707458496094 994.6715698242188 152.93552051211557 66.51667785644531 +sub-ds205s09_task-view_acq-LR_run-02_bold 0.005322328767123288 0.012026137260273973 0 33.41256753611111 1.7174510266666663 1.7695175424999998 0.5052536160578726 1231.6690673828125 0.3345943377548133 47 64.38356164383562 2.0866866666666666 2.0591175 2.3078075 1.893135 0.0393612 -0.0014911222970113158 0.021568207070231438 73 53 53 27 5.293464753553317 2.200000047683716 4.0 4.0 4.0 33.28261057706237 6.199152583273201 39.20170865624716 23.102066040039062 58361.0 1.1750245094299316 131.959716796875 62.4980652302799 2.7186948010130694 115.65690863557867 792.8781950379683 810.7703247070312 17482.0 484.3922119140625 996.912109375 153.16001398016576 59.753475189208984 +sub-ds205s09_task-view_acq-RL_run-01_bold 0.015073150684931508 0.015836369041095893 0 34.205267455972226 1.4370176472222222 1.2820233763888884 0.5044353080190683 1274.420654296875 0.5699059645824405 53 72.6027397260274 2.08056 2.068795 2.288455 1.88443 0.0560156 -0.0013055995805189013 0.018262900412082672 73 53 53 27 5.229556365898414 2.200000047683716 4.0 4.0 4.0 31.57766555712935 5.840893620522815 38.89056185642749 22.814455032348633 58350.0 3.888171911239624 132.03933715820312 61.89358205002823 2.599912401497959 117.24294190485038 794.526159305144 814.4539184570312 17493.0 471.96441650390625 999.9562377929688 155.73608576273267 52.76657485961914 diff --git a/mriqc/data/tests/ds000005/sub-01/func/sub-01_task-mixedgamblestask_run-01_bold.json b/mriqc/data/tests/ds000005/sub-01/func/sub-01_task-mixedgamblestask_run-01_bold.json new file mode 100644 index 000000000..3a93d7ce1 --- /dev/null +++ b/mriqc/data/tests/ds000005/sub-01/func/sub-01_task-mixedgamblestask_run-01_bold.json @@ -0,0 +1,3 @@ +{ + "EchoTime": 0.030 +} \ No newline at end of file diff --git a/mriqc/data/tests/ds000005/sub-01/func/sub-01_task-mixedgamblestask_run-02_bold.json b/mriqc/data/tests/ds000005/sub-01/func/sub-01_task-mixedgamblestask_run-02_bold.json new file mode 120000 index 000000000..22bc9c6f2 --- /dev/null +++ b/mriqc/data/tests/ds000005/sub-01/func/sub-01_task-mixedgamblestask_run-02_bold.json @@ -0,0 +1 @@ +sub-01_task-mixedgamblestask_run-01_bold.json \ No newline at end of file diff --git a/mriqc/data/tests/ds000005/sub-01/func/sub-01_task-mixedgamblestask_run-03_bold.json b/mriqc/data/tests/ds000005/sub-01/func/sub-01_task-mixedgamblestask_run-03_bold.json new file mode 120000 index 000000000..22bc9c6f2 --- /dev/null +++ b/mriqc/data/tests/ds000005/sub-01/func/sub-01_task-mixedgamblestask_run-03_bold.json @@ -0,0 +1 @@ +sub-01_task-mixedgamblestask_run-01_bold.json \ No newline at end of file diff --git a/mriqc/interfaces/bids.py b/mriqc/interfaces/bids.py index eb48122ae..d4ae9de42 100644 --- a/mriqc/interfaces/bids.py +++ b/mriqc/interfaces/bids.py @@ -49,6 +49,7 @@ class IQMFileSinkInputSpec(DynamicTraitedSpec, BaseInterfaceInputSpec): rec_id = traits.Either(None, Str, usedefault=True) run_id = traits.Either(None, traits.Int, usedefault=True) dataset = Str(desc="dataset identifier") + dismiss_entities = traits.List(["part"], usedefault=True) metadata = traits.Dict() provenance = traits.Dict() @@ -111,6 +112,17 @@ def _gen_outfile(self): break in_file = str(path.relative_to(bids_root)) + if ( + isdefined(self.inputs.dismiss_entities) + and (dismiss := self.inputs.dismiss_entities) + ): + for entity in dismiss: + bids_chunks = [ + chunk for chunk in path.name.split("_") + if not chunk.startswith(f"{entity}-") + ] + path = path.parent / "_".join(bids_chunks) + # Build path and ensure directory exists bids_path = out_dir / in_file.replace("".join(Path(in_file).suffixes), ".json") bids_path.parent.mkdir(parents=True, exist_ok=True) diff --git a/mriqc/interfaces/functional.py b/mriqc/interfaces/functional.py index 38670745a..dd72af0f6 100644 --- a/mriqc/interfaces/functional.py +++ b/mriqc/interfaces/functional.py @@ -20,6 +20,7 @@ # # https://www.nipreps.org/community/licensing/ # +from __future__ import annotations from os import path as op import nibabel as nb @@ -30,9 +31,12 @@ from nipype.interfaces.base import ( BaseInterfaceInputSpec, File, + InputMultiObject, + isdefined, SimpleInterface, TraitedSpec, traits, + Undefined, ) @@ -267,6 +271,44 @@ def _run_interface(self, runtime): return runtime +class _SelectEchoInputSpec(BaseInterfaceInputSpec): + in_files = InputMultiObject(File(exists=True), mandatory=True, desc="input EPI file(s)") + metadata = InputMultiObject(traits.Dict(), desc="sidecar JSON files corresponding to in_files") + te_reference = traits.Float(0.030, usedefault=True, desc="reference SE-EPI echo time") + + +class _SelectEchoOutputSpec(TraitedSpec): + out_file = File(desc="selected echo") + echo_index = traits.Int(desc="index of the selected echo") + is_multiecho = traits.Bool(desc="whether it is a multiecho dataset") + + +class SelectEcho(SimpleInterface): + """ + Computes anatomical :abbr:`QC (Quality Control)` measures on the + structural image given as input + + """ + + input_spec = _SelectEchoInputSpec + output_spec = _SelectEchoOutputSpec + + def _run_interface(self, runtime): + ( + self._results["out_file"], + self._results["echo_index"], + ) = select_echo( + self.inputs.in_files, + te_echos=( + _get_echotime(self.inputs.metadata) if isdefined(self.inputs.metadata) + else None + ), + te_reference=self.inputs.te_reference, + ) + self._results["is_multiecho"] = self._results["echo_index"] != -1 + return runtime + + def find_peaks(data): t_z = [data[:, :, i, :].mean(axis=0).mean(axis=0) for i in range(data.shape[2])] return t_z @@ -292,3 +334,95 @@ def find_spikes(data, spike_thresh): def _robust_zscore(data): return (data - np.atleast_2d(np.median(data, axis=1)).T) / np.atleast_2d(data.std(axis=1)).T + + +def select_echo( + in_files: str | list[str], + te_echos: list[float | Undefined | None] | None = None, + te_reference: float = 0.030, +) -> str: + """ + Select the echo file with the closest echo time to the reference echo time. + + Used to grab the echo file when processing multi-echo data through workflows + that only accept a single file. + + Parameters + ---------- + in_files : :obj:`str` or :obj:`list` + A single filename or a list of filenames. + te_echos : :obj:`list` of :obj:`float` + List of echo times corresponding to each file. + If not a number (typically, a :obj:`~nipype.interfaces.base.Undefined`), + the function selects the second echo. + te_reference : float, optional + Reference echo time used to find the closest echo time. + + Returns + ------- + str + The selected echo file. + + Examples + -------- + >>> select_echo("single-echo.nii.gz") + ('single-echo.nii.gz', -1) + + >>> select_echo(["single-echo.nii.gz"]) + ('single-echo.nii.gz', -1) + + >>> select_echo( + ... [f"echo{n}.nii.gz" for n in range(1,7)], + ... ) + ('echo2.nii.gz', 1) + + >>> select_echo( + ... [f"echo{n}.nii.gz" for n in range(1,7)], + ... te_echos=[12.5, 28.5, 34.2, 45.0, 56.1, 68.4], + ... te_reference=33.1, + ... ) + ('echo3.nii.gz', 2) + + >>> select_echo( + ... [f"echo{n}.nii.gz" for n in range(1,7)], + ... te_echos=[12.5, 28.5, 34.2, 45.0, 56.1], + ... te_reference=33.1, + ... ) + ('echo2.nii.gz', 1) + + >>> select_echo( + ... [f"echo{n}.nii.gz" for n in range(1,7)], + ... te_echos=[12.5, 28.5, 34.2, 45.0, 56.1, None], + ... te_reference=33.1, + ... ) + ('echo2.nii.gz', 1) + + """ + if not isinstance(in_files, (list, tuple)): + return in_files, -1 + + if len(in_files) == 1: + return in_files[0], -1 + + import numpy as np + + n_echos = len(in_files) + if te_echos is not None and len(te_echos) == n_echos: + try: + index = np.argmin(np.abs(np.array(te_echos) - te_reference)) + return in_files[index], index + except TypeError: + pass + + return in_files[1], 1 + + +def _get_echotime(inlist): + if isinstance(inlist, list): + retval = [_get_echotime(el) for el in inlist] + return retval[0] if len(retval) == 1 else retval + + echo_time = inlist.get("EchoTime", None) if inlist else None + + if echo_time: + return float(echo_time) diff --git a/mriqc/reports/individual.py b/mriqc/reports/individual.py index 1e2350dde..8619f1ec2 100644 --- a/mriqc/reports/individual.py +++ b/mriqc/reports/individual.py @@ -47,11 +47,13 @@ def _single_report(in_file): from mriqc import config # Ensure it's a Path - in_file = Path(in_file) + in_file = Path(in_file if not isinstance(in_file, list) else in_file[0]) # Extract BIDS entities entities = config.execution.layout.get_file(in_file).get_entities() entities.pop("extension", None) + entities.pop("echo", None) + entities.pop("part", None) report_type = entities.pop("datatype", None) # Read output file: diff --git a/mriqc/utils/bids.py b/mriqc/utils/bids.py index 4ab142655..40bbfed8d 100644 --- a/mriqc/utils/bids.py +++ b/mriqc/utils/bids.py @@ -21,6 +21,8 @@ # https://www.nipreps.org/community/licensing/ # """PyBIDS tooling.""" +from __future__ import annotations + import json import os from pathlib import Path @@ -90,3 +92,99 @@ def write_derivative_description(bids_dir, deriv_dir): desc["License"] = orig_desc["License"] Path.write_text(deriv_dir / "dataset_description.json", json.dumps(desc, indent=4)) + + +def derive_bids_fname( + orig_path: str | Path, + entity: str | None = None, + newsuffix: str | None = None, + newpath: str | Path | None = None, + newext: str | None = None, + position: int = -1, + absolute: bool = True, +) -> Path | str: + """ + Derive a new file name from a BIDS-formatted path. + + Parameters + ---------- + orig_path : :obj:`str` or :obj:`os.pathlike` + A filename (may or may not include path). + entity : :obj:`str`, optional + A new BIDS-like key-value pair. + newsuffix : :obj:`str`, optional + Replace the BIDS suffix. + newpath : :obj:`str` or :obj:`os.pathlike`, optional + Path to replace the path of the input orig_path. + newext : :obj:`str`, optional + Replace the extension of the file. + position : :obj:`int`, optional + Position to insert the entity in the filename. + absolute : :obj:`bool`, optional + If True (default), returns the absolute path of the modified filename. + + Returns + ------- + Absolute path of the modified filename + + Examples + -------- + >>> derive_bids_fname( + ... 'sub-001/ses-01/anat/sub-001_ses-01_T1w.nii.gz', + ... entity='desc-preproc', + ... absolute=False, + ... ) + PosixPath('sub-001/ses-01/anat/sub-001_ses-01_desc-preproc_T1w.nii.gz') + + >>> derive_bids_fname( + ... 'sub-001/ses-01/anat/sub-001_ses-01_T1w.nii.gz', + ... entity='desc-brain', + ... newsuffix='mask', + ... newext=".nii", + ... absolute=False, + ... ) # doctest: +ELLIPSIS + PosixPath('sub-001/ses-01/anat/sub-001_ses-01_desc-brain_mask.nii') + + >>> derive_bids_fname( + ... 'sub-001/ses-01/anat/sub-001_ses-01_T1w.nii.gz', + ... entity='desc-brain', + ... newsuffix='mask', + ... newext=".nii", + ... newpath="/output/node", + ... absolute=True, + ... ) # doctest: +ELLIPSIS + PosixPath('/output/node/sub-001_ses-01_desc-brain_mask.nii') + + >>> derive_bids_fname( + ... 'sub-001/ses-01/anat/sub-001_ses-01_T1w.nii.gz', + ... entity='desc-brain', + ... newsuffix='mask', + ... newext=".nii", + ... newpath=".", + ... absolute=False, + ... ) # doctest: +ELLIPSIS + PosixPath('sub-001_ses-01_desc-brain_mask.nii') + + """ + + orig_path = Path(orig_path) + newpath = orig_path.parent if newpath is None else Path(newpath) + + ext = "".join(orig_path.suffixes) + newext = newext if newext is not None else ext + orig_stem = orig_path.name.replace(ext, "") + + suffix = orig_stem.rsplit("_", maxsplit=1)[-1].strip("_") + newsuffix = newsuffix.strip("_") if newsuffix is not None else suffix + + orig_stem = orig_stem.replace(suffix, "").strip("_") + bidts = [bit for bit in orig_stem.split("_") if bit] + if entity: + if position == -1: + bidts.append(entity) + else: + bidts.insert(position, entity.strip("_")) + + retval = (newpath / f"{'_'.join(bidts)}_{newsuffix}.{newext.strip('.')}") + + return retval.absolute() if absolute else retval diff --git a/mriqc/workflows/anatomical/base.py b/mriqc/workflows/anatomical/base.py index 400dee068..0751ee327 100644 --- a/mriqc/workflows/anatomical/base.py +++ b/mriqc/workflows/anatomical/base.py @@ -86,6 +86,7 @@ def anat_qc_workflow(name="anatMRIQC"): wf = anat_qc_workflow() """ + from mriqc.workflows.shared import synthstrip_wf dataset = config.workflow.inputs.get("t1w", []) + config.workflow.inputs.get("t2w", []) @@ -682,101 +683,6 @@ def airmsk_wf(name="AirMaskWorkflow"): return workflow -def synthstrip_wf(name="synthstrip_wf", omp_nthreads=None): - """Create a brain-extraction workflow using SynthStrip.""" - from nipype.interfaces.ants import N4BiasFieldCorrection - from niworkflows.interfaces.nibabel import IntensityClip, ApplyMask - from mriqc.interfaces.synthstrip import SynthStrip - - inputnode = pe.Node(niu.IdentityInterface(fields=["in_files"]), name="inputnode") - outputnode = pe.Node( - niu.IdentityInterface(fields=["out_corrected", "out_brain", "bias_image", "out_mask"]), - name="outputnode", - ) - - # truncate target intensity for N4 correction - pre_clip = pe.Node(IntensityClip(p_min=10, p_max=99.9), name="pre_clip") - - pre_n4 = pe.Node( - N4BiasFieldCorrection( - dimension=3, - num_threads=omp_nthreads, - rescale_intensities=True, - copy_header=True, - ), - name="pre_n4", - ) - - post_n4 = pe.Node( - N4BiasFieldCorrection( - dimension=3, - save_bias=True, - num_threads=omp_nthreads, - n_iterations=[50] * 4, - copy_header=True, - ), - name="post_n4", - ) - - synthstrip = pe.Node( - SynthStrip(num_threads=omp_nthreads), - name="synthstrip", - num_threads=omp_nthreads, - ) - - final_masked = pe.Node(ApplyMask(), name="final_masked") - final_inu = pe.Node(niu.Function(function=_apply_bias_correction), name="final_inu") - - workflow = pe.Workflow(name=name) - # fmt: off - workflow.connect([ - (inputnode, final_inu, [("in_files", "in_file")]), - (inputnode, pre_clip, [("in_files", "in_file")]), - (pre_clip, pre_n4, [("out_file", "input_image")]), - (pre_n4, synthstrip, [("output_image", "in_file")]), - (synthstrip, post_n4, [("out_mask", "weight_image")]), - (synthstrip, final_masked, [("out_mask", "in_mask")]), - (pre_clip, post_n4, [("out_file", "input_image")]), - (post_n4, final_inu, [("bias_image", "bias_image")]), - (post_n4, final_masked, [("output_image", "in_file")]), - (final_masked, outputnode, [("out_file", "out_brain")]), - (post_n4, outputnode, [("bias_image", "bias_image")]), - (synthstrip, outputnode, [("out_mask", "out_mask")]), - (post_n4, outputnode, [("output_image", "out_corrected")]), - ]) - # fmt: on - return workflow - - -def _apply_bias_correction(in_file, bias_image, out_file=None): - import os.path as op - - import numpy as np - import nibabel as nb - - img = nb.load(in_file) - data = np.clip( - img.get_fdata() * nb.load(bias_image).get_fdata(), - a_min=0, - a_max=None, - ) - out_img = img.__class__( - data.astype(img.get_data_dtype()), - img.affine, - img.header, - ) - - if out_file is None: - fname, ext = op.splitext(op.basename(in_file)) - if ext == ".gz": - fname, ext2 = op.splitext(fname) - ext = ext2 + ext - out_file = op.abspath(f"{fname}_inu{ext}") - - out_img.to_filename(out_file) - return out_file - - def _binarize(in_file, threshold=0.5, out_file=None): import os.path as op diff --git a/mriqc/workflows/diffusion/base.py b/mriqc/workflows/diffusion/base.py index 470f02fa7..a81d08520 100644 --- a/mriqc/workflows/diffusion/base.py +++ b/mriqc/workflows/diffusion/base.py @@ -76,6 +76,7 @@ def dmri_qc_workflow(name="dwiMRIQC"): ReadDWIMetadata, WeightedStat, ) + from mriqc.workflows.shared import synthstrip_wf as dmri_bmsk_workflow from mriqc.messages import BUILDING_WORKFLOW workflow = pe.Workflow(name=name) @@ -349,70 +350,6 @@ def compute_iqms(name="ComputeIQMs"): return workflow -def dmri_bmsk_workflow(name="dmri_brainmask", omp_nthreads=None): - """ - Compute a brain mask for the input :abbr:`dMRI (diffusion MRI)` dataset. - - .. workflow:: - - from mriqc.workflows.diffusion.base import dmri_bmsk_workflow - from mriqc.testing import mock_config - with mock_config(): - wf = dmri_bmsk_workflow() - - - """ - - from nipype.interfaces.ants import N4BiasFieldCorrection - from niworkflows.interfaces.nibabel import ApplyMask - from mriqc.interfaces.synthstrip import SynthStrip - from mriqc.workflows.anatomical.base import _apply_bias_correction - - inputnode = pe.Node(niu.IdentityInterface(fields=["in_file"]), name="inputnode") - outputnode = pe.Node( - niu.IdentityInterface(fields=["out_corrected", "out_brain", "bias_image", "out_mask"]), - name="outputnode", - ) - - post_n4 = pe.Node( - N4BiasFieldCorrection( - dimension=3, - save_bias=True, - num_threads=omp_nthreads, - n_iterations=[50] * 4, - copy_header=True, - ), - name="post_n4", - ) - - synthstrip = pe.Node( - SynthStrip(num_threads=omp_nthreads), - name="synthstrip", - num_threads=omp_nthreads, - ) - - final_masked = pe.Node(ApplyMask(), name="final_masked") - final_inu = pe.Node(niu.Function(function=_apply_bias_correction), name="final_inu") - - workflow = pe.Workflow(name=name) - # fmt: off - workflow.connect([ - (inputnode, final_inu, [("in_file", "in_file")]), - (inputnode, synthstrip, [("in_file", "in_file")]), - (inputnode, post_n4, [("in_file", "input_image")]), - (synthstrip, post_n4, [("out_mask", "weight_image")]), - (synthstrip, final_masked, [("out_mask", "in_mask")]), - (post_n4, final_inu, [("bias_image", "bias_image")]), - (post_n4, final_masked, [("output_image", "in_file")]), - (final_masked, outputnode, [("out_file", "out_brain")]), - (post_n4, outputnode, [("bias_image", "bias_image")]), - (synthstrip, outputnode, [("out_mask", "out_mask")]), - (post_n4, outputnode, [("output_image", "out_corrected")]), - ]) - # fmt: on - return workflow - - def init_dmriref_wf( in_file=None, name="init_dmriref_wf", diff --git a/mriqc/workflows/functional/base.py b/mriqc/workflows/functional/base.py index 1ede04f61..3db1cb224 100644 --- a/mriqc/workflows/functional/base.py +++ b/mriqc/workflows/functional/base.py @@ -45,6 +45,8 @@ from mriqc import config from nipype.interfaces import utility as niu from nipype.pipeline import engine as pe +from niworkflows.utils.connections import pop_file as _pop + from mriqc.interfaces.datalad import DataladIdentityInterface from mriqc.workflows.functional.output import init_func_report_wf @@ -64,8 +66,10 @@ def fmri_qc_workflow(name="funcMRIQC"): """ from nipype.algorithms.confounds import TSNR, NonSteadyStateDetector from nipype.interfaces.afni import TStat + from niworkflows.interfaces.bids import ReadSidecarJSON from niworkflows.interfaces.header import SanitizeImage from mriqc.messages import BUILDING_WORKFLOW + from mriqc.interfaces.functional import SelectEcho workflow = pe.Workflow(name=name) @@ -76,7 +80,7 @@ def fmri_qc_workflow(name="funcMRIQC"): message = BUILDING_WORKFLOW.format( modality="functional", detail=( - f"for {len(dataset)} NIfTI files." + f"for {len(dataset)} BOLD runs." if len(dataset) > 2 else f"({' and '.join(('<%s>' % v for v in dataset))})." ), @@ -88,9 +92,10 @@ def fmri_qc_workflow(name="funcMRIQC"): inputnode = pe.Node(niu.IdentityInterface(fields=["in_file"]), name="inputnode") inputnode.iterables = [("in_file", dataset)] - datalad_get = pe.Node( + datalad_get = pe.MapNode( DataladIdentityInterface(fields=["in_file"], dataset_path=config.execution.bids_dir), name="datalad_get", + iterfield=["in_file"], ) outputnode = pe.Node( @@ -98,32 +103,49 @@ def fmri_qc_workflow(name="funcMRIQC"): name="outputnode", ) + # Get metadata + meta = pe.MapNode(ReadSidecarJSON( + index_db=config.execution.bids_database_dir + ), name="metadata", iterfield=["in_file"]) + + pick_echo = pe.Node(SelectEcho(), name="pick_echo") + non_steady_state_detector = pe.Node(NonSteadyStateDetector(), name="non_steady_state_detector") - sanitize = pe.Node(SanitizeImage(), name="sanitize", mem_gb=mem_gb * 4.0) - sanitize.inputs.max_32bit = config.execution.float32 + sanitize = pe.MapNode( + SanitizeImage(max_32bit=config.execution.float32), + name="sanitize", + mem_gb=mem_gb * 4.0, + iterfield=["in_file"], + ) # Workflow -------------------------------------------------------- # 1. HMC: head motion correct - hmcwf = hmc() + hmcwf = hmc(omp_nthreads=config.nipype.omp_nthreads) # Set HMC settings hmcwf.inputs.inputnode.fd_radius = config.workflow.fd_radius # 2. Compute mean fmri - mean = pe.Node( + mean = pe.MapNode( TStat(options="-mean", outputtype="NIFTI_GZ"), name="mean", mem_gb=mem_gb * 1.5, + iterfield=["in_file"], + ) + + # Compute TSNR using nipype implementation + tsnr = pe.MapNode( + TSNR(), + name="compute_tsnr", + mem_gb=mem_gb * 2.5, + iterfield=["in_file"], ) # EPI to MNI registration ema = epi_mni_align() - # Compute TSNR using nipype implementation - tsnr = pe.Node(TSNR(), name="compute_tsnr", mem_gb=mem_gb * 2.5) - # 7. Compute IQMs iqmswf = compute_iqms() # Reports @@ -133,22 +155,35 @@ def fmri_qc_workflow(name="funcMRIQC"): workflow.connect([ (inputnode, datalad_get, [("in_file", "in_file")]), - (inputnode, func_report_wf, [ - ("in_file", "inputnode.name_source"), - ]), - (datalad_get, iqmswf, [("in_file", "inputnode.in_file")]), + (datalad_get, meta, [("in_file", "in_file")]), + (datalad_get, pick_echo, [("in_file", "in_files")]), (datalad_get, sanitize, [("in_file", "in_file")]), - (datalad_get, non_steady_state_detector, [("in_file", "in_file")]), + (meta, pick_echo, [("out_dict", "metadata")]), + (pick_echo, non_steady_state_detector, [("out_file", "in_file")]), (non_steady_state_detector, sanitize, [("n_volumes_to_discard", "n_volumes_to_discard")]), (sanitize, hmcwf, [("out_file", "inputnode.in_file")]), (hmcwf, mean, [("outputnode.out_file", "in_file")]), (hmcwf, tsnr, [("outputnode.out_file", "in_file")]), - (mean, ema, [("out_file", "inputnode.epi_mean")]), + (mean, ema, [(("out_file", _pop), "inputnode.epi_mean")]), + # Feed IQMs computation + (meta, iqmswf, [("out_dict", "inputnode.metadata"), + ("subject", "inputnode.subject"), + ("session", "inputnode.session"), + ("task", "inputnode.task"), + ("acquisition", "inputnode.acquisition"), + ("reconstruction", "inputnode.reconstruction"), + ("run", "inputnode.run")]), + (datalad_get, iqmswf, [("in_file", "inputnode.in_file")]), (sanitize, iqmswf, [("out_file", "inputnode.in_ras")]), (mean, iqmswf, [("out_file", "inputnode.epi_mean")]), (hmcwf, iqmswf, [("outputnode.out_file", "inputnode.hmc_epi"), ("outputnode.out_fd", "inputnode.hmc_fd")]), (tsnr, iqmswf, [("tsnr_file", "inputnode.in_tsnr")]), + (non_steady_state_detector, iqmswf, [("n_volumes_to_discard", "inputnode.exclude_index")]), + # Feed reportlet generation + (inputnode, func_report_wf, [ + ("in_file", "inputnode.name_source"), + ]), (sanitize, func_report_wf, [("out_file", "inputnode.in_ras")]), (mean, func_report_wf, [("out_file", "inputnode.epi_mean")]), (tsnr, func_report_wf, [("stddev_file", "inputnode.in_stddev")]), @@ -160,13 +195,12 @@ def fmri_qc_workflow(name="funcMRIQC"): ("outputnode.epi_parc", "inputnode.epi_parc"), ("outputnode.report", "inputnode.mni_report"), ]), - (non_steady_state_detector, iqmswf, [("n_volumes_to_discard", "inputnode.exclude_index")]), (iqmswf, func_report_wf, [ ("outputnode.out_file", "inputnode.in_iqms"), ("outputnode.out_dvars", "inputnode.in_dvars"), ("outputnode.outliers", "inputnode.outliers"), - ("outputnode.meta_sidecar", "inputnode.meta_sidecar"), ]), + (meta, func_report_wf, [("out_dict", "inputnode.meta_sidecar")]), (hmcwf, outputnode, [("outputnode.out_fd", "out_fd")]), ]) # fmt: on @@ -183,13 +217,15 @@ def fmri_qc_workflow(name="funcMRIQC"): # population specific changes to brain masking if config.workflow.species == "human": - skullstrip_epi = fmri_bmsk_workflow() + from mriqc.workflows.shared import synthstrip_wf as fmri_bmsk_workflow + + skullstrip_epi = fmri_bmsk_workflow(omp_nthreads=config.nipype.omp_nthreads) # fmt: off workflow.connect([ - (mean, skullstrip_epi, [("out_file", "inputnode.in_file")]), - (skullstrip_epi, ema, [("outputnode.out_file", "inputnode.epi_mask")]), - (skullstrip_epi, iqmswf, [("outputnode.out_file", "inputnode.brainmask")]), - (skullstrip_epi, func_report_wf, [("outputnode.out_file", "inputnode.brainmask")]), + (mean, skullstrip_epi, [(("out_file", _pop), "inputnode.in_files")]), + (skullstrip_epi, ema, [("outputnode.out_mask", "inputnode.epi_mask")]), + (skullstrip_epi, iqmswf, [("outputnode.out_mask", "inputnode.brainmask")]), + (skullstrip_epi, func_report_wf, [("outputnode.out_mask", "inputnode.brainmask")]), ]) # fmt: on else: @@ -216,11 +252,11 @@ def fmri_qc_workflow(name="funcMRIQC"): if not config.execution.no_sub: from mriqc.interfaces.webapi import UploadIQMs - upldwf = pe.Node(UploadIQMs( + upldwf = pe.MapNode(UploadIQMs( endpoint=config.execution.webapi_url, auth_token=config.execution.webapi_token, strict=config.execution.upload_strict, - ), name="UploadMetrics") + ), name="UploadMetrics", iterfield=["in_iqms"]) # fmt: off workflow.connect([ @@ -245,7 +281,6 @@ def compute_iqms(name="ComputeIQMs"): """ from nipype.algorithms.confounds import ComputeDVARS from nipype.interfaces.afni import OutlierCount, QualityIndex - from niworkflows.interfaces.bids import ReadSidecarJSON from mriqc.interfaces import FunctionalQC, IQMFileSink from mriqc.interfaces.reports import AddProvenance @@ -268,6 +303,12 @@ def compute_iqms(name="ComputeIQMs"): "in_tsnr", "metadata", "exclude_index", + "subject", + "session", + "task", + "acquisition", + "reconstruction", + "run", ] ), name="inputnode", @@ -280,7 +321,6 @@ def compute_iqms(name="ComputeIQMs"): "outliers", "out_spikes", "out_fft", - "meta_sidecar", ] ), name="outputnode", @@ -290,32 +330,40 @@ def compute_iqms(name="ComputeIQMs"): inputnode.inputs.fd_thres = config.workflow.fd_thres # Compute DVARS - dvnode = pe.Node( + dvnode = pe.MapNode( ComputeDVARS(save_plot=False, save_all=True), name="ComputeDVARS", mem_gb=mem_gb * 3, + iterfield=["in_file"], ) # AFNI quality measures - fwhm_interface = get_fwhmx() - fwhm = pe.Node(fwhm_interface, name="smoothness") - # fwhm.inputs.acf = True # add when AFNI >= 16 - outliers = pe.Node( + fwhm = pe.MapNode(get_fwhmx(), name="smoothness", iterfield=["in_file"]) + fwhm.inputs.acf = True # Only AFNI >= 16 + + outliers = pe.MapNode( OutlierCount(fraction=True, out_file="outliers.out"), name="outliers", mem_gb=mem_gb * 2.5, + iterfield=["in_file"], ) - quality = pe.Node( + quality = pe.MapNode( QualityIndex(automask=True), out_file="quality.out", name="quality", mem_gb=mem_gb * 3, + iterfield=["in_file"], ) - gcor = pe.Node(GCOR(), name="gcor", mem_gb=mem_gb * 2) + gcor = pe.MapNode(GCOR(), name="gcor", mem_gb=mem_gb * 2, iterfield=["in_file"]) - measures = pe.Node(FunctionalQC(), name="measures", mem_gb=mem_gb * 3) + measures = pe.MapNode( + FunctionalQC(), + name="measures", + mem_gb=mem_gb * 3, + iterfield=["in_epi", "in_hmc", "in_tsnr", "in_dvars", "in_fwhm"], + ) # fmt: off workflow.connect([ @@ -341,19 +389,15 @@ def compute_iqms(name="ComputeIQMs"): ]) # fmt: on - # Add metadata - meta = pe.Node(ReadSidecarJSON( - index_db=config.execution.bids_database_dir - ), name="metadata") - - addprov = pe.Node( + addprov = pe.MapNode( AddProvenance(modality="bold"), name="provenance", run_without_submitting=True, + iterfield=["in_file"], ) # Save to JSON file - datasink = pe.Node( + datasink = pe.MapNode( IQMFileSink( modality="bold", out_dir=str(config.execution.output_dir), @@ -361,28 +405,27 @@ def compute_iqms(name="ComputeIQMs"): ), name="datasink", run_without_submitting=True, + iterfield=["in_file", "root", "metadata", "provenance"], ) # fmt: off workflow.connect([ - (inputnode, datasink, [("in_file", "in_file"), - ("exclude_index", "dummy_trs")]), - (inputnode, meta, [("in_file", "in_file")]), (inputnode, addprov, [("in_file", "in_file")]), - (meta, datasink, [("subject", "subject_id"), - ("session", "session_id"), - ("task", "task_id"), - ("acquisition", "acq_id"), - ("reconstruction", "rec_id"), - ("run", "run_id"), - ("out_dict", "metadata")]), + (inputnode, datasink, [("in_file", "in_file"), + ("exclude_index", "dummy_trs"), + (("subject", _pop), "subject_id"), + (("session", _pop), "session_id"), + (("task", _pop), "task_id"), + (("acquisition", _pop), "acq_id"), + (("reconstruction", _pop), "rec_id"), + (("run", _pop), "run_id"), + ("metadata", "metadata")]), (addprov, datasink, [("out_prov", "provenance")]), (outliers, datasink, [(("out_file", _parse_tout), "aor")]), (gcor, datasink, [(("out", _tofloat), "gcor")]), (quality, datasink, [(("out_file", _parse_tqual), "aqi")]), (measures, datasink, [("out_qc", "root")]), (datasink, outputnode, [("out_file", "out_file")]), - (meta, outputnode, [("out_dict", "meta_sidecar")]), ]) # fmt: on @@ -390,13 +433,14 @@ def compute_iqms(name="ComputeIQMs"): if config.workflow.fft_spikes_detector: from mriqc.workflows.utils import slice_wise_fft - spikes_fft = pe.Node( + spikes_fft = pe.MapNode( niu.Function( input_names=["in_file"], output_names=["n_spikes", "out_spikes", "out_fft"], function=slice_wise_fft, ), name="SpikesFinderFFT", + iterfield=["in_file"], ) # fmt: off @@ -441,7 +485,7 @@ def fmri_bmsk_workflow(name="fMRIBrainMask"): return workflow -def hmc(name="fMRI_HMC"): +def hmc(name="fMRI_HMC", omp_nthreads=None): """ Create a :abbr:`HMC (head motion correction)` workflow for fMRI. @@ -468,9 +512,9 @@ def hmc(name="fMRI_HMC"): outputnode = pe.Node(niu.IdentityInterface(fields=["out_file", "out_fd"]), name="outputnode") # calculate hmc parameters - hmc = pe.Node( + estimate_hm = pe.Node( Volreg(args="-Fourier -twopass", zpad=4, outputtype="NIFTI_GZ"), - name="motion_correct", + name="estimate_hm", mem_gb=mem_gb * 2.5, ) @@ -480,49 +524,73 @@ def hmc(name="fMRI_HMC"): name="ComputeFD", ) + # Apply transforms to other echos + apply_hmc = pe.MapNode( + niu.Function(function=_apply_transforms, input_names=["in_file", "in_xfm"]), + name="apply_hmc", + iterfield=["in_file"], + # NiTransforms is a memory hog, so ensure only one process is running at a time + num_threads=config.environment.cpu_count, + ) + # fmt: off workflow.connect([ (inputnode, fdnode, [("fd_radius", "radius")]), - (hmc, outputnode, [("out_file", "out_file")]), - (hmc, fdnode, [("oned_file", "in_file")]), + (estimate_hm, apply_hmc, [("oned_matrix_save", "in_xfm")]), + (apply_hmc, outputnode, [("out", "out_file")]), + (estimate_hm, fdnode, [("oned_file", "in_file")]), (fdnode, outputnode, [("out_file", "out_fd")]), ]) # fmt: on - # despiking, and deoblique - - deoblique_node = pe.Node(Refit(deoblique=True), name="deoblique") - - despike_node = pe.Node(Despike(outputtype="NIFTI_GZ"), name="despike") + if not (config.workflow.despike or config.workflow.deoblique): + # fmt: off + workflow.connect([ + (inputnode, estimate_hm, [(("in_file", _pop), "in_file")]), + (inputnode, apply_hmc, [("in_file", "in_file")]), + ]) + # fmt: on + return workflow + # despiking, and deoblique + deoblique_node = pe.MapNode( + Refit(deoblique=True), + name="deoblique", + iterfield=["in_file"], + ) + despike_node = pe.MapNode( + Despike(outputtype="NIFTI_GZ"), + name="despike", + iterfield=["in_file"], + ) if config.workflow.despike and config.workflow.deoblique: # fmt: off workflow.connect([ (inputnode, despike_node, [("in_file", "in_file")]), (despike_node, deoblique_node, [("out_file", "in_file")]), - (deoblique_node, hmc, [("out_file", "in_file")]), + (deoblique_node, estimate_hm, [(("out_file", _pop), "in_file")]), + (deoblique_node, apply_hmc, [("out_file", "in_file")]), ]) # fmt: on elif config.workflow.despike: # fmt: off workflow.connect([ (inputnode, despike_node, [("in_file", "in_file")]), - (despike_node, hmc, [("out_file", "in_file")]), + (despike_node, estimate_hm, [(("out_file", _pop), "in_file")]), + (despike_node, apply_hmc, [("out_file", "in_file")]), ]) # fmt: on elif config.workflow.deoblique: # fmt: off workflow.connect([ (inputnode, deoblique_node, [("in_file", "in_file")]), - (deoblique_node, hmc, [("out_file", "in_file")]), + (deoblique_node, estimate_hm, [(("out_file", _pop), "in_file")]), + (deoblique_node, apply_hmc, [("out_file", "in_file")]), ]) # fmt: on else: - # fmt: off - workflow.connect([ - (inputnode, hmc, [("in_file", "in_file")]), - ]) - # fmt: on + raise NotImplementedError + return workflow @@ -665,13 +733,13 @@ def epi_mni_align(name="SpatialNormalization"): return workflow -def _mean(inlist): - import numpy as np - - return np.mean(inlist) - - def _parse_tqual(in_file): + if isinstance(in_file, (list, tuple)): + return ( + [_parse_tqual(f) for f in in_file] if len(in_file) > 1 + else _parse_tqual(in_file[0]) + ) + import numpy as np with open(in_file, "r") as fin: @@ -680,7 +748,30 @@ def _parse_tqual(in_file): def _parse_tout(in_file): + if isinstance(in_file, (list, tuple)): + return ( + [_parse_tout(f) for f in in_file] if len(in_file) > 1 + else _parse_tout(in_file[0]) + ) + import numpy as np data = np.loadtxt(in_file) # pylint: disable=no-member return data.mean() + + +def _apply_transforms(in_file, in_xfm): + from pathlib import Path + from nitransforms.linear import load + from mriqc.utils.bids import derive_bids_fname + + realigned = load(in_xfm, fmt="afni", reference=in_file, moving=in_file).apply(in_file) + out_file = derive_bids_fname( + in_file, + entity="desc-realigned", + newpath=Path.cwd(), + absolute=True, + ) + + realigned.to_filename(out_file) + return str(out_file) diff --git a/mriqc/workflows/functional/output.py b/mriqc/workflows/functional/output.py index fc6a421ce..f620cb2a1 100644 --- a/mriqc/workflows/functional/output.py +++ b/mriqc/workflows/functional/output.py @@ -81,7 +81,7 @@ def init_func_report_wf(name="func_report_wf"): # Set FD threshold inputnode.inputs.fd_thres = config.workflow.fd_thres - spmask = pe.Node( + spmask = pe.MapNode( niu.Function( input_names=["in_file", "in_mask"], output_names=["out_file", "out_plot"], @@ -89,12 +89,14 @@ def init_func_report_wf(name="func_report_wf"): ), name="SpikesMask", mem_gb=mem_gb * 3.5, + iterfield=["in_file"], ) - spikes_bg = pe.Node( + spikes_bg = pe.MapNode( Spikes(no_zscore=True, detrend=False), name="SpikesFinderBgMask", mem_gb=mem_gb * 2.5, + iterfield=["in_file", "in_mask"], ) # Generate crown mask @@ -103,7 +105,12 @@ def init_func_report_wf(name="func_report_wf"): subtract_mask = pe.Node(BinarySubtraction(), name="subtract_mask") parcels = pe.Node(niu.Function(function=_carpet_parcellation), name="parcels") - bigplot = pe.Node(FMRISummary(), name="BigPlot", mem_gb=mem_gb * 3.5) + bigplot = pe.MapNode( + FMRISummary(), + name="BigPlot", + mem_gb=mem_gb * 3.5, + iterfield=["in_func", "dvars", "outliers", "in_spikes_bg"], + ) # fmt: off workflow.connect([ @@ -118,58 +125,66 @@ def init_func_report_wf(name="func_report_wf"): (inputnode, parcels, [("epi_parc", "segmentation")]), (inputnode, dilated_mask, [("brainmask", "in_mask")]), (inputnode, subtract_mask, [("brainmask", "in_subtract")]), + (spmask, spikes_bg, [("out_file", "in_mask")]), (dilated_mask, subtract_mask, [("out_mask", "in_base")]), (subtract_mask, parcels, [("out_mask", "crown_mask")]), (parcels, bigplot, [("out", "in_segm")]), (spikes_bg, bigplot, [("out_tsz", "in_spikes_bg")]), - (spmask, spikes_bg, [("out_file", "in_mask")]), ]) # fmt: on - mosaic_mean = pe.Node( + mosaic_mean = pe.MapNode( PlotMosaic( out_file="plot_func_mean_mosaic1.svg", cmap="Greys_r", ), name="PlotMosaicMean", + iterfield=["in_file"], ) - mosaic_stddev = pe.Node( + mosaic_stddev = pe.MapNode( PlotMosaic( out_file="plot_func_stddev_mosaic2_stddev.svg", cmap="viridis", ), name="PlotMosaicSD", + iterfield=["in_file"], ) - ds_report_mean = pe.Node( + ds_report_mean = pe.MapNode( DerivativesDataSink( base_directory=reportlets_dir, desc="mean", datatype="figures", + dismiss_entities=("part",) ), name="ds_report_mean", run_without_submitting=True, + iterfield=["in_file", "source_file"], ) - ds_report_stdev = pe.Node( + ds_report_stdev = pe.MapNode( DerivativesDataSink( base_directory=reportlets_dir, desc="stdev", datatype="figures", + dismiss_entities=("part",) ), name="ds_report_stdev", run_without_submitting=True, + iterfield=["in_file", "source_file"], ) - ds_report_carpet = pe.Node( + ds_report_carpet = pe.MapNode( DerivativesDataSink( base_directory=reportlets_dir, desc="carpet", datatype="figures", + dismiss_entities=("part",) ), name="ds_report_carpet", run_without_submitting=True, + iterfield=["in_file", "source_file"], ) # fmt: off @@ -201,6 +216,7 @@ def init_func_report_wf(name="func_report_wf"): base_directory=reportlets_dir, desc="spikes", datatype="figures", + dismiss_entities=("part",) ), name="ds_report_spikes", run_without_submitting=True, @@ -220,21 +236,24 @@ def init_func_report_wf(name="func_report_wf"): return workflow # Verbose-reporting goes here + from niworkflows.utils.connections import pop_file as _pop from nireports.interfaces import PlotContours - mosaic_zoom = pe.Node( + mosaic_zoom = pe.MapNode( PlotMosaic( cmap="Greys_r", ), name="PlotMosaicZoomed", + iterfield=["in_file"], ) - mosaic_noise = pe.Node( + mosaic_noise = pe.MapNode( PlotMosaic( only_noise=True, cmap="viridis_r", ), name="PlotMosaicNoise", + iterfield=["in_file"] ) if config.workflow.species.lower() in ("rat", "mouse"): @@ -254,24 +273,28 @@ def init_func_report_wf(name="func_report_wf"): name="PlotBrainmask", ) - ds_report_zoomed = pe.Node( + ds_report_zoomed = pe.MapNode( DerivativesDataSink( base_directory=reportlets_dir, desc="zoomed", datatype="figures", + dismiss_entities=("part",) ), name="ds_report_zoomed", run_without_submitting=True, + iterfield=["in_file", "source_file"], ) - ds_report_background = pe.Node( + ds_report_background = pe.MapNode( DerivativesDataSink( base_directory=reportlets_dir, desc="background", datatype="figures", + dismiss_entities=("part",) ), name="ds_report_background", run_without_submitting=True, + iterfield=["in_file", "source_file"], ) ds_report_bmask = pe.Node( @@ -279,6 +302,7 @@ def init_func_report_wf(name="func_report_wf"): base_directory=reportlets_dir, desc="brainmask", datatype="figures", + dismiss_entities=("part", "echo"), ), name="ds_report_bmask", run_without_submitting=True, @@ -289,6 +313,7 @@ def init_func_report_wf(name="func_report_wf"): base_directory=reportlets_dir, desc="norm", datatype="figures", + dismiss_entities=("part", "echo"), ), name="ds_report_norm", run_without_submitting=True, @@ -298,7 +323,7 @@ def init_func_report_wf(name="func_report_wf"): workflow.connect([ (inputnode, ds_report_norm, [("mni_report", "in_file"), ("name_source", "source_file")]), - (inputnode, plot_bmask, [("epi_mean", "in_file"), + (inputnode, plot_bmask, [(("epi_mean", _pop), "in_file"), ("brainmask", "in_contours")]), (inputnode, mosaic_zoom, [("epi_mean", "in_file"), ("brainmask", "bbox_mask_file")]), @@ -308,7 +333,7 @@ def init_func_report_wf(name="func_report_wf"): (inputnode, ds_report_bmask, [("name_source", "source_file")]), (mosaic_zoom, ds_report_zoomed, [("out_file", "in_file")]), (mosaic_noise, ds_report_background, [("out_file", "in_file")]), - (plot_bmask, ds_report_bmask, [("out_file", "in_file")]), + (plot_bmask, ds_report_bmask, [(("out_file", _pop), "in_file")]), ]) # fmt: on @@ -398,4 +423,7 @@ def _carpet_parcellation(segmentation, crown_mask): def _get_tr(meta_dict): + if isinstance(meta_dict, (list, tuple)): + meta_dict = meta_dict[0] + return meta_dict.get("RepetitionTime", None) diff --git a/mriqc/workflows/shared.py b/mriqc/workflows/shared.py new file mode 100644 index 000000000..c2c778a07 --- /dev/null +++ b/mriqc/workflows/shared.py @@ -0,0 +1,88 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +# +# Copyright 2023 The NiPreps Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# We support and encourage derived works from this project, please read +# about our expectations at +# +# https://www.nipreps.org/community/licensing/ +# +"""Shared workflows.""" +from nipype.interfaces import utility as niu +from nipype.pipeline import engine as pe + + +def synthstrip_wf(name="synthstrip_wf", omp_nthreads=None): + """Create a brain-extraction workflow using SynthStrip.""" + from nipype.interfaces.ants import N4BiasFieldCorrection + from niworkflows.interfaces.nibabel import IntensityClip, ApplyMask + from mriqc.interfaces.synthstrip import SynthStrip + + inputnode = pe.Node(niu.IdentityInterface(fields=["in_files"]), name="inputnode") + outputnode = pe.Node( + niu.IdentityInterface(fields=["out_corrected", "out_brain", "bias_image", "out_mask"]), + name="outputnode", + ) + + # truncate target intensity for N4 correction + pre_clip = pe.Node(IntensityClip(p_min=10, p_max=99.9), name="pre_clip") + + pre_n4 = pe.Node( + N4BiasFieldCorrection( + dimension=3, + num_threads=omp_nthreads, + rescale_intensities=True, + copy_header=True, + ), + name="pre_n4", + ) + + post_n4 = pe.Node( + N4BiasFieldCorrection( + dimension=3, + save_bias=True, + num_threads=omp_nthreads, + n_iterations=[50] * 4, + copy_header=True, + ), + name="post_n4", + ) + + synthstrip = pe.Node( + SynthStrip(num_threads=omp_nthreads), + name="synthstrip", + num_threads=omp_nthreads, + ) + + final_masked = pe.Node(ApplyMask(), name="final_masked") + + workflow = pe.Workflow(name=name) + # fmt: off + workflow.connect([ + (inputnode, pre_clip, [("in_files", "in_file")]), + (pre_clip, pre_n4, [("out_file", "input_image")]), + (pre_n4, synthstrip, [("output_image", "in_file")]), + (synthstrip, post_n4, [("out_mask", "weight_image")]), + (synthstrip, final_masked, [("out_mask", "in_mask")]), + (pre_clip, post_n4, [("out_file", "input_image")]), + (post_n4, final_masked, [("output_image", "in_file")]), + (final_masked, outputnode, [("out_file", "out_brain")]), + (post_n4, outputnode, [("bias_image", "bias_image")]), + (synthstrip, outputnode, [("out_mask", "out_mask")]), + (post_n4, outputnode, [("output_image", "out_corrected")]), + ]) + # fmt: on + return workflow diff --git a/mriqc/workflows/utils.py b/mriqc/workflows/utils.py index 690c88a67..10c30f100 100644 --- a/mriqc/workflows/utils.py +++ b/mriqc/workflows/utils.py @@ -25,9 +25,11 @@ def _tofloat(inlist): if isinstance(inlist, (list, tuple)): - return [float(el) for el in inlist] - else: - return float(inlist) + return ( + [_tofloat(el) for el in inlist] if len(inlist) > 1 + else _tofloat(inlist[0]) + ) + return float(inlist) def fwhm_dict(fwhm):