Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

ENH: Drop utilization of "head" mask from template #1104

Merged
merged 7 commits into from
Mar 30, 2023
Merged

ENH: Drop utilization of "head" mask from template #1104

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
eec1622
enh: drop utilization of "head" mask from template
oesteban Mar 29, 2023
4194206
sty: pacify flake8, run black
oesteban Mar 29, 2023
aadefc5
sty: run black
oesteban Mar 29, 2023
2ed3022
fix: standardization of background distribution
oesteban Mar 29, 2023
4059dd6
maint: pin feature branch of nitransforms
oesteban Mar 30, 2023
e9dee0e
fix: update metrics oracle
oesteban Mar 30, 2023
2332c7b
doc: better document the new feature in the air mask workflow
oesteban Mar 30, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions mriqc/data/testdata/group_T1w.tsv
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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.37313632227722987 3.511490654414502 0.6960134252231683 7605.10513033907 4.202703333333333 4.04437 4.57298 3.99076 0.18062183162155526 0.45129552684595725 0.3680826415324875 0.6880311369895935 0.3810786575078964 0.0 0.001930232337521182 18.825247702549973 8.39017106314948 11.649061977664282 208 256 176 1.6993239139994756 9.60548868379397 9.580703779775568 17.43729874153326 14.688511153161976 44.63016146481638 42.10615608188988 66.99979562769128 1.0 1.0 1.0 31.012418369962887 9.778582763117985 18.452746423366403 6.59555384516716 1001410.0 0.0 74.66818365454674 30.701479513750634 19.004358184769764 99.59667629101651 240.17309068576688 219.24109540879726 23305.0 76.80156033039093 460.92335920929895 129.01389182901477 0.027295873315861918 68.99333393977689 666.3344565866706 666.1509383618832 48318.0 551.0137391388416 780.3517317920923 69.35035445376488 0.37152398970452927 55.593053856985584 1001.987130865773 1000.0406734496355 249478.0 911.7335358560085 1099.2182608991861 57.350549761438494 0.16935039516742867 0.46407615668059243 0.5031231842500415 0.41705010441957996
sub-50152_T1w 0.31516709879412075 4.201130171979436 0.6407575000929343 3312.974194636532 3.538223004967018 3.5068090149010525 3.72369 3.38417 0.14842499625730998 0.5250892530167288 0.3264857507259612 0.6055912971496582 0.2994886189699173 0.0 0.006674976684377104 24.0285265653838 7.9353643561617515 13.148145318541953 160 239 200 2.510834587303088 11.943084741826837 11.615337183684636 20.39209222192398 22.324232151069324 46.41076606716801 45.741405159382374 68.48921725990981 1.100000023841858 1.0 1.0 3.606360415844245 9.565905773394173 12.380584157317957 6.452105395495892 1134850.0 0.0 45.89012239873409 16.24558734327338 0.7961230412523719 114.96067825895601 334.6061178320976 325.96494595706463 13917.0 121.97151667177677 557.9582219704985 129.8186811509053 0.08243483568350429 55.301123908911876 678.5507979221064 677.661957219243 56655.0 585.84353428334 774.0235786288977 56.740447818984975 0.5065511281865795 46.301248062937475 1000.1093290204658 1000.0381581634283 147620.0 920.0549581423402 1080.6341173063959 49.040322104779726 0.1462579099851338 0.4817811633325324 0.47714719815037343 0.49068974560713347
sub-50785_T1w 0.3739190780358043 3.279949774261422 0.7530480076783987 -1.0 3.499472057939973 3.46213 3.7065461738199192 3.32974 0.189876038761227 0.44213808013778255 0.36798588110099045 0.9646123945713043 0.2733038604259488 0.0 0.000767147494457557 26.060490432885494 11.501433085474199 15.490402541642467 256 180 256 1.6733744484964148 7.939790311775997 6.544267894218604 10.019638922383399 2.8588400295341465 8.487622615748716 8.734310428505255 14.856468640232904 1.0 0.9999875426292419 1.0 14.337529878058128 0.0 20.32412254580298 0.0 103748.0 0.0 111.94158869981766 44.09772058964646 0.9433749341342765 118.2763764450988 205.39404429515787 192.4306584596634 20108.0 39.272444665431976 405.72885352373123 114.99271645101418 0.1807781348962867 69.97938871717868 571.6115923516921 571.3082200586796 21002.0 452.9962001532316 692.0828685075045 71.95336352499645 1.2097531340465286 90.31664813841826 1009.8780921704935 1000.0000046491623 178124.0 863.9937826395035 1189.4063824415207 99.80371601794556 0.1668555373238035 0.4587601573604803 0.4864809220911667 0.38380923954498725
sub-51187_T1w 0.33882862703385835 3.9544575680768377 0.6068733970705135 4679.158051897638 3.5394978356516567 4.25443345257453 2.965211919410108 3.398848134970332 0.16404671986716 0.4935345249617393 0.3424187551711007 0.8655245900154114 0.6595467209815977 0.0 0.00011954362363694744 26.046560208106577 9.5216518888818 14.574225774565697 256 132 256 0.5218688172129379 7.963006853091069 7.5318240199687905 14.110596389602364 6.657180161772777 14.980777564628841 15.760659092162987 25.644019550087343 0.8593999743461609 1.5000007152557373 0.8593999743461609 9.847905534202837 0.0 16.162994984246847 0.0 913891.0 0.0 67.68180740252137 25.547336056604777 3.8791461166165755 170.97945856656423 434.6798289708983 259.59970897994936 5336.0 45.10847321804613 1694.362357551232 497.3958473249557 0.004220801228355775 73.69218439433283 582.6507497875704 584.1821013651788 47874.0 459.6302606094629 698.9381590856239 73.36123286903359 0.6090564532154255 67.19881568745514 1002.0922130810112 999.9999775439501 138601.0 889.7433086000383 1121.053142476827 70.86846951393275 0.1331905213170387 0.4737172245107991 0.4812258086341682 0.340853963466229
sub-50137_T1w 0.37313632227722987 3.5844663252515097 0.6960134252231683 7605.10513033907 4.202703333333333 4.04437 4.57298 3.99076 0.18062183162155526 0.45129552684595725 0.3680826415324875 0.6880311369895935 0.3810786575078964 0.0 0.0016246831072444018 18.825247702549973 8.39017106314948 11.649061977664282 208 256 176 1.6993239139994756 9.60548868379397 9.580703779775568 17.43729874153326 5.973337476873543 18.14962818878457 17.123197682080704 27.246627380583995 1.0 1.0 1.0 367.95295392615037 0.0 7.789915131389547 0.0 2969163.0 0.0 48.733814522624016 24.04565664241587 19.004358184769764 99.59667629101651 240.17309068576688 219.24109540879726 23305.0 76.80156033039093 460.92335920929895 129.01389182901477 0.027295873315861918 68.99333393977689 666.3344565866706 666.1509383618832 48318.0 551.0137391388416 780.3517317920923 69.35035445376488 0.37152398970452927 55.593053856985584 1001.987130865773 1000.0406734496355 249478.0 911.7335358560085 1099.2182608991861 57.350549761438494 0.1707051203342759 0.4671890606536107 0.5031855306067675 0.41705010441957996
sub-50152_T1w 0.31516709879412075 4.172440280943563 0.6407575000929343 3312.974194636532 3.538223004967018 3.5068090149010525 3.72369 3.38417 0.14842499625730998 0.5250892530167288 0.3264857507259612 0.6055912971496582 0.2994886189699173 0.0 0.00504991530578035 24.0285265653838 7.9353643561617515 13.148145318541953 160 239 200 2.510834587303088 11.943084741826837 11.615337183684636 20.39209222192398 11.494317323697395 23.896009896419393 23.551368852279484 35.26377933672167 1.100000023841858 1.0 1.0 630.6396765272589 0.0 9.778584125865612 0.0 2189774.0 0.0 45.05020335316658 18.57887643142015 0.7961230412523719 114.96067825895601 334.6061178320976 325.96494595706463 13917.0 121.97151667177677 557.9582219704985 129.8186811509053 0.08243483568350429 55.301123908911876 678.5507979221064 677.661957219243 56655.0 585.84353428334 774.0235786288977 56.740447818984975 0.5065511281865795 46.301248062937475 1000.1093290204658 1000.0381581634283 147620.0 920.0549581423402 1080.6341173063959 49.040322104779726 0.14675306814451586 0.4834100383925821 0.47801106843313007 0.49068974560713347
sub-50785_T1w 0.3739190780358043 3.460710620103788 0.7530480076783987 -1.0 3.499472057939973 3.46213 3.7065461738199192 3.32974 0.189876038761227 0.44213808013778255 0.36798588110099045 0.9646123945713043 0.2733038604259488 0.0 0.0018900993982801935 26.060490432885494 11.501433085474199 15.490402541642467 256 180 256 1.6733744484964148 7.939790311775997 6.544267894218604 10.019638922383399 8.769148428721659 26.034762965244912 26.791448202358822 45.5704332131099 1.0 0.9999875426292419 1.0 226.8740565908444 0.0 1.9238560057587766 0.0 5103243.0 0.0 0.0 14.376347926781678 0.9433749341342765 118.2763764450988 205.39404429515787 192.4306584596634 20108.0 39.272444665431976 405.72885352373123 114.99271645101418 0.1807781348962867 69.97938871717868 571.6115923516921 571.3082200586796 21002.0 452.9962001532316 692.0828685075045 71.95336352499645 1.2097531340465286 90.31664813841826 1009.8780921704935 1000.0000046491623 178124.0 863.9937826395035 1189.4063824415207 99.80371601794556 0.1680471213771674 0.4634107232994584 0.4858448244116643 0.38380923954498725
sub-51187_T1w 0.33882862703385835 3.9445394143408588 0.6068733970705135 4679.158051897638 3.5394978356516567 4.25443345257453 2.965211919410108 3.398848134970332 0.16404671986716 0.4935345249617393 0.3424187551711007 0.8655245900154114 0.6595467209815977 0.0 6.414677477794047e-05 26.046560208106577 9.5216518888818 14.574225774565697 256 132 256 0.5218688172129379 7.963006853091069 7.5318240199687905 14.110596389602364 6.3902127171011305 14.38001570924317 15.128622286499123 24.61563843315307 0.8593999743461609 1.5000007152557373 0.8593999743461609 20.953484796959334 0.0 11.339826321213017 0.0 2337817.0 0.0 64.88082966953516 26.61464121953738 3.8791461166165755 170.97945856656423 434.6798289708983 259.59970897994936 5336.0 45.10847321804613 1694.362357551232 497.3958473249557 0.004220801228355775 73.69218439433283 582.6507497875704 584.1821013651788 47874.0 459.6302606094629 698.9381590856239 73.36123286903359 0.6090564532154255 67.19881568745514 1002.0922130810112 999.9999775439501 138601.0 889.7433086000383 1121.053142476827 70.86846951393275 0.13692615298141483 0.47923963498730343 0.4788423111786441 0.340853963466229
144 changes: 63 additions & 81 deletions mriqc/interfaces/anatomical.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,8 +21,7 @@
# https://www.nipreps.org/community/licensing/
#
"""Nipype interfaces to support anatomical workflow."""
import os.path as op
from builtins import zip
from pathlib import Path

import nibabel as nb
import numpy as np
Expand All @@ -48,7 +47,6 @@
InputMultiPath,
SimpleInterface,
TraitedSpec,
isdefined,
traits,
)
from nipype.utils.filemanip import fname_presuffix
Expand All @@ -70,9 +68,7 @@ class StructuralQCInputSpec(BaseInterfaceInputSpec):
)
in_tpms = InputMultiPath(File(), desc="tissue probability maps from FSL FAST")
mni_tpms = InputMultiPath(File(), desc="tissue probability maps from FSL FAST")
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")
human = traits.Bool(True, usedefault=True, desc="human workflow")


Expand Down Expand Up @@ -137,9 +133,7 @@ def _run_interface(self, runtime): # pylint: disable=R0914,E1101

# Load air, artifacts and head masks
airdata = np.asanyarray(nb.load(self.inputs.air_msk).dataobj).astype(np.uint8)
artdata = np.asanyarray(nb.load(self.inputs.artifact_msk).dataobj).astype(
np.uint8
)
artdata = np.asanyarray(nb.load(self.inputs.artifact_msk).dataobj).astype(np.uint8)

headdata = np.asanyarray(nb.load(self.inputs.head_msk).dataobj).astype(np.uint8)
if np.sum(headdata > 0) < 100:
Expand Down Expand Up @@ -222,9 +216,7 @@ def _run_interface(self, runtime): # pylint: disable=R0914,E1101
)

# FWHM
fwhm = np.array(self.inputs.in_fwhm[:3]) / np.array(
imnii.header.get_zooms()[:3]
)
fwhm = np.array(self.inputs.in_fwhm[:3]) / np.array(imnii.header.get_zooms()[:3])
self._results["fwhm"] = {
"x": float(fwhm[0]),
"y": float(fwhm[1]),
Expand Down Expand Up @@ -261,9 +253,7 @@ def _run_interface(self, runtime): # pylint: disable=R0914,E1101
# Bias
bias = nb.load(self.inputs.in_bias).get_fdata()[segdata > 0]
self._results["inu"] = {
"range": float(
np.abs(np.percentile(bias, 95.0) - np.percentile(bias, 5.0))
),
"range": float(np.abs(np.percentile(bias, 95.0) - np.percentile(bias, 5.0))),
"med": float(np.median(bias)),
} # pylint: disable=E1101

Expand All @@ -281,18 +271,26 @@ def _run_interface(self, runtime): # pylint: disable=R0914,E1101
return runtime


class ArtifactMaskInputSpec(BaseInterfaceInputSpec):
class _ArtifactMaskInputSpec(BaseInterfaceInputSpec):
in_file = File(exists=True, mandatory=True, desc="File to be plotted")
head_mask = File(exists=True, mandatory=True, desc="head mask")
rot_mask = File(exists=True, desc="a rotation mask")
nasion_post_mask = File(
exists=True,
mandatory=True,
desc="nasion to posterior of cerebellum mask",
glabella_xyz = traits.Tuple(
(0.0, 90.0, -14.0),
types=(traits.Float, traits.Float, traits.Float),
usedefault=True,
desc="position of the top of the glabella in standard coordinates",
)
inion_xyz = traits.Tuple(
(0.0, -120.0, -14.0),
types=(traits.Float, traits.Float, traits.Float),
usedefault=True,
desc="position of the top of the inion in standard coordinates",
)
ind2std_xfm = File(exists=True, mandatory=True, desc="individual to standard affine transform")
zscore = traits.Float(10.0, usedefault=True, desc="z-score to consider artifacts")


class ArtifactMaskOutputSpec(TraitedSpec):
class _ArtifactMaskOutputSpec(TraitedSpec):
out_hat_msk = File(exists=True, desc='output "hat" mask')
out_art_msk = File(exists=True, desc="output artifacts mask")
out_air_msk = File(exists=True, desc='output "hat" mask, without artifacts')
Expand All @@ -303,60 +301,57 @@ class ArtifactMask(SimpleInterface):
Computes the artifact mask using the method described in [Mortamet2009]_.
"""

input_spec = ArtifactMaskInputSpec
output_spec = ArtifactMaskOutputSpec
input_spec = _ArtifactMaskInputSpec
output_spec = _ArtifactMaskOutputSpec

def _run_interface(self, runtime):
imnii = nb.load(self.inputs.in_file)
from nibabel.affines import apply_affine
from nitransforms.linear import Affine

in_file = Path(self.inputs.in_file)
imnii = nb.as_closest_canonical(nb.load(in_file))
imdata = np.nan_to_num(imnii.get_fdata().astype(np.float32))

# Remove negative values
imdata[imdata < 0] = 0
xfm = Affine.from_filename(self.inputs.ind2std_xfm, fmt="itk")

hmdata = np.asanyarray(nb.load(self.inputs.head_mask).dataobj)
npdata = np.asanyarray(nb.load(self.inputs.nasion_post_mask).dataobj)
ras2ijk = np.linalg.inv(imnii.affine)
glabella_ijk, inion_ijk = apply_affine(
ras2ijk, xfm.map([self.inputs.glabella_xyz, self.inputs.inion_xyz])
)

# Invert head mask
airdata = np.ones_like(hmdata, dtype=np.uint8)
airdata[hmdata == 1] = 0
hmdata = np.bool_(nb.load(self.inputs.head_mask).dataobj)

# Calculate distance to border
dist = nd.morphology.distance_transform_edt(airdata)
dist = nd.morphology.distance_transform_edt(~hmdata)

# Apply nasion-to-posterior mask
airdata[npdata == 1] = 0
dist[npdata == 1] = 0
dist /= dist.max()
hmdata[:, :, : int(inion_ijk[2])] = 1
hmdata[:, (hmdata.shape[1] // 2) :, : int(glabella_ijk[2])] = 1

# Apply rotation mask (if supplied)
if isdefined(self.inputs.rot_mask):
rotmskdata = np.asanyarray(nb.load(self.inputs.rot_mask).dataobj)
airdata[rotmskdata == 1] = 0
dist[~hmdata] = 0
dist /= dist.max()

# Run the artifact detection
qi1_img = artifact_mask(imdata, airdata, dist)
qi1_img = artifact_mask(imdata, (~hmdata), dist, zscore=self.inputs.zscore)

fname, ext = op.splitext(op.basename(self.inputs.in_file))
if ext == ".gz":
fname, ext2 = op.splitext(fname)
ext = ext2 + ext
fname = in_file.relative_to(in_file.parent).stem
ext = "".join(in_file.suffixes)

self._results["out_hat_msk"] = op.abspath("{}_hat{}".format(fname, ext))
self._results["out_art_msk"] = op.abspath("{}_art{}".format(fname, ext))
self._results["out_air_msk"] = op.abspath("{}_air{}".format(fname, ext))
outdir = Path(runtime.cwd).absolute()
self._results["out_hat_msk"] = str(outdir / f"{fname}_hat{ext}")
self._results["out_art_msk"] = str(outdir / f"{fname}_art{ext}")
self._results["out_air_msk"] = str(outdir / f"{fname}_air{ext}")

hdr = imnii.header.copy()
hdr.set_data_dtype(np.uint8)
nb.Nifti1Image(qi1_img, imnii.affine, hdr).to_filename(
imnii.__class__(qi1_img.astype(np.uint8), imnii.affine, hdr).to_filename(
self._results["out_art_msk"]
)

nb.Nifti1Image(airdata, imnii.affine, hdr).to_filename(
self._results["out_hat_msk"]
)
airdata = (~hmdata).astype(np.uint8)
imnii.__class__(airdata, imnii.affine, hdr).to_filename(self._results["out_hat_msk"])

airdata[qi1_img > 0] = 0
nb.Nifti1Image(airdata, imnii.affine, hdr).to_filename(
imnii.__class__(airdata.astype(np.uint8), imnii.affine, hdr).to_filename(
self._results["out_air_msk"]
)
return runtime
Expand Down Expand Up @@ -390,9 +385,7 @@ def _run_interface(self, runtime):


class HarmonizeInputSpec(BaseInterfaceInputSpec):
in_file = File(
exists=True, mandatory=True, desc="input data (after bias correction)"
)
in_file = File(exists=True, mandatory=True, desc="input data (after bias correction)")
wm_mask = File(exists=True, mandatory=True, desc="white-matter mask")
erodemsk = traits.Bool(True, usedefault=True, desc="erode mask")
thresh = traits.Float(0.9, usedefault=True, desc="WM probability threshold")
Expand Down Expand Up @@ -427,9 +420,7 @@ def _run_interface(self, runtime):
data = in_file.get_fdata()
data *= 1000.0 / np.median(data[wm_mask > 0])

out_file = fname_presuffix(
self.inputs.in_file, suffix="_harmonized", newpath="."
)
out_file = fname_presuffix(self.inputs.in_file, suffix="_harmonized", newpath=".")
in_file.__class__(data, in_file.affine, in_file.header).to_filename(out_file)

self._results["out_file"] = out_file
Expand Down Expand Up @@ -490,36 +481,27 @@ def _run_interface(self, runtime):


def artifact_mask(imdata, airdata, distance, zscore=10.0):
"""Computes a mask of artifacts found in the air region"""
"""Compute a mask of artifacts found in the air region."""
from statsmodels.robust.scale import mad

if not np.issubdtype(airdata.dtype, np.integer):
airdata[airdata < 0.95] = 0
airdata[airdata > 0.0] = 1

bg_img = imdata * airdata
if np.sum((bg_img > 0).astype(np.uint8)) < 100:
return np.zeros_like(airdata)
qi1_msk = np.zeros(imdata.shape, dtype=bool)
bg_data = imdata[airdata]
if (bg_data > 0).sum() < 10:
return qi1_msk

# Find the background threshold (the most frequently occurring value
# excluding 0)
bg_location = np.median(bg_img[bg_img > 0])
bg_spread = mad(bg_img[bg_img > 0])
bg_img[bg_img > 0] -= bg_location
bg_img[bg_img > 0] /= bg_spread
# Standardize the distribution of the background
bg_spread = mad(bg_data[bg_data > 0])
bg_data[bg_data > 0] = bg_data[bg_data > 0] / bg_spread

# Apply this threshold to the background voxels to identify voxels
# contributing artifacts.
qi1_img = np.zeros_like(bg_img)
qi1_img[bg_img > zscore] = 1
qi1_img[distance < 0.10] = 0
qi1_msk[airdata] = bg_data > zscore
qi1_msk[distance < 0.10] = False

# Create a structural element to be used in an opening operation.
struct = nd.generate_binary_structure(3, 1)
qi1_img = nd.binary_opening(qi1_img, struct).astype(np.uint8)
qi1_img[airdata <= 0] = 0

return qi1_img
qi1_msk = nd.binary_opening(qi1_msk, struct).astype(np.uint8)
return qi1_msk


def fuzzy_jaccard(in_tpms, in_mni_tpms):
Expand Down
Loading