Add --masktype option to control adaptive mask method

tedana/tests/test_decay.py

@@ -16,7 +16,7 @@ def testdata1():
     tes = np.array([14.5, 38.5, 62.5])
     in_files = [op.join(get_test_data_path(), f"echo{i + 1}.nii.gz") for i in range(3)]
     data, _ = io.load_data(in_files, n_echos=len(tes))
-    mask, adaptive_mask = utils.make_adaptive_mask(data, getsum=True)
+    mask, adaptive_mask = utils.make_adaptive_mask(data)
     fittype = "loglin"
     data_dict = {
         "data": data,

tedana/tests/test_metrics.py

@@ -17,7 +17,7 @@ def testdata1():
     tes = np.array([14.5, 38.5, 62.5])
     in_files = [op.join(get_test_data_path(), f"echo{i + 1}.nii.gz") for i in range(3)]
     data_cat, ref_img = io.load_data(in_files, n_echos=len(tes))
-    _, adaptive_mask = utils.make_adaptive_mask(data_cat, getsum=True)
+    _, adaptive_mask = utils.make_adaptive_mask(data_cat)
     data_optcom = np.mean(data_cat, axis=1)
     mixing = np.random.random((data_optcom.shape[1], 50))
     io_generator = io.OutputGenerator(ref_img)

tedana/tests/test_utils.py

@@ -77,10 +77,8 @@ def test_reshape_niimg():
 def test_make_adaptive_mask():
     # load data make masks
     data = io.load_data(fnames, n_echos=len(tes))[0]
-    mask, masksum = utils.make_adaptive_mask(data, getsum=True, threshold=1)
+    mask, masksum = utils.make_adaptive_mask(data, threshold=1)
 
-    # getsum doesn't change mask values
-    assert np.allclose(mask, utils.make_adaptive_mask(data))
     # shapes are all the same
     assert mask.shape == masksum.shape == (64350,)
     assert np.allclose(mask, (masksum >= 1).astype(bool))
@@ -89,13 +87,11 @@ def test_make_adaptive_mask():
     # masksum has correct values
     vals, counts = np.unique(masksum, return_counts=True)
     assert np.allclose(vals, np.array([0, 1, 2, 3]))
-    assert np.allclose(counts, np.array([13564, 3977, 5060, 41749]))
+    assert np.allclose(counts, np.array([13564, 4959, 5349, 40478]))
 
     # test user-defined mask
     # TODO: Add mask file with no bad voxels to test against
-    mask, masksum = utils.make_adaptive_mask(
-        data, mask=pjoin(datadir, "mask.nii.gz"), getsum=True, threshold=3
-    )
+    mask, masksum = utils.make_adaptive_mask(data, mask=pjoin(datadir, "mask.nii.gz"), threshold=3)
     assert np.allclose(mask, (masksum >= 3).astype(bool))
 
 
@@ -127,7 +123,7 @@ def test_smoke_make_adaptive_mask():
 
     in the correct format.
 
-    Note: make_adaptive_mask has optional paramters - mask and getsum.
+    Note: make_adaptive_mask has optional paramters - mask and threshold.
     """
     n_samples = 100
     n_echos = 5
@@ -137,7 +133,6 @@ def test_smoke_make_adaptive_mask():
 
     assert utils.make_adaptive_mask(data) is not None
     assert utils.make_adaptive_mask(data, mask=mask) is not None  # functions with mask
-    assert utils.make_adaptive_mask(data, getsum=True) is not None  # functions when getsumis true
 
 
 def test_smoke_unmask():

tedana/utils.py

@@ -49,37 +49,70 @@ def reshape_niimg(data):
     return fdata
 
 
-def make_adaptive_mask(data, mask=None, getsum=False, threshold=1):
-    """
-    Make map of `data` specifying longest echo a voxel can be sampled with.
+def make_adaptive_mask(data, mask=None, threshold=1):
+    """Make map of `data` specifying longest echo a voxel can be sampled with.
 
     Parameters
     ----------
     data : (S x E x T) array_like
-        Multi-echo data array, where `S` is samples, `E` is echos, and `T` is
-        time
+        Multi-echo data array, where `S` is samples, `E` is echos, and `T` is time.
     mask : :obj:`str` or img_like, optional
         Binary mask for voxels to consider in TE Dependent ANAlysis. Default is
         to generate mask from data with good signal across echoes
-    getsum : :obj:`bool`, optional
-        Return `masksum` in addition to `mask`. Default: False
     threshold : :obj:`int`, optional
         Minimum echo count to retain in the mask. Default is 1, which is
         equivalent not thresholding.
 
     Returns
     -------
     mask : (S,) :obj:`numpy.ndarray`
-        Boolean array of voxels that have sufficient signal in at least one
-        echo
+        Boolean array of voxels that have sufficient signal in at least ``threshold`` echos.
     masksum : (S,) :obj:`numpy.ndarray`
-        Valued array indicating the number of echos with sufficient signal in a
-        given voxel. Only returned if `getsum = True`
+        Valued array indicating the number of echos with sufficient signal in a given voxel.
+
+    Notes
+    -----
+    The adaptive mask is constructed from two methods:
+
+    1.  Count the total number of echoes in each voxel that have "good" data.
+        This method assumes that an exemplar voxel's signal at later echoes is also reasonable,
+        and that any voxels whose values at a given echo are less than 1/3 of the exemplar voxel's
+        values at that echo are affected by dropout.
+
+        This method uses distributions of values across the mask.
+        Therefore, it is sensitive to the quality of the mask;
+        a bad mask may result in a bad adaptive mask.
+
+        This method is implemented as follows:
+
+        a.  Calculate the 33rd percentile of values in the first echo,
+            based on voxel-wise mean over time.
+        b.  Identify the voxel where the first echo's mean value is equal to the 33rd percentile.
+            Basically, this identifies "exemplar" voxel reflecting the 33rd percentile.
+
+            -   The 33rd percentile is arbitrary.
+            -   If more than one voxel has a value exactly equal to the 33rd percentile,
+                keep all of them.
+        c.  Calculate 1/3 of the mean value of the exemplar voxel for each echo.
+
+            -   This is the threshold for "good" data.
+            -   The 1/3 value is arbitrary.
+            -   If there was more than one exemplar voxel,
+                retain the the highest value for each echo.
+        d.  For each voxel, count the number of echoes that have a mean value greater than the
+            corresponding echo's threshold.
+    2.  Determine the echo at which the signal stops decreasing for each voxel.
+        If a voxel's signal stops decreasing as echo time increases, then we can infer that the
+        voxel has either fully dephased (i.e., "bottomed out") or been contaminated by noise.
+        This essentially identifies the last echo with "good" data.
+
+    The element-wise minimum value between the two methods is used to construct the adaptive mask.
     """
     RepLGR.info(
         "An adaptive mask was then generated, in which each voxel's "
         "value reflects the number of echoes with 'good' data."
     )
+    n_samples, n_echos, _ = data.shape
 
     # take temporal mean of echos and extract non-zero values in first echo
     echo_means = data.mean(axis=-1)  # temporal mean of echos
@@ -108,6 +141,15 @@ def make_adaptive_mask(data, mask=None, getsum=False, threshold=1):
     # and count # of echos that pass criterion
     masksum = (np.abs(echo_means) > lthrs).sum(axis=-1)
 
+    # Determine where voxels stop decreasing in signal from echo to echo
+    echo_diffs = np.hstack((np.full((n_samples, 1), -1), np.diff(echo_means, axis=1)))
+    diff_mask = echo_diffs >= 0  # flag where signal is not decreasing
+    last_decreasing_echo = diff_mask.argmax(axis=1)
+    last_decreasing_echo[last_decreasing_echo == 0] = n_echos  # if no increase, set to n_echos
+
+    # Retain the more conservative of the two adaptive mask estimates
+    masksum = np.minimum(masksum, last_decreasing_echo)
+
     if mask is None:
         # make it a boolean mask to (where we have at least `threshold` echoes with good signal)
         mask = (masksum >= threshold).astype(bool)
@@ -127,10 +169,7 @@ def make_adaptive_mask(data, mask=None, getsum=False, threshold=1):
             masksum[masksum < threshold] = 0
             mask = masksum.astype(bool)
 
-    if getsum:
-        return mask, masksum
-
-    return mask
+    return mask, masksum
 
 
 def unmask(data, mask):

tedana/workflows/t2smap.py

@@ -278,7 +278,7 @@ def t2smap_workflow(
         LGR.info("Computing adaptive mask")
     else:
         LGR.info("Using user-defined mask")
-    mask, masksum = utils.make_adaptive_mask(catd, mask=mask, getsum=True, threshold=1)
+    mask, masksum = utils.make_adaptive_mask(catd, mask=mask, threshold=1)
 
     LGR.info("Computing adaptive T2* map")
     if fitmode == "all":

tedana/workflows/tedana.py

@@ -586,12 +586,7 @@ def tedana_workflow(
         )
 
     # Create an adaptive mask with at least 1 good echo, for denoising
-    mask_denoise, masksum_denoise = utils.make_adaptive_mask(
-        catd,
-        mask=mask,
-        getsum=True,
-        threshold=1,
-    )
+    mask_denoise, masksum_denoise = utils.make_adaptive_mask(catd, mask=mask, threshold=1)
     LGR.debug(f"Retaining {mask_denoise.sum()}/{n_samp} samples for denoising")
     io_generator.save_file(masksum_denoise, "adaptive mask img")