Reuse saved ICA component maps when generating component figures

tedana/io.py

@@ -769,7 +769,7 @@ def write_split_ts(data, mixing, mask, component_table, io_generator, echo=0):
     LGR.info(f"Writing denoised time series: {fout}")
 
 
-def writeresults(data_optcom, mask, component_table, mixing, io_generator):
+def writeresults(data_optcom, mask, component_table, mixing, io_generator, mixing_orig=None):
     """Denoise `ts` and save all resulting files to disk.
 
     Parameters
@@ -782,11 +782,14 @@ def writeresults(data_optcom, mask, component_table, mixing, io_generator):
         Component metric table. One row for each component, with a column for
         each metric. Requires at least two columns: "component" and
         "classification".
-    mixing : (C x T) array_like
+    mixing : (T x C) array_like
         Mixing matrix for converting input data to component space, where `C`
         is components and `T` is the same as in `data`
-    ref_img : :obj:`str` or img_like
+    io_generator : :obj:`tedana.io.OutputGenerator`
         Reference image to dictate how outputs are saved to disk
+    mixing_orig : (T x C) array_like or None, optional
+        Non-orthogonalized mixing matrix.
+        If defined, it will be used for the component maps, but not for denoising.
 
     See Also
     --------
@@ -807,17 +810,30 @@ def writeresults(data_optcom, mask, component_table, mixing, io_generator):
     desc-ICAAccepted_components.nii.gz           Spatial component maps for accepted components.
     desc-ICAAccepted_stat-z_components.nii.gz    Z-normalized spatial component maps
                                                  for accepted components.
+    desc-ICAOrth_components.nii.gz               Spatial component maps for all components from
+                                                 orthogonalized mixing matrix. Only created if
+                                                 tedort and verbose are enabled.
     =========================================    ===============================================
     """
     acc = component_table[component_table.classification == "accepted"].index.values
     write_split_ts(data_optcom, mixing, mask, component_table, io_generator)
 
-    ts_pes = get_coeffs(data_optcom, mixing)
+    if mixing_orig is not None:
+        mixing_for_denoising = mixing
+        mixing_for_components = mixing_orig
+    else:
+        mixing_for_denoising = mixing_for_components = mixing
+
+    ts_pes = get_coeffs(data_optcom, mixing_for_components)
     fout = io_generator.save_file(ts_pes, "ICA components img")
     LGR.info(f"Writing full ICA coefficient feature set: {fout}")
+    if io_generator.verbose and (mixing_orig is not None):
+        ts_pes_orth = get_coeffs(data_optcom, mixing_for_denoising)
+        fout = io_generator.save_file(ts_pes_orth, "ICA orthogonalized components img")
+        LGR.info(f"Writing orthogonalized ICA coefficient feature set: {fout}")
 
     data_optcom_z = stats.zscore(data_optcom[mask, :], axis=-1)
-    mixing_z = stats.zscore(mixing, axis=0)
+    mixing_z = stats.zscore(mixing_for_components, axis=0)
     betas_oc = get_coeffs(data_optcom_z, mixing_z)
     fout = io_generator.save_file(betas_oc, "z-scored ICA components img", mask=mask)
     del data_optcom_z, mixing_z

tedana/reporting/static_figures.py

@@ -14,7 +14,7 @@
 import matplotlib.pyplot as plt
 from nilearn import image, masking, plotting
 
-from tedana import io, stats, utils
+from tedana import io, utils
 
 LGR = logging.getLogger("GENERAL")
 MPL_LGR = logging.getLogger("matplotlib")
@@ -333,95 +333,121 @@ def plot_component(
     plt.close(fig)
 
 
-def comp_figures(ts, component_table, mixing, io_generator, png_cmap):
+def _generate_single_component_figure(
+    compnum,
+    component_table,
+    mixing,
+    component_img,
+    tr,
+    png_cmap,
+    out_dir,
+    prefix,
+):
+    """Generate a figure for a single component.
+
+    Parameters
+    ----------
+    compnum : int
+        The component number to plot.
+    component_table : (C x M) :obj:`pandas.DataFrame`
+        Component metric table.
+    mixing : (T x C) array_like
+        Mixing matrix for converting input data to component space.
+    component_img : :obj:`nibabel.spatialimages.SpatialImage`
+        Spatial map for the component.
+    tr : float
+        Repetition time of the time series.
+    png_cmap : str
+        Colormap to use for the spatial map.
+    out_dir : str
+        Output directory path.
+    prefix : str
+        Prefix for the output file name.
+    """
+    classification = component_table.loc[compnum, "classification"]
+    classification_tags = str(component_table.loc[compnum, "classification_tags"])
+
+    color_map = {"accepted": "g", "rejected": "r", "ignored": "k"}
+    line_color = color_map.get(classification, "0.75")
+
+    if classification in color_map:
+        expl_text = f"{classification} reason(s): {classification_tags}"
+    else:
+        expl_text = "other classification"
+
+    # Title will include variance from component_table
+    comp_var = f"{component_table.loc[compnum, 'variance explained']:.2f}"
+    comp_kappa = f"{component_table.loc[compnum, 'kappa']:.2f}"
+    comp_rho = f"{component_table.loc[compnum, 'rho']:.2f}"
+
+    plt_title = (
+        f"Comp. {compnum}: variance: {comp_var}%, kappa: {comp_kappa}, "
+        f"rho: {comp_rho}, {expl_text}"
+    )
+
+    component_timeseries = mixing[:, compnum]
+
+    # Get fft and freqs for this component
+    # adapted from @dangom
+    spectrum, freqs = utils.get_spectrum(component_timeseries, tr)
+
+    plot_name = f"{prefix}comp_{str(compnum).zfill(3)}.png"
+    compplot_name = os.path.join(out_dir, "figures", plot_name)
+
+    plot_component(
+        stat_img=component_img,
+        component_timeseries=component_timeseries,
+        power_spectrum=spectrum,
+        frequencies=freqs,
+        tr=tr,
+        classification_color=line_color,
+        png_cmap=png_cmap,
+        title=plt_title,
+        out_file=compplot_name,
+    )
+
+
+def comp_figures(component_table, mixing, io_generator, png_cmap):
     """Create static figures that highlight certain aspects of tedana processing.
 
     This includes a figure for each component showing the component time course,
     the spatial weight map and a fast Fourier transform of the time course.
 
     Parameters
     ----------
-    ts : (Mb x T) array_like
-        Time series from which to derive ICA betas, where `Mb` is samples in base mask,
-        and `T` is time
     component_table : (C x X) :obj:`pandas.DataFrame`
         Component metric table. One row for each component, with a column for
         each metric. The index should be the component number.
-    mixing : (C x T) array_like
+    mixing : (T x C) array_like
         Mixing matrix for converting input data to component space, where `C`
         is components and `T` is the same as in `data`
     io_generator : :obj:`tedana.io.OutputGenerator`
         Output Generator object to use for this workflow
+    png_cmap : str
+        Colormap to use for the spatial map.
     """
-    # regenerate the beta images
-    component_maps_arr = stats.get_coeffs(ts, mixing)
-    component_maps_arr = masking.unmask(component_maps_arr.T, io_generator.mask)
-    component_maps_arr = component_maps_arr.get_fdata()
+    component_betas_file = io_generator.get_name("ICA components img")
+    component_betas_img = nb.load(component_betas_file)
+    component_betas_arr = np.asanyarray(component_betas_img.dataobj)
 
     # Get repetition time from reference image
     tr = io_generator.reference_img.header.get_zooms()[-1]
 
-    # Remove trailing ';' from rationale column
-    # component_table["rationale"] = component_table["rationale"].str.rstrip(";")
     for compnum in component_table.index.values:
-        if component_table.loc[compnum, "classification"] == "accepted":
-            line_color = "g"
-            expl_text = "accepted reason(s): " + str(
-                component_table.loc[compnum, "classification_tags"]
-            )
-
-        elif component_table.loc[compnum, "classification"] == "rejected":
-            line_color = "r"
-            expl_text = "rejected reason(s): " + str(
-                component_table.loc[compnum, "classification_tags"]
-            )
-
-        elif component_table.loc[compnum, "classification"] == "ignored":
-            line_color = "k"
-            expl_text = "ignored reason(s): " + str(
-                component_table.loc[compnum, "classification_tags"]
-            )
-
-        else:
-            # Classification not added
-            # If new, this will keep code running
-            line_color = "0.75"
-            expl_text = "other classification"
-
-        # Title will include variance from component_table
-        comp_var = f"{component_table.loc[compnum, 'variance explained']:.2f}"
-        comp_kappa = f"{component_table.loc[compnum, 'kappa']:.2f}"
-        comp_rho = f"{component_table.loc[compnum, 'rho']:.2f}"
-
-        plt_title = (
-            f"Comp. {compnum}: variance: {comp_var}%, kappa: {comp_kappa}, "
-            f"rho: {comp_rho}, {expl_text}"
-        )
         component_img = nb.Nifti1Image(
-            component_maps_arr[:, :, :, compnum],
-            affine=io_generator.reference_img.affine,
-            header=io_generator.reference_img.header,
+            component_betas_arr[..., compnum],
+            affine=component_betas_img.affine,
+            header=component_betas_img.header,
         )
-
-        component_timeseries = mixing[:, compnum]
-
-        # Get fft and freqs for this component
-        # adapted from @dangom
-        spectrum, freqs = utils.get_spectrum(component_timeseries, tr)
-
-        plot_name = f"{io_generator.prefix}comp_{str(compnum).zfill(3)}.png"
-        compplot_name = os.path.join(io_generator.out_dir, "figures", plot_name)
-
-        plot_component(
-            stat_img=component_img,
-            component_timeseries=component_timeseries,
-            power_spectrum=spectrum,
-            frequencies=freqs,
+        _generate_single_component_figure(
+            compnum=compnum,
+            component_table=component_table,
+            mixing=mixing,
+            component_img=component_img,
             tr=tr,
-            classification_color=line_color,
             png_cmap=png_cmap,
-            title=plt_title,
-            out_file=compplot_name,
+            out_dir=io_generator.out_dir,
+            prefix=io_generator.prefix,
         )
 
 

tedana/resources/config/outputs.json

@@ -43,6 +43,10 @@
         "orig": "ica_components",
         "bidsv1.5.0": "desc-ICA_components"
     },
+    "ICA orthogonalized components img": {
+        "orig": "ica_orth_components",
+        "bidsv1.5.0": "desc-ICAOrth_components"
+    },
     "z-scored PCA components img": {
         "orig": "pca_components",
         "bidsv1.5.0": "desc-PCA_stat-z_components"

tedana/tests/data/nih_five_echo_outputs_verbose.txt

@@ -14,6 +14,7 @@ sub-01_desc-tedana_registry.json
 sub-01_desc-ICACrossComponent_metrics.json
 sub-01_desc-ICA_status_table.tsv
 sub-01_desc-ICA_decision_tree.json
+sub-01_desc-ICAOrth_components.nii.gz
 sub-01_desc-ICAS0_stat-F_statmap.nii.gz
 sub-01_desc-ICAT2_stat-F_statmap.nii.gz
 sub-01_desc-ICA_mixing.tsv

tedana/tests/test_io.py

@@ -118,6 +118,49 @@ def test_smoke_write_split_ts():
     os.remove(fname)
 
 
+def test_writeresults_uses_original_mixing_for_component_maps():
+    """Component maps should use the original mixing matrix when one is provided."""
+    rng = np.random.default_rng(42)
+    n_samples, n_times, n_components = 5, 8, 3
+    data = rng.normal(size=(n_samples, n_times))
+    mask = np.ones(n_samples, dtype=bool)
+    mixing = rng.normal(size=(n_times, n_components))
+    mixing_orig = rng.normal(size=(n_times, n_components))
+    component_table = pd.DataFrame(
+        {"classification": ["rejected", "ignored", "rejected"]},
+    )
+    saved = {}
+
+    class StubOutputGenerator:
+        verbose = True
+
+        def save_file(self, data, description, **kwargs):
+            del kwargs
+            saved[description] = np.asarray(data).copy()
+            return description
+
+    with mock.patch.object(me, "write_split_ts") as write_split_ts:
+        me.writeresults(
+            data_optcom=data,
+            mask=mask,
+            component_table=component_table,
+            mixing=mixing,
+            io_generator=StubOutputGenerator(),
+            mixing_orig=mixing_orig,
+        )
+
+    write_split_ts.assert_called_once()
+    np.testing.assert_allclose(
+        saved["ICA components img"],
+        me.get_coeffs(data, mixing_orig),
+    )
+    np.testing.assert_allclose(
+        saved["ICA orthogonalized components img"],
+        me.get_coeffs(data, mixing),
+    )
+    assert "z-scored ICA components img" in saved
+
+
 def test_load_data_nilearn_multi_echo_fastpath(tmp_path):
     """`load_data_nilearn` should return (Mb, E, T) for multi-echo files."""
     affine = np.eye(4)