Skip to content

Commit

Permalink
Merge pull request #533 from oesteban/fix/528
Browse files Browse the repository at this point in the history 
[FIX] Datasets where air (hat) mask is empty
  • Loading branch information
@oesteban
oesteban committed Jun 1, 2017
2 parents 91368d1 + 36eae0e commit 8ba2576
Show file tree
Hide file tree
Showing 2 changed files with 67 additions and 16 deletions.
10 changes: 5 additions & 5 deletions mriqc/data/testdata/T1w.csv
View file
@@ -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
73 changes: 62 additions & 11 deletions mriqc/qc/anatomical.py
View file
Expand Up @@ -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

Expand Down Expand Up @@ -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"""
Expand Down Expand Up @@ -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:
Expand All @@ -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
Expand All @@ -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):
Expand Down

0 comments on commit 8ba2576

Please sign in to comment.