Write out T2* & S0 variance + T2*-S0 covariance

tedana/reporting/data/html/report_body_template.html

@@ -216,6 +216,15 @@ <h2>T2*</h2>
       style="width: 95%; max-width: 1200px; height: auto"
     />
   </div>
+  {% endif %} {% if varianceExists %}
+  <h3>T2* Estimate Variance</h3>
+  <div class="carpet-plots-image">
+    <img
+      id="t2starVariancePlot"
+      src="{{ t2starVariancePlot }}"
+      style="width: 95%; max-width: 1200px; height: auto"
+    />
+  </div>
   {% endif %} {% if s0Exists %}
   <h2>S0</h2>
   <div class="carpet-plots-image">
@@ -232,6 +241,23 @@ <h2>S0</h2>
       style="width: 95%; max-width: 1200px; height: auto"
     />
   </div>
+  {% endif %} {% if varianceExists %}
+  <h3>S0 Estimate Variance</h3>
+  <div class="carpet-plots-image">
+    <img
+      id="s0VariancePlot"
+      src="{{ s0VariancePlot }}"
+      style="width: 95%; max-width: 1200px; height: auto"
+    />
+  </div>
+  <h3>T2* and S0 Estimate Covariance</h3>
+  <div class="carpet-plots-image">
+    <img
+      id="t2sS0CovariancePlot"
+      src="{{ t2sS0CovariancePlot }}"
+      style="width: 95%; max-width: 1200px; height: auto"
+    />
+  </div>
   {% endif %} {% if rmseExists %}
   <h2>
     T2* and S0 model fit (RMSE). (Scaled between 2nd and 98th percentiles)

tedana/reporting/html_report.py

@@ -170,12 +170,18 @@ def _update_template_bokeh(
     t2star_histogram_filename = f"{prefix}t2star_histogram.svg"
     t2star_brain = f"./figures/{t2star_brain_filename}"
     t2star_histogram = f"./figures/{t2star_histogram_filename}"
+    t2star_variance_filename = f"{prefix}stat-variance_desc-t2star_statmap.svg"
+    t2star_variance = f"./figures/{t2star_variance_filename}"
 
     # Check for S0 images
     s0_brain_filename = f"{prefix}s0_brain.svg"
     s0_histogram_filename = f"{prefix}s0_histogram.svg"
     s0_brain = f"./figures/{s0_brain_filename}"
     s0_histogram = f"./figures/{s0_histogram_filename}"
+    s0_variance_filename = f"{prefix}stat-variance_desc-s0_statmap.svg"
+    s0_variance = f"./figures/{s0_variance_filename}"
+    t2s_s0_covariance_filename = f"{prefix}stat-covariance_desc-t2star+s0_statmap.svg"
+    t2s_s0_covariance = f"./figures/{t2s_s0_covariance_filename}"
 
     # Check for RMSE images
     rmse_brain_filename = f"{prefix}rmse_brain.svg"
@@ -205,6 +211,13 @@ def _update_template_bokeh(
     LGR.info(f"S0 files exist: {s0_exists}")
     LGR.info(f"RMSE files exist: {rmse_exists}")
 
+    variance_exists = (
+        t2star_variance_filename in files_in_figures
+        and s0_variance_filename in files_in_figures
+        and t2s_s0_covariance_filename in files_in_figures
+    )
+    LGR.info(f"Variance files exist: {variance_exists}")
+
     # Check for external regressors
     erc_filename = f"{prefix}confound_correlations.svg"
     external_regressors_exist = erc_filename in files_in_figures
@@ -231,9 +244,13 @@ def _update_template_bokeh(
         t2starBrainPlot=t2star_brain,
         t2starHistogram=t2star_histogram,
         t2starExists=t2star_exists,
+        t2starVariancePlot=t2star_variance,
         s0BrainPlot=s0_brain,
         s0Histogram=s0_histogram,
         s0Exists=s0_exists,
+        s0VariancePlot=s0_variance,
+        t2sS0CovariancePlot=t2s_s0_covariance,
+        varianceExists=variance_exists,
         rmseBrainPlot=rmse_brain,
         rmseTimeseries=rmse_timeseries,
         rmseExists=rmse_exists,

tedana/reporting/static_figures.py

@@ -626,6 +626,7 @@ def plot_t2star_and_s0(
             vmax=t2s_p98,
             annotate=False,
             output_file=os.path.join(io_generator.out_dir, "figures", t2star_plot),
+            resampling_interpolation="nearest",
         )
 
     # Only plot S0 map if the file exists
@@ -646,6 +647,7 @@ def plot_t2star_and_s0(
                 vmax=s0_p98,
                 annotate=False,
                 output_file=os.path.join(io_generator.out_dir, "figures", s0_plot),
+                resampling_interpolation="nearest",
             )
 
 
@@ -734,6 +736,7 @@ def plot_rmse(
             vmax=rmse_p98,
             annotate=False,
             output_file=rmse_brain_plot,
+            resampling_interpolation="nearest",
         )
 
 
@@ -852,6 +855,7 @@ def plot_gscontrol(
                 cmap="coolwarm",
                 annotate=False,
                 output_file=os.path.join(io_generator.out_dir, "figures", gsr_plot),
+                resampling_interpolation="nearest",
             )
 
     if "mir" in gscontrol:
@@ -872,6 +876,7 @@ def plot_gscontrol(
                 cmap="coolwarm",
                 annotate=False,
                 output_file=os.path.join(io_generator.out_dir, "figures", mir_plot),
+                resampling_interpolation="nearest",
             )
 
     if "gsr" in gscontrol or "mir" in gscontrol:
@@ -1057,3 +1062,55 @@ def _correlate_dataframes(df1, df2):
     # Convert the result back to a DataFrame with appropriate labels
     correlation_df = pd.DataFrame(cross_corr_matrix, index=df1.columns, columns=df2.columns)
     return correlation_df
+
+
+def plot_decay_variance(
+    *,
+    io_generator: io.OutputGenerator,
+    adaptive_mask: np.ndarray,
+):
+    """Plot the variance of the T2* and S0 estimates.
+
+    Parameters
+    ----------
+    io_generator : :obj:`~tedana.io.OutputGenerator`
+        The output generator for this workflow.
+    adaptive_mask : (S,) :obj:`numpy.ndarray`
+        Array where each value indicates the number of echoes with good signal
+        for that voxel. This mask may be thresholded; for example, with values
+        less than 3 set to 0.
+        For more information on thresholding, see `make_adaptive_mask`.
+    """
+    # Mask that only includes values >=2 (i.e. at least 2 good echoes)
+    mask_img = io.new_nii_like(io_generator.reference_img, (adaptive_mask >= 2).astype(np.int32))
+
+    names = [
+        "stat-variance_desc-t2star_statmap",
+        "stat-variance_desc-s0_statmap",
+        "stat-covariance_desc-t2star+s0_statmap",
+    ]
+    imgs = ["t2star variance img", "s0 variance img", "t2star-s0 covariance img"]
+    for name, img in zip(names, imgs):
+        in_file = io_generator.get_name(img)
+        data = masking.apply_mask(in_file, mask_img)
+        data_p02, data_p98 = np.percentile(data, [2, 98])
+        plot_name = f"{io_generator.prefix}{name}.svg"
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                "ignore",
+                message="A non-diagonal affine.*",
+                category=UserWarning,
+            )
+            plotting.plot_stat_map(
+                in_file,
+                bg_img=None,
+                display_mode="mosaic",
+                symmetric_cbar=False,
+                black_bg=True,
+                cmap="Reds",
+                vmin=data_p02,
+                vmax=data_p98,
+                annotate=False,
+                output_file=os.path.join(io_generator.out_dir, "figures", plot_name),
+                resampling_interpolation="nearest",
+            )

tedana/resources/config/outputs.json

@@ -15,6 +15,18 @@
         "orig": "fit_failures",
         "bidsv1.5.0": "desc-decayFitFailures_mask"
     },
+    "t2star variance img": {
+        "orig": "t2s_var",
+        "bidsv1.5.0": "stat-variance_desc-t2star_statmap"
+    },
+    "s0 variance img": {
+        "orig": "s0_var",
+        "bidsv1.5.0": "stat-variance_desc-s0_statmap"
+    },
+    "t2star-s0 covariance img": {
+        "orig": "t2s_s0_covar",
+        "bidsv1.5.0": "stat-covariance_desc-t2star+s0_statmap"
+    },
     "combined img": {
         "orig": "ts_OC",
         "bidsv1.5.0": "desc-optcom_bold"

tedana/tests/data/nih_five_echo_outputs_verbose.txt

@@ -83,6 +83,9 @@ sub-01_echo-5_desc-PCA_components.nii.gz
 sub-01_echo-5_desc-Rejected_bold.nii.gz
 sub-01_references.bib
 sub-01_report.txt
+sub-01_stat-covariance_desc-t2star+s0_statmap.nii.gz
+sub-01_stat-variance_desc-s0_statmap.nii.gz
+sub-01_stat-variance_desc-t2star_statmap.nii.gz
 sub-01_tedana_report.html
 figures
 figures/sub-01_adaptive_mask.svg
@@ -129,5 +132,8 @@ figures/sub-01_rmse_brain.svg
 figures/sub-01_rmse_timeseries.svg
 figures/sub-01_s0_brain.svg
 figures/sub-01_s0_histogram.svg
+figures/sub-01_stat-covariance_desc-t2star+s0_statmap.svg
+figures/sub-01_stat-variance_desc-s0_statmap.svg
+figures/sub-01_stat-variance_desc-t2star_statmap.svg
 figures/sub-01_t2star_brain.svg
 figures/sub-01_t2star_histogram.svg

tedana/tests/test_decay.py

@@ -35,32 +35,38 @@ def testdata1():
 
 def test_fit_decay(testdata1):
     """Fit_decay should return data in (samples,) shape."""
-    t2sv, s0v, t2svg, s0vg, _ = me.fit_decay(
-        testdata1["data"],
+    masked_data = testdata1["data"][testdata1["mask"], ...]
+    masked_adaptive_mask = testdata1["adaptive_mask"][testdata1["mask"]]
+    t2s, s0, failures, t2s_var, s0_var, t2s_s0_covar = me.fit_decay(
+        masked_data,
         testdata1["tes"],
-        testdata1["mask"],
-        testdata1["adaptive_mask"],
+        masked_adaptive_mask,
         testdata1["fittype"],
     )
-    assert t2sv.ndim == 1
-    assert s0v.ndim == 1
-    assert t2svg.ndim == 1
-    assert s0vg.ndim == 1
+    assert t2s.ndim == 1
+    assert s0.ndim == 1
+    assert failures is None
+    assert t2s_var is None
+    assert s0_var is None
+    assert t2s_s0_covar is None
 
 
 def test_fit_decay_ts(testdata1):
     """Fit_decay_ts should return data in samples x time shape."""
-    t2sv, s0v, t2svg, s0vg, _ = me.fit_decay_ts(
-        testdata1["data"],
+    masked_data = testdata1["data"][testdata1["mask"], ...]
+    masked_adaptive_mask = testdata1["adaptive_mask"][testdata1["mask"]]
+    t2s, s0, failures, t2s_var, s0_var, t2s_s0_covar = me.fit_decay_ts(
+        masked_data,
         testdata1["tes"],
-        testdata1["mask"],
-        testdata1["adaptive_mask"],
+        masked_adaptive_mask,
         testdata1["fittype"],
     )
-    assert t2sv.ndim == 2
-    assert s0v.ndim == 2
-    assert t2svg.ndim == 2
-    assert s0vg.ndim == 2
+    assert t2s.ndim == 2
+    assert s0.ndim == 2
+    assert failures is None
+    assert t2s_var is None
+    assert s0_var is None
+    assert t2s_s0_covar is None
 
 
 def test__apply_t2s_floor():
@@ -110,21 +116,27 @@ def test_smoke_fit_decay():
     mask[n_samples // 2 :] = 0
     adaptive_mask = np.random.randint(2, n_echos, size=(n_samples)) * mask
     fittype = "loglin"
-    t2s_limited, s0_limited, t2s_full, s0_full, _ = me.fit_decay(
-        data, tes, mask, adaptive_mask, fittype
+    masked_data = data[mask, ...]
+    masked_adaptive_mask = adaptive_mask[mask]
+    t2s, s0, failures, t2s_var, s0_var, t2s_s0_covar = me.fit_decay(
+        masked_data,
+        tes,
+        masked_adaptive_mask,
+        fittype,
     )
-    assert t2s_limited is not None
-    assert s0_limited is not None
-    assert t2s_full is not None
-    assert s0_full is not None
+    assert t2s.ndim == 1
+    assert s0.ndim == 1
+    assert failures is None
+    assert t2s_var is None
+    assert s0_var is None
+    assert t2s_s0_covar is None
 
 
 def test_smoke_fit_decay_curvefit():
     """
     Test_smoke_fit_decay tests that the function fit_decay returns reasonable.
 
     objects with random inputs in the correct format when using the direct.
-
     monoexponetial approach.
     """
     n_samples = 100
@@ -136,13 +148,20 @@ def test_smoke_fit_decay_curvefit():
     mask[n_samples // 2 :] = 0
     adaptive_mask = np.random.randint(2, n_echos, size=(n_samples)) * mask
     fittype = "curvefit"
-    t2s_limited, s0_limited, t2s_full, s0_full, _ = me.fit_decay(
-        data, tes, mask, adaptive_mask, fittype
+    masked_data = data[mask, ...]
+    masked_adaptive_mask = adaptive_mask[mask]
+    t2s, s0, failures, t2s_var, s0_var, t2s_s0_covar = me.fit_decay(
+        masked_data,
+        tes,
+        masked_adaptive_mask,
+        fittype,
     )
-    assert t2s_limited is not None
-    assert s0_limited is not None
-    assert t2s_full is not None
-    assert s0_full is not None
+    assert t2s.ndim == 1
+    assert s0.ndim == 1
+    assert failures.ndim == 1
+    assert t2s_var.ndim == 1
+    assert s0_var.ndim == 1
+    assert t2s_s0_covar.ndim == 1
 
 
 def test_smoke_fit_decay_ts():
@@ -160,21 +179,27 @@ def test_smoke_fit_decay_ts():
     mask[n_samples // 2 :] = 0
     adaptive_mask = np.random.randint(2, n_echos, size=(n_samples)) * mask
     fittype = "loglin"
-    t2s_limited_ts, s0_limited_ts, t2s_full_ts, s0_full_ts, _ = me.fit_decay_ts(
-        data, tes, mask, adaptive_mask, fittype
+    masked_data = data[mask, ...]
+    masked_adaptive_mask = adaptive_mask[mask]
+    t2s, s0, failures, t2s_var, s0_var, t2s_s0_covar = me.fit_decay_ts(
+        masked_data,
+        tes,
+        masked_adaptive_mask,
+        fittype,
     )
-    assert t2s_limited_ts is not None
-    assert s0_limited_ts is not None
-    assert t2s_full_ts is not None
-    assert s0_full_ts is not None
+    assert t2s.ndim == 2
+    assert s0.ndim == 2
+    assert failures is None
+    assert t2s_var is None
+    assert s0_var is None
+    assert t2s_s0_covar is None
 
 
 def test_smoke_fit_decay_curvefit_ts():
     """
     Test_smoke_fit_decay_ts tests that the function fit_decay_ts returns reasonable.
 
     objects with random inputs in the correct format when using the direct.
-
     monoexponetial approach.
     """
     n_samples = 100
@@ -186,13 +211,20 @@ def test_smoke_fit_decay_curvefit_ts():
     mask[n_samples // 2 :] = 0
     adaptive_mask = np.random.randint(2, n_echos, size=(n_samples)) * mask
     fittype = "curvefit"
-    t2s_limited_ts, s0_limited_ts, t2s_full_ts, s0_full_ts, _ = me.fit_decay_ts(
-        data, tes, mask, adaptive_mask, fittype
+    masked_data = data[mask, ...]
+    masked_adaptive_mask = adaptive_mask[mask]
+    t2s, s0, failures, t2s_var, s0_var, t2s_s0_covar = me.fit_decay_ts(
+        masked_data,
+        tes,
+        masked_adaptive_mask,
+        fittype,
     )
-    assert t2s_limited_ts is not None
-    assert s0_limited_ts is not None
-    assert t2s_full_ts is not None
-    assert s0_full_ts is not None
+    assert t2s.ndim == 2
+    assert s0.ndim == 2
+    assert failures.ndim == 2
+    assert t2s_var.ndim == 2
+    assert s0_var.ndim == 2
+    assert t2s_s0_covar.ndim == 2
 
 
 # TODO: BREAK AND UNIT TESTS