Merge pull request #1099 from nipreps/enh/drop-fsl-fast

ENH: Replace FSL FAST with ANTs Atropos for brain tissue segmentation
nipreps · Apr 4, 2023 · b3c4813 · b3c4813
2 parents ad44d58 + eef5ed6
commit b3c4813
Show file tree

Hide file tree

Showing 5 changed files with 171 additions and 197 deletions.
diff --git a/mriqc/data/testdata/group_T1w.tsv b/mriqc/data/testdata/group_T1w.tsv
@@ -1,5 +1,5 @@
 bids_name	cjv	cnr	efc	fber	fwhm_avg	fwhm_x	fwhm_y	fwhm_z	icvs_csf	icvs_gm	icvs_wm	inu_med	inu_range	qi_1	qi_2	rpve_csf	rpve_gm	rpve_wm	size_x	size_y	size_z	snr_csf	snr_gm	snr_total	snr_wm	snrd_csf	snrd_gm	snrd_total	snrd_wm	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_csf_k	summary_csf_mad	summary_csf_mean	summary_csf_median	summary_csf_n	summary_csf_p05	summary_csf_p95	summary_csf_stdv	summary_gm_k	summary_gm_mad	summary_gm_mean	summary_gm_median	summary_gm_n	summary_gm_p05	summary_gm_p95	summary_gm_stdv	summary_wm_k	summary_wm_mad	summary_wm_mean	summary_wm_median	summary_wm_n	summary_wm_p05	summary_wm_p95	summary_wm_stdv	tpm_overlap_csf	tpm_overlap_gm	tpm_overlap_wm	wm2max
-sub-50137_T1w	0.3733171053791084	3.5963188952818133	0.6960134252231683	7632.404958677686	4.202703333333333	4.04437	4.57298	3.99076	0.1806197891725539	0.4512949612224343	0.3680852496050118	0.6880311369895935	0.3810786575078964	0.0	0.0016184252287961387	18.825199358805268	8.39010791744813	11.649075350809115	208	256	176	1.699438922766317	9.605216897014863	9.58064136132585	17.43726826419637	6.293807552347417	19.123357320614378	18.04185872789903	28.708411310735297	1.0	1.0	1.0	186.37771714301948	0.0	6.5942646190456005	0.0	2639142.0	0.0	42.62682022154331	22.821292323214568	19.010185895313597	99.59667629101651	240.167928369194	219.24109540879726	23305.0	76.80156033039093	460.7849340051412	129.0051608717524	0.026768181310648842	69.05369556177145	666.3318587159577	666.1509383618832	48307.0	551.0137391388416	780.3517317920923	69.35231661036728	0.3715126186254438	55.593053856985584	1001.9872130846554	1000.0406734496355	249477.0	911.7335358560085	1099.2182608991861	57.35064999984544	0.17070114849123164	0.46718816687526277	0.5031832025227133	0.41705010441957996
-sub-50152_T1w	0.31520442759890505	4.226655470805463	0.6407575000929343	6400.0	3.538223004967018	3.5068090149010525	3.72369	3.38417	0.14843869455707587	0.5250711585882605	0.3264901468546637	0.6055912971496582	0.2994886189699173	0.0	0.006535249887689515	24.026460235258366	7.935436527754754	13.147359602746743	160	239	200	2.5109487264748225	11.946763921301619	11.61657785576157	20.39202091950827	15.362286564812278	31.939089811663475	31.477273591826457	47.130444399003615	1.100000023841858	1.0	1.0	6.5072174689779505	0.0	7.6793006426825645	0.0	1890843.0	0.0	37.758178912103176	13.901023153001923	0.796359241530658	114.96067825895601	334.59021023421326	325.96494595706463	13920.0	121.97724339254201	557.9410418082027	129.8127810412912	0.08289860847729758	55.301123908911876	678.5826404873823	677.7001353576779	56601.0	585.8817124217749	774.0235786288977	56.726168956794645	0.5065109326668424	46.301248062937475	1000.1088627206178	1000.0381581634283	147620.0	920.0549581423402	1080.6341173063959	49.04049357838929	0.1467541776658981	0.48341101783423046	0.4780148060252424	0.49068974560713347
-sub-50785_T1w	0.37382553166391375	3.3598989753118196	0.7530480076783987	-1.0	3.499472057939973	3.46213	3.7065461738199192	3.32974	0.1898804810921801	0.4421335942413909	0.36798592466642893	0.9646123945713043	0.2733038604259488	0.0	0.0010024604792166635	26.059428378074543	11.501432093376257	15.490493983987355	256	180	256	1.6733726682386914	7.93970382721396	6.544238472612016	10.019638922383399	3.7285399580690743	11.069259960832369	11.391272766670696	19.376018381110644	1.0	0.9999875426292419	1.0	953.6855162280638	0.0	1.7893970473216114	0.0	4819101.0	0.0	0.0	33.8117145720979	0.9429151469101642	118.30840362528973	205.39651663696313	192.4306584596634	20111.0	39.272444665431976	405.72885352373123	114.99283921533977	0.1813288386977976	69.94736153698776	571.6101821171164	571.2866180539131	20999.0	452.994039952755	692.0850287079811	71.95142635217684	1.2097531340465286	90.31664813841826	1009.8780921704935	1000.0000046491623	178124.0	863.9937826395035	1189.4063824415207	99.80371601794556	0.16805108018430504	0.46340815878321173	0.48584510178683277	0.38380923954498725
-sub-51187_T1w	0.33878713442707	3.868874272644968	0.6068733970705135	4668.127902703063	3.5394978356516567	4.25443345257453	2.965211919410108	3.398848134970332	0.16404346098394018	0.4935346261940047	0.34242191282205514	0.8655245900154114	0.6595467209815977	0.0	0.002621822353053779	26.046542016972328	9.521649650671295	14.57417618406793	256	132	256	0.5218688172129379	7.961581995404261	7.531157292291048	14.110021064255942	5.018050930340288	11.291229852084964	11.87974542138791	19.329955481738477	0.8593999743461609	1.5000007152557373	0.8593999743461609	341.30063924084027	0.0	11.137474288870767	0.0	2250419.0	0.0	60.959460843354464	33.89228629663162	3.8791461166165755	170.97945856656423	434.6798289708983	259.59970897994936	5336.0	45.10847321804613	1694.362357551232	497.3958473249557	0.00554510270675701	73.69218439433283	582.6331838506209	584.1311744973063	47879.0	459.5589629944414	698.9254273686557	73.36796565376363	0.6089280266844308	67.19881568745514	1002.0900624234341	999.9999775439501	138605.0	889.6923817321658	1121.053142476827	70.87135912910031	0.13692845333498696	0.47923757751556184	0.4788353229114439	0.340853963466229
+sub-50137_T1w	0.5946364030347385	1.6705988220313406	0.6960134252231684	7632.404958677686	4.202703333333333	4.04437	4.57298	3.99076	0.2607070742917811	0.38396092314479024	0.3553320025634285	0.6877774000167847	0.38117922544479355	0.0	0.002295668304182044	5.081432626162044	4.991782430776381	5.114091002400782	208	256	176	2.2951325004992604	4.936353477427729	5.434675057986866	9.07253919603361	10.199045584302324	14.572319191529168	15.232386344514095	20.925794257710795	1.0	1.0	1.0	338.4370104980457	0.0	7.295875806489531	0.0	2905640.0	0.0	43.55016700923443	29.428849400342717	1.903743377653055	147.99168233119116	442.87693678889445	458.14303178340197	468168.6393911941	108.32414958626032	741.2191173434258	199.61485552173647	-0.43557073499837085	93.58211707175519	663.6434283040721	654.591298699379	689503.5873361042	459.91496443748474	911.4427809789777	132.60614885246827	2.0117351281983806	76.12668787371726	923.2319039154182	939.9905848503113	638092.7737545212	724.9960894882679	1062.2933142632246	103.60824329073633	0.2406811374858812	0.5351042169358428	0.549711129303519	0.405317758120108
+sub-50152_T1w	0.5389448521321344	1.951832110632953	0.6407575000929345	6400.0	3.538223004967018	3.5068090149010525	3.72369	3.38417	0.24581070515371412	0.39994911569375136	0.3542401791525344	0.6055886149406433	0.29949661493301377	0.0	0.007373185025089607	5.3317928277705136	5.260562667625682	5.373010998628089	160	239	200	2.859771605234868	7.654488968704144	6.401109684774777	8.68906848038532	25.078592655941325	31.375238888142178	33.547277021352265	44.18799951997329	1.100000023841858	1.0	1.0	6.414540030448748	0.0	7.886303632707387	0.0	2096437.0	0.0	38.88384077697992	14.240638115831892	0.04184774685015302	176.6526273457668	525.277949354461	545.1310024410486	343748.8186203772	191.50386236608028	807.6631854474545	190.62019096981507	-0.5752850710339215	70.2793420517636	686.4136575956873	682.0006075128913	559300.441947877	546.0775418952107	840.110557936132	89.09804437765503	-0.3790885336436278	79.82206421534173	933.7240992117161	960.5103765055537	495379.73952520755	712.4791779369116	1076.101774647832	110.5422761029612	0.24458913017627987	0.5240078354587275	0.534175557705883	0.47523509797309954
+sub-50785_T1w	0.7462348936737846	1.3384000692729365	0.7530480076783985	-1.0	3.499472057939973	3.46213	3.7065461738199192	3.32974	0.26652116229801986	0.3772997886870599	0.35617904901492026	0.9645677804946899	0.27349820435047145	0.0	0.0009944221999575347	7.277619985914051	7.159163031906821	7.30470686003244	256	180	256	2.3315218639714548	3.6975246632030516	3.8881798392340765	5.635492990527722	6.229615920210668	8.968627139000166	9.61189448660999	13.637440400619134	1.0	0.9999875426292419	1.0	579.384989482046	0.0	2.7983003899529795	0.0	5091436.0	0.0	0.0	41.35426782549021	0.42896562817136985	119.23098496445634	386.216226576264	393.23292861506343	422567.70370231196	97.09253230690956	657.3769222311676	168.65913607823262	-0.21950911347902347	105.7205633958213	571.9650741568756	566.1279219612479	598206.5511727771	329.8122172765434	835.4616258479655	153.109850600149	0.8381371141597027	114.20244240752042	846.6359659267068	860.8380831554532	564719.74514658	578.6322632431984	1077.089632384479	152.75279774464403	0.24122428647139732	0.5360100491712251	0.5573704873348334	0.34931847577355485
+sub-51187_T1w	0.7051518353945694	1.215726053025216	0.6068733970705135	4668.127902703064	3.5394978356516567	4.25443345257453	2.965211919410108	3.398848134970332	0.2607189196731251	0.3850319288883145	0.35424915143856045	0.8645452558994293	0.659174335002899	0.0	0.002621822353053779	6.231531509930732	6.101578640123206	6.23860156813424	256	132	256	1.9068147296421873	3.490092089821521	3.4894459153244077	5.071430926509515	5.401224044545881	7.612318350917183	8.178972191748775	11.523374179783266	0.8593999743461609	1.5000007152557373	0.8593999743461609	643.420315158402	0.0	12.502871617173104	0.0	2453296.0	0.0	63.80114300921559	49.546339902331766	17.383575114497404	203.01221375156962	394.8127140206175	408.4811688065529	352493.0222813282	37.886436089873314	722.4866435565054	214.2214359050029	19.093366757493296	103.00303682268894	585.9201437803881	575.7007433250546	520564.7079959175	332.4198212251067	879.8980544172227	164.95272203442693	10.128701373968426	105.56920244595317	841.1839217378181	871.4841885343194	478946.2697518447	541.1317743547261	1057.166188761592	171.84169347257236	0.23964798460924375	0.5492138840500247	0.5613043528545612	0.3146424394881272
diff --git a/mriqc/interfaces/anatomical.py b/mriqc/interfaces/anatomical.py
@@ -105,11 +105,6 @@ class StructuralQC(SimpleInterface):
 
     def _run_interface(self, runtime):  # pylint: disable=R0914,E1101
         imnii = nb.load(self.inputs.in_noinu)
-        erode = (
-            np.all(np.array(imnii.header.get_zooms()[:3], dtype=np.float32) < 1.9)
-            if self.inputs.human
-            else False
-        )
 
         # Load image corrected for INU
         inudata = np.nan_to_num(imnii.get_fdata())
@@ -144,18 +139,18 @@ def _run_interface(self, runtime):  # pylint: disable=R0914,E1101
 
         rotdata = np.asanyarray(nb.load(self.inputs.rot_msk).dataobj).astype(np.uint8)
 
-        # Load Partial Volume Maps (pvms) from FSL FAST
-        pvmdata = []
-        for fname in self.inputs.in_pvms:
-            pvmdata.append(nb.load(fname).get_fdata(dtype="float32"))
-            if np.sum(pvmdata[-1] > 1e-4) < 10:
-                raise RuntimeError(
-                    "Detected less than 10 voxels belonging to one tissue prob. map. "
-                    "MRIQC failed to process this dataset."
-                )
+        # Load brain tissue probability maps from GMM segmentation
+        pvms = {
+            label: nb.load(fname).get_fdata()
+            for label, fname in zip(("csf", "gm", "wm"), self.inputs.in_pvms)
+        }
+        pvmdata = list(pvms.values())
+
+        # Add probability maps
+        pvms["bg"] = airdata
 
         # Summary stats
-        stats = summary_stats(inudata, pvmdata, airdata, erode=erode)
+        stats = summary_stats(inudata, pvms)
         self._results["summary"] = stats
 
         # SNR

diff --git a/mriqc/interfaces/functional.py b/mriqc/interfaces/functional.py
@@ -56,13 +56,9 @@ class FunctionalQCInputSpec(BaseInterfaceInputSpec):
         mandatory=True,
         desc="motion parameters for FD computation",
     )
-    fd_thres = traits.Float(
-        0.2, usedefault=True, desc="motion threshold for FD computation"
-    )
+    fd_thres = traits.Float(0.2, usedefault=True, desc="motion threshold for FD computation")
     in_dvars = File(exists=True, mandatory=True, desc="input file containing DVARS")
-    in_fwhm = traits.List(
-        traits.Float, mandatory=True, desc="smoothness estimated with AFNI"
-    )
+    in_fwhm = traits.List(traits.Float, mandatory=True, desc="smoothness estimated with AFNI")
 
 
 class FunctionalQCOutputSpec(TraitedSpec):
@@ -105,23 +101,21 @@ def _run_interface(self, runtime):
         # Get brain mask data
         msknii = nb.load(self.inputs.in_mask)
         mskdata = np.asanyarray(msknii.dataobj) > 0
-        mskdata = mskdata.astype(np.uint8)
         if np.sum(mskdata) < 100:
             raise RuntimeError(
                 "Detected less than 100 voxels belonging to the brain mask. "
                 "MRIQC failed to process this dataset."
             )
 
         # Summary stats
-        stats = summary_stats(epidata, mskdata, erode=True)
+        rois = {"fg": mskdata.astype(np.uint8), "bg": (~mskdata).astype(np.uint8)}
+        stats = summary_stats(epidata, rois)
         self._results["summary"] = stats
 
         # SNR
-        self._results["snr"] = snr(
-            stats["fg"]["median"], stats["fg"]["stdv"], stats["fg"]["n"]
-        )
+        self._results["snr"] = snr(stats["fg"]["median"], stats["fg"]["stdv"], stats["fg"]["n"])
         # FBER
-        self._results["fber"] = fber(epidata, mskdata)
+        self._results["fber"] = fber(epidata, mskdata.astype(np.uint8))
         # EFC
         self._results["efc"] = efc(epidata)
         # GSR
@@ -132,20 +126,18 @@ def _run_interface(self, runtime):
             epidir = [self.inputs.direction]
 
         for axis in epidir:
-            self._results["gsr"][axis] = gsr(epidata, mskdata, direction=axis)
+            self._results["gsr"][axis] = gsr(epidata, mskdata.astype(np.uint8), direction=axis)
 
         # DVARS
-        dvars_avg = np.loadtxt(
-            self.inputs.in_dvars, skiprows=1, usecols=list(range(3))
-        ).mean(axis=0)
+        dvars_avg = np.loadtxt(self.inputs.in_dvars, skiprows=1, usecols=list(range(3))).mean(
+            axis=0
+        )
         dvars_col = ["std", "nstd", "vstd"]
-        self._results["dvars"] = {
-            key: float(val) for key, val in zip(dvars_col, dvars_avg)
-        }
+        self._results["dvars"] = {key: float(val) for key, val in zip(dvars_col, dvars_avg)}
 
         # tSNR
         tsnr_data = nb.load(self.inputs.in_tsnr).get_fdata()
-        self._results["tsnr"] = float(np.median(tsnr_data[mskdata > 0]))
+        self._results["tsnr"] = float(np.median(tsnr_data[mskdata]))
 
         # FD
         fd_data = np.loadtxt(self.inputs.in_fd, skiprows=1)
@@ -157,9 +149,7 @@ def _run_interface(self, runtime):
         }
 
         # FWHM
-        fwhm = np.array(self.inputs.in_fwhm[:3]) / np.array(
-            hmcnii.header.get_zooms()[:3]
-        )
+        fwhm = np.array(self.inputs.in_fwhm[:3]) / np.array(hmcnii.header.get_zooms()[:3])
         self._results["fwhm"] = {
             "x": float(fwhm[0]),
             "y": float(fwhm[1]),
@@ -302,6 +292,4 @@ 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
+    return (data - np.atleast_2d(np.median(data, axis=1)).T) / np.atleast_2d(data.std(axis=1)).T
diff --git a/mriqc/qc/anatomical.py b/mriqc/qc/anatomical.py
@@ -289,9 +289,7 @@ def cnr(mu_wm, mu_gm, sigma_air, sigma_wm, sigma_gm):
     :return: the computed CNR
 
     """
-    return float(
-        abs(mu_wm - mu_gm) / sqrt(sigma_air**2 + sigma_gm**2 + sigma_wm**2)
-    )
+    return float(abs(mu_wm - mu_gm) / sqrt(sigma_air**2 + sigma_gm**2 + sigma_wm**2))
 
 
 def cjv(mu_wm, mu_gm, sigma_wm, sigma_gm):
@@ -392,10 +390,7 @@ def efc(img, framemask=None):
     # Calculate EFC (add 1e-16 to the image data to keep log happy)
     return float(
         (1.0 / efc_max)
-        * np.sum(
-            (img[framemask == 0] / b_max)
-            * np.log((img[framemask == 0] + 1e-16) / b_max)
-        )
+        * np.sum((img[framemask == 0] / b_max) * np.log((img[framemask == 0] + 1e-16) / b_max))
     )
 
 
@@ -557,13 +552,11 @@ def rpve(pvms, seg):
         loth = np.percentile(pvmap[pvmap > 0], 2)
         pvmap[pvmap < loth] = 0
         pvmap[pvmap > upth] = 0
-        pvfs[k] = (
-            pvmap[pvmap > 0.5].sum() + (1.0 - pvmap[pvmap <= 0.5]).sum()
-        ) / totalvol
+        pvfs[k] = (pvmap[pvmap > 0.5].sum() + (1.0 - pvmap[pvmap <= 0.5]).sum()) / totalvol
     return {k: float(v) for k, v in list(pvfs.items())}
 
 
-def summary_stats(img, pvms, airmask=None, erode=True):
+def summary_stats(data, pvms, airmask=None, erode=True):
     r"""
     Estimates the mean, the median, the standard deviation,
     the kurtosis,the median absolute deviation (mad), the 95\%
@@ -584,67 +577,29 @@ def summary_stats(img, pvms, airmask=None, erode=True):
 
 
     """
+    from statsmodels.stats.weightstats import DescrStatsW
     from statsmodels.robust.scale import mad
 
-    from .. import config
-
-    # Check type of input masks
-    dims = np.squeeze(np.array(pvms)).ndim
-    if dims == 4:
-        # If pvms is from FSL FAST, create the bg mask
-        stats_pvms = [np.zeros_like(img)] + pvms
-    elif dims == 3:
-        stats_pvms = [np.ones_like(pvms) - pvms, pvms]
-    else:
-        raise RuntimeError(
-            "Incorrect image dimensions ({0:d})".format(np.array(pvms).ndim)
-        )
-
-    if airmask is not None:
-        stats_pvms[0] = airmask
-
-    labels = list(FSL_FAST_LABELS.items())
-    if len(stats_pvms) == 2:
-        labels = list(zip(["bg", "fg"], list(range(2))))
-
     output = {}
-    for k, lid in labels:
-        mask = stats_pvms[lid] > 0.85
-        nvox = mask.sum()
-
-        if erode and nvox > 2e3:
-            struct = nd.generate_binary_structure(3, 2)
-            mask = nd.binary_erosion(mask, structure=struct)
-            nvox = mask.sum()
-
-        # Initialize with zeros.
-        output[k] = {
-            "mean": 0.0,
-            "median": 0.0,
-            "p95": 0.0,
-            "p05": 0.0,
-            "k": 0.0,
-            "stdv": 0.0,
-            "mad": 0.0,
+    for label, probmap in pvms.items():
+        wstats = DescrStatsW(data=data.reshape(-1), weights=probmap.reshape(-1))
+        nvox = probmap.sum()
+        p05, median, p95 = wstats.quantile(
+            np.array([0.05, 0.50, 0.95]),
+            return_pandas=False,
+        )
+        thresholded = data[probmap > (0.5 * probmap.max())]
+
+        output[label] = {
+            "mean": float(wstats.mean),
+            "median": float(median),
+            "p95": float(p95),
+            "p05": float(p05),
+            "k": float(kurtosis(thresholded)),
+            "stdv": float(wstats.std),
+            "mad": float(mad(thresholded, center=median)),
             "n": float(nvox),
         }
-        if nvox > 100:
-            data = img[mask]
-            output[k] = {
-                "mean": float(data.mean()),
-                "stdv": float(data.std()),
-                "median": float(np.median(data)),
-                "mad": float(mad(data)),
-                "p95": float(np.percentile(data, 95)),
-                "p05": float(np.percentile(data, 5)),
-                "k": float(kurtosis(data)),
-                "n": float(nvox),
-            }
-        else:
-            config.loggers.interface.warning(
-                f'calculating summary stats of label "{k}" in a very small '
-                f"mask ({nvox} voxels)",
-            )
 
     return output