[FIX] Datasets where air (hat) mask is empty

mriqc/data/testdata/T1w.csv

@@ -1,5 +1,5 @@
-subject_id,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
-50137,0.38345367902512617,3.3390848252812737,0.7364744617160943,143.3969268798828,3.97176,3.84915,4.30615,3.75998,0.2133978898855149,0.42686969810291575,0.35973241201156936,1.2152518033981323,0.45196907520294216,0.0,0.0015283866670782967,16.44412079704881,8.760761342520016,12.31291744093443,208,256,176,1.229885161866144,9.780251289684715,9.71615429906627,18.13832644564795,7.453227749245331,14.339182824122451,14.269061689796928,21.014774496022998,1.0,1.0,1.0,9.003927826475195,31.175771713256836,57.12042999267578,52.41670608520508,931665.0,12.3324462890625,115.81640624999993,33.820404052734375,51.22421869843451,131.04031372070312,380.7245788574219,354.6744384765625,25681.0,158.8266143798828,552.8338623046875,288.37451171875,0.026355223848861353,68.94453430175781,684.1145629882812,682.354248046875,32575.0,571.0976440429687,801.4391540527344,69.76750946044922,0.39190896169503864,52.8669319152832,999.591064453125,1000.0235595703125,253466.0,908.4280853271484,1089.0757141113281,55.1330680847168,0.23382176458835602,0.514546811580658,0.5550331473350525,0.4346056835375569
-50152,0.3039255032198539,4.138389104114073,0.6948864997316359,141.25128173828125,3.4652784600909166,3.4267453802727488,3.65075,3.31834,0.16868135537307122,0.4885912301853383,0.3427274144415905,0.9292771816253662,0.3650368571281435,0.0,0.00640636674822209,20.768276536956023,8.227114386622855,12.576371650369857,160,239,200,2.9823042042332055,11.791153070535714,11.597888009517725,20.020206753784255,14.23166929691522,24.44159809193212,24.945725357356576,36.163908683222374,1.100000023841858,1.0,1.0,0.890990236303058,18.115840911865234,53.94853973388672,51.86589431762695,1026679.0,27.001422309875487,88.03534851074218,18.83237075805664,7.371940288316537,119.54910278320312,401.8689270019531,393.5343322753906,15271.0,200.7536392211914,618.4881591796875,131.9521484375,0.5707499379506284,52.85783767700195,677.1936645507812,675.8594360351562,44244.0,585.3509674072266,775.0797363281249,57.31855010986328,0.8734180296762428,45.65824890136719,997.28564453125,1000.0049438476562,162045.0,909.9402221679688,1073.4810058593748,49.94962692260742,0.2079515904188156,0.5339961051940918,0.5221186280250549,0.5738063926509048
-50785,0.3881841572763748,3.139893975583801,0.7756479415644941,596.0380859375,3.352164997200236,3.199,3.611644991600707,3.24585,0.20414998904652737,0.4384699518998948,0.35738005905357784,1.3165277242660522,0.1928600192070007,0.0,0.001401353360826457,25.763511237624744,12.151800259177817,17.077029805928472,256,180,256,2.212405051008006,8.842233294683378,7.124761956850802,10.319647524861022,7.88793153943957,18.604621983602517,19.090215448458384,30.778092822333075,1.0,0.9999875426292419,1.0,8.30652321640877,21.286128997802734,35.83394241333008,19.925045013427734,54911.0,0.0,120.24365234375,42.48101043701172,-0.1380615749962093,126.91279602050781,254.24244689941406,256.28790283203125,14368.0,67.84258117675782,435.20629577636714,115.8372802734375,0.11097776186414565,66.68869018554688,604.939453125,604.4854125976562,16249.0,493.00067138671875,717.1786499023436,68.361328125,1.874603456607331,86.84977722167969,1009.2825317382812,1000.015380859375,170016.0,868.9075317382812,1180.7556762695312,96.90373992919922,0.20489205420017242,0.49960723519325256,0.5370237827301025,0.46003933468878877
-51187,0.3422398269200667,3.7929999294559202,0.6250856477873908,357.52484130859375,3.7801451255770178,4.346218421569919,3.1878984798915484,3.8063184752695856,0.18837399720708245,0.4525009740272195,0.359125028765698,1.7580251693725586,0.6602846145629884,0.0,0.00011083487655092947,24.788834614717615,10.992975441201002,16.14208395065447,256,132,256,2.1998750852762163,8.231342517618168,7.723079514521209,12.73802094066924,9.244710429990599,18.467250758809488,19.722918095870515,31.456793098811456,0.8593999743461609,1.5000007152557373,0.8593999743461609,4.483646425309813,20.82635498046875,30.774946212768555,25.417274475097656,890628.0,1.623040348291398,77.399446105957,24.54763412475586,99.15187259411111,73.18208312988281,299.7581787109375,293.8832702636719,4277.0,153.52493896484376,419.7971374511718,133.5752716064453,0.24931478286346787,69.55706787109375,583.7290649414062,587.061767578125,20450.0,461.2132232666016,695.972280883789,71.31855010986328,0.8288354047164588,71.76368713378906,998.3399047851562,999.9907836914062,148399.0,864.3146362304688,1123.3076538085938,78.50414276123047,0.17126674950122833,0.48971468210220337,0.5144017338752747,0.35421248708597003
+subject_id,cjv,cnr,efc,fber,fwhm_avg,fwhm_x,fwhm_y,fwhm_z,icvs_csf,icvs_gm,icvs_wm,inu_med,inu_range,provenance,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
+50137,0.38345367902512617,3.3390848252812737,0.7364744617160943,143.3969268798828,3.97176,3.84915,4.30615,3.75998,0.2133978898855149,0.42686969810291575,0.35973241201156936,1.2152518033981323,0.45196907520294216,"{'md5sum': '57cd35190da3c813a3c3dadccd8a4ad7', 'mriqc_pred': 0, 'settings': {'testing': False}, 'software': 'mriqc', 'version': '4a892e1e481901c378356d934bb2a27108868511', 'warnings': {'large_rot_frame': True, 'small_air_mask': False}}",0.0,0.0015283866670782967,16.44412079704881,8.760761342520016,12.31291744093443,208,256,176,1.229885161866144,9.780251289684715,9.71615429906627,18.13832644564795,7.453227749245331,14.339182824122453,14.269061689796928,21.014774496023,1.0,1.0,1.0,9.003927826475195,31.175771713256836,57.12042999267578,52.41670608520508,931665.0,12.3324462890625,115.81640624999993,33.820404052734375,51.22421869843451,131.04031372070312,380.7245788574219,354.6744384765625,25681.0,158.8266143798828,552.8338623046875,288.37451171875,0.026355223848861353,68.94453430175781,684.1145629882812,682.354248046875,32575.0,571.0976440429687,801.4391540527344,69.76750946044922,0.39190896169503864,52.8669319152832,999.591064453125,1000.0235595703125,253466.0,908.4280853271484,1089.0757141113281,55.1330680847168,0.23382176458835602,0.514546811580658,0.5550331473350525,0.4346056835375569
+50152,0.3039255032198539,4.138389104114073,0.6948864997316359,141.25128173828125,3.4652784600909166,3.4267453802727488,3.65075,3.31834,0.16868135537307122,0.4885912301853383,0.3427274144415905,0.9292771816253662,0.3650368571281435,"{'md5sum': '27439008c5ad203d1fe19deef665bb11', 'mriqc_pred': 0, 'settings': {'testing': False}, 'software': 'mriqc', 'version': '4a892e1e481901c378356d934bb2a27108868511', 'warnings': {'large_rot_frame': True, 'small_air_mask': False}}",0.0,0.00640636674822209,20.768276536956023,8.227114386622855,12.576371650369857,160,239,200,2.9823042042332055,11.791153070535714,11.597888009517725,20.020206753784255,14.23166929691522,24.44159809193212,24.945725357356576,36.163908683222374,1.100000023841858,1.0,1.0,0.890990236303058,18.115840911865234,53.94853973388672,51.86589431762695,1026679.0,27.001422309875487,88.03534851074218,18.83237075805664,7.371940288316537,119.54910278320312,401.8689270019531,393.5343322753906,15271.0,200.7536392211914,618.4881591796875,131.9521484375,0.5707499379506284,52.85783767700195,677.1936645507812,675.8594360351562,44244.0,585.3509674072266,775.0797363281249,57.31855010986328,0.8734180296762428,45.65824890136719,997.28564453125,1000.0049438476562,162045.0,909.9402221679688,1073.4810058593748,49.94962692260742,0.2079515904188156,0.5339961051940918,0.5221186280250549,0.5738063926509048
+50785,0.3881841572763748,3.139893975583801,0.7756479415644941,596.0380859375,3.352164997200236,3.199,3.611644991600707,3.24585,0.20414998904652737,0.4384699518998948,0.35738005905357784,1.3165277242660522,0.1928600192070007,"{'md5sum': 'fd16630ba43611e54db04f8914016ef6', 'mriqc_pred': 0, 'settings': {'testing': False}, 'software': 'mriqc', 'version': '4a892e1e481901c378356d934bb2a27108868511', 'warnings': {'large_rot_frame': True, 'small_air_mask': True}}",0.0,0.001401353360826457,25.763511237624744,12.151800259177817,17.077029805928472,256,180,256,2.212405051008006,8.842233294683378,7.124761956850802,10.319647524861022,7.88793153943957,18.60462198360252,19.09021544845839,30.77809282233308,1.0,0.9999875426292419,1.0,8.30652321640877,21.286128997802734,35.83394241333008,19.925045013427734,54911.0,0.0,120.24365234375,42.48101043701172,-0.1380615749962093,126.91279602050781,254.24244689941406,256.28790283203125,14368.0,67.84258117675782,435.20629577636714,115.8372802734375,0.11097776186414565,66.68869018554688,604.939453125,604.4854125976562,16249.0,493.00067138671875,717.1786499023436,68.361328125,1.874603456607331,86.84977722167969,1009.2825317382812,1000.015380859375,170016.0,868.9075317382812,1180.7556762695312,96.90373992919922,0.20489205420017242,0.49960723519325256,0.5370237827301025,0.46003933468878877
+51187,0.3422398269200667,3.7929999294559202,0.6250856477873908,357.52484130859375,3.7801451255770178,4.346218421569919,3.1878984798915484,3.8063184752695856,0.18837399720708245,0.4525009740272195,0.359125028765698,1.7580251693725586,0.6602846145629884,"{'md5sum': '6db5cf300438e2b6303c2d054370c091', 'mriqc_pred': 0, 'settings': {'testing': False}, 'software': 'mriqc', 'version': '4a892e1e481901c378356d934bb2a27108868511', 'warnings': {'large_rot_frame': True, 'small_air_mask': False}}",0.0,0.00011083487655092947,24.788834614717615,10.992975441201002,16.14208395065447,256,132,256,2.1998750852762163,8.231342517618168,7.723079514521209,12.73802094066924,9.244710429990599,18.467250758809488,19.722918095870515,31.456793098811456,0.8593999743461609,1.5000007152557373,0.8593999743461609,4.483646425309813,20.82635498046875,30.774946212768555,25.417274475097656,890628.0,1.623040348291398,77.399446105957,24.54763412475586,99.15187259411111,73.18208312988281,299.7581787109375,293.8832702636719,4277.0,153.52493896484376,419.7971374511718,133.5752716064453,0.24931478286346787,69.55706787109375,583.7290649414062,587.061767578125,20450.0,461.2132232666016,695.972280883789,71.31855010986328,0.8288354047164588,71.76368713378906,998.3399047851562,999.9907836914062,148399.0,864.3146362304688,1123.3076538085938,78.50414276123047,0.17126674950122833,0.48971468210220337,0.5144017338752747,0.35421248708597003

mriqc/qc/anatomical.py

@@ -176,6 +176,7 @@
 from builtins import zip, range, str, bytes  # pylint: disable=W0622
 from six import string_types
 
+DIETRICH_FACTOR = 1.0 / sqrt(2/(4 - pi))
 FSL_FAST_LABELS = {'csf': 1, 'gm': 2, 'wm': 3, 'bg': 0}
 PY3 = version_info[0] > 2
 
@@ -225,7 +226,12 @@ def snr_dietrich(mu_fg, sigma_air):
     :return: the computed SNR for the foreground segmentation
 
     """
-    return float(mu_fg / (sigma_air * sqrt(2/(4 - pi))))
+    if sigma_air < 1.0:
+        from mriqc import MRIQC_LOG
+        MRIQC_LOG.warn('SNRd - background sigma is too small (%f)', sigma_air)
+        sigma_air += 1.0
+
+    return float(DIETRICH_FACTOR * mu_fg / sigma_air)
 
 def cnr(mu_wm, mu_gm, sigma_air):
     r"""
@@ -506,13 +512,28 @@ def summary_stats(img, pvms, airmask=None, erode=True):
     r"""
     Estimates the mean, the standard deviation, the 95\%
     and the 5\% percentiles of each tissue distribution.
+
+    .. warning ::
+
+        Sometimes (with datasets that have been partially processed), the air
+        mask will be empty. In those cases, the background stats will be zero
+        for the mean, median, percentiles and kurtosis, the sum of voxels in
+        the other remaining labels for ``n``, and finally the MAD and the
+        :math:`\sigma` will be calculated as:
+
+        .. math ::
+
+            \sigma_\text{BG} = \sqrt{\sum \sigma_\text{i}^2}
+
+
     """
+    from mriqc import MRIQC_LOG
 
     # 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.array(pvms).sum(axis=0)] + pvms
+        stats_pvms = [np.zeros_like(img)] + pvms
     elif dims == 3:
         stats_pvms = [np.ones_like(pvms) - pvms, pvms]
     else:
@@ -526,7 +547,7 @@ def summary_stats(img, pvms, airmask=None, erode=True):
     if len(stats_pvms) == 2:
         labels = list(zip(['bg', 'fg'], list(range(2))))
 
-    output = {k: {} for k, _ in labels}
+    output = {}
     for k, lid in labels:
         mask = np.zeros_like(img, dtype=np.uint8)
         mask[stats_pvms[lid] > 0.85] = 1
@@ -536,14 +557,44 @@ def summary_stats(img, pvms, airmask=None, erode=True):
             mask = nd.binary_erosion(
                 mask, structure=struc).astype(np.uint8)
 
-        output[k]['mean'] = float(img[mask == 1].mean())
-        output[k]['stdv'] = float(img[mask == 1].std())
-        output[k]['median'] = float(np.median(img[mask == 1]))
-        output[k]['mad'] = float(mad(img[mask == 1]))
-        output[k]['p95'] = float(np.percentile(img[mask == 1], 95))
-        output[k]['p05'] = float(np.percentile(img[mask == 1], 5))
-        output[k]['k'] = float(kurtosis(img[mask == 1]))
-        output[k]['n'] = float(mask.sum())
+        nvox = float(mask.sum())
+        if nvox < 1e3:
+            MRIQC_LOG.warn('calculating summary stats of label "%s" in a very small '
+                           'mask (%d voxels)', k, int(nvox))
+            if k == 'bg':
+                continue
+
+        output[k] = {
+            'mean': float(img[mask == 1].mean()),
+            'stdv': float(img[mask == 1].std()),
+            'median': float(np.median(img[mask == 1])),
+            'mad': float(mad(img[mask == 1])),
+            'p95': float(np.percentile(img[mask == 1], 95)),
+            'p05': float(np.percentile(img[mask == 1], 5)),
+            'k': float(kurtosis(img[mask == 1])),
+            'n': nvox,
+        }
+
+    if 'bg' not in output:
+        output['bg'] = {
+            'mean': 0.,
+            'median': 0.,
+            'p95': 0.,
+            'p05': 0.,
+            'k': 0.,
+            'stdv': sqrt(sum(val['stdv']**2
+                             for _, val in list(output.items()))),
+            'mad': sqrt(sum(val['mad']**2
+                            for _, val in list(output.items()))),
+            'n': sum(val['n'] for _, val in list(output.items()))
+        }
+
+    if 'bg' in output and output['bg']['mad'] == 0.0 and output['bg']['stdv'] > 1.0:
+        MRIQC_LOG.warn('estimated MAD in the background was too small ('
+                       'MAD=%f)', output['bg']['mad'])
+        output['bg']['mad'] = output['bg']['stdv'] / DIETRICH_FACTOR
+
+
     return output
 
 def _prepare_mask(mask, label, erode=True):