Add QC summary to HTML report

CLAUDE.md

@@ -0,0 +1,18 @@
+# tedana
+
+## Environment
+
+Development uses the micromamba environment **`tedenv`** (python 3.12.13; numpy,
+scipy, nibabel, nilearn, pandas, scikit-learn, mapca, matplotlib). `tedana` is
+installed editable from this checkout, so source edits take effect immediately
+with no reinstall.
+
+Run everything through it:
+
+    micromamba run -n tedenv pytest tedana/tests/test_decay.py -q
+    micromamba run -n tedenv flake8 tedana/decay.py
+    micromamba run -n tedenv black tedana/decay.py
+
+Use `tedenv`, not `tedanapy`. Both environments have tedana and past sessions
+used them interchangeably, but only `tedenv` has the editable install and the
+lint toolchain (`flake8`, `black`, `isort`). This project does not use ruff.

docs/outputs.rst

@@ -58,6 +58,12 @@ tedana_report.html                                                           The
 "s0 img": S0map.nii.gz                                                       Full S0 3D map. If a voxel has at least 1 good
                                                                              echo then the first two echoes will be used to estimate
                                                                              a value
+"decay metrics json": desc-decay_metrics.json                                Precomputed decay-fit QC summary read by the report's
+                                                                             QC summary card: mean/median T2*, S0, and fit RMSE;
+                                                                             base-mask and fit-mask voxel counts; a good-echo voxel
+                                                                             histogram; and, when curve-fitting is used, first-pass
+                                                                             and post-interpolation fit-failure counts. Only written
+                                                                             when fitmode is "all".
 "PCA mixing tsv": desc-PCA_mixing.tsv                                        Mixing matrix (component time series) from PCA
                                                                              decomposition in a tab-delimited file. Each column is
                                                                              a different component, and the column name is the
@@ -99,6 +105,10 @@ tedana_report.html                                                           The
 "ICA cross component metrics json": desc-ICACrossComponent_metrics.json      Metric names and values that are each a single number
                                                                              calculated across components. For example, kappa and
                                                                              rho elbows.
+                                                                             Also includes aggregate variance measures used by the
+                                                                             QC summary card: accepted_variance, rejected_variance,
+                                                                             ignored_variance, unmodeled_variance (relative to the
+                                                                             raw optimally-combined data), and retained_variance.
 "ICA decision tree json": desc-ICA_decision_tree                             A copy of the inputted decision tree specification with
                                                                              an added "output" field for each node. The output field
                                                                              contains information about what happened during

tedana/decay.py

@@ -815,6 +815,99 @@ def rmse_of_fit_decay_ts(
     return rmse_map, rmse_df
 
 
+def generate_decay_metrics(
+    *,
+    t2star,
+    s0,
+    rmse_map,
+    adaptive_mask,
+    n_fit_failures=None,
+    n_fit_failures_after_interpolation=None,
+):
+    """Summarize decay-fit quality into a dict for the QC report.
+
+    All array inputs are in base-mask sample space (one entry per base-mask voxel)
+    and must be 1D (``(S,)``). This function is only meaningful for scalar,
+    per-voxel summaries; time-varying (``fitmode == "ts"``) maps are not supported.
+
+    Parameters
+    ----------
+    t2star, s0, rmse_map : (S,) array_like
+        Full T2*, S0, and fit-RMSE maps in base-mask sample space. Must be 1D.
+    adaptive_mask : (S,) array_like
+        Integer count of good echoes per base-mask voxel (0 means no good echo).
+        Must be 1D.
+    n_fit_failures, n_fit_failures_after_interpolation : int or None
+        Curve-fit failure counts. ``None`` when curve-fitting was not used, in which
+        case the corresponding keys are omitted from the returned dict.
+
+    Returns
+    -------
+    dict
+
+    Raises
+    ------
+    ValueError
+        If any of ``adaptive_mask``, ``t2star``, ``s0``, or ``rmse_map`` is not 1D.
+        This most often happens when ``fitmode == "ts"`` produces 2D
+        (voxels x time) maps, for which a scalar summary is ill-defined.
+    """
+    adaptive_mask = np.asarray(adaptive_mask)
+    for name, arr in (
+        ("adaptive_mask", adaptive_mask),
+        ("t2star", t2star),
+        ("s0", s0),
+        ("rmse_map", rmse_map),
+    ):
+        arr = np.asarray(arr)
+        if arr.ndim != 1:
+            raise ValueError(
+                "generate_decay_metrics expects 1D base-mask-space arrays; "
+                f"got shape {arr.shape} for '{name}'."
+            )
+
+    fit_mask = adaptive_mask >= 1
+
+    def _finite_mean_median(arr):
+        arr = np.asarray(arr, dtype=float)[fit_mask]
+        # Exclude zeros: modify_t2s_s0_maps fills NaN S0 values with 0.0, and
+        # those placeholders would otherwise bias the mean/median.
+        arr = arr[np.isfinite(arr) & (arr != 0)]
+        if arr.size == 0:
+            return None, None
+        return float(np.mean(arr)), float(np.median(arr))
+
+    t2star_mean, t2star_median = _finite_mean_median(t2star)
+    s0_mean, s0_median = _finite_mean_median(s0)
+
+    rmse = np.asarray(rmse_map, dtype=float)[fit_mask]
+    rmse = rmse[np.isfinite(rmse)]
+    rmse_mean = float(np.mean(rmse)) if rmse.size else None
+    rmse_median = float(np.median(rmse)) if rmse.size else None
+
+    good_echo_voxel_counts = {
+        int(n): int((adaptive_mask == n).sum()) for n in np.unique(adaptive_mask[fit_mask])
+    }
+
+    metrics = {
+        "t2star_mean": t2star_mean,
+        "t2star_median": t2star_median,
+        "s0_mean": s0_mean,
+        "s0_median": s0_median,
+        "rmse_mean": rmse_mean,
+        "rmse_median": rmse_median,
+        "n_voxels_base_mask": int(adaptive_mask.size),
+        "n_voxels_fit_mask": int(fit_mask.sum()),
+        "good_echo_voxel_counts": good_echo_voxel_counts,
+    }
+    if n_fit_failures is not None:
+        metrics["n_fit_failures"] = int(n_fit_failures)
+    if n_fit_failures_after_interpolation is not None:
+        metrics["n_fit_failures_after_interpolation"] = int(n_fit_failures_after_interpolation)
+
+    return metrics
+
+
 def t2smap_subworkflow(
     data_cat,
     tes,
@@ -861,6 +954,9 @@ def t2smap_subworkflow(
     t2s_full : (Mb x T) :obj:`numpy.ndarray`
         The full T2* map.
     """
+    n_fit_failures = None
+    n_fit_failures_after_interpolation = None
+
     data_for_fit = (
         data_without_excluded_vols if data_without_excluded_vols is not None else data_cat
     )
@@ -992,6 +1088,11 @@ def t2smap_subworkflow(
         failures_map = first_pass_failures.astype(np.uint8)
         if interpolate_failing_voxels and first_pass_failures.any():
             failures_map += failures.astype(np.uint8)
+
+        n_fit_failures = int(first_pass_failures.sum())
+        if interpolate_failing_voxels and first_pass_failures.any():
+            n_fit_failures_after_interpolation = int(failures.sum())
+
         io_generator.save_file(failures_map, "fit failures img", mask=mask_denoise)
 
         if io_generator.verbose:
@@ -1040,7 +1141,21 @@ def t2smap_subworkflow(
         s0=s0_full,
         fitmode=fitmode,
     )
-    del s0_full
     io_generator.save_file(rmse_map, "rmse img")
     io_generator.save_file(rmse_df, "confounds tsv")
+
+    if fitmode == "all":
+        # In fitmode == "ts", the T2*/S0/RMSE maps are time-varying (Mb x T), so a
+        # single scalar summary per voxel is ill-defined and decay metrics are
+        # simply not written.
+        decay_metrics = generate_decay_metrics(
+            t2star=t2s_full,
+            s0=s0_full,
+            rmse_map=rmse_map,
+            adaptive_mask=masksum_denoise,
+            n_fit_failures=n_fit_failures,
+            n_fit_failures_after_interpolation=n_fit_failures_after_interpolation,
+        )
+        io_generator.save_file(decay_metrics, "decay metrics json")
+
     return t2s_full

tedana/reporting/data/html/report_body_template.html

@@ -424,6 +424,17 @@ <h2 id="ricaModalTitle">Open in Rica</h2>
 
 <!-- Info tab -->
 <div class="tab-pane" id="pane-info" role="tabpanel" aria-labelledby="tab-info" tabindex="0">
+  <div class="qc-summary-card">
+    <h1>tedana QC summary</h1>
+    <table>
+      {% for row in qcCard %}
+      <tr>
+        <th style="text-align: left; padding-right: 1em; white-space: nowrap;">{{ row.label }}</th>
+        <td>{{ row.value }}</td>
+      </tr>
+      {% endfor %}
+    </table>
+  </div>
   <div class="info">
     <h1>Info</h1>
     {{ info }}

tedana/reporting/html_report.py

@@ -129,6 +129,7 @@ def _update_template_bokeh(
     tsne,
     tree_table,
     status_table,
+    qc_card,
 ):
     """
     Populate a report with content.
@@ -155,6 +156,8 @@ def _update_template_bokeh(
         HTML table of decision tree nodes created by _generate_tree_tables()
     status_table : str or None
         HTML table of component statuses created by _generate_tree_tables()
+    qc_card : list of dict
+        Display rows created by _generate_qc_card()
 
     Returns
     -------
@@ -289,6 +292,7 @@ def _update_template_bokeh(
         treeExists=tree_exists,
         treeTable=tree_table,
         statusTable=status_table,
+        qcCard=qc_card,
     )
     return body
 
@@ -390,6 +394,98 @@ def _generate_tree_tables(io_generator):
     return tree_table, status_table
 
 
+def _fmt_num(value, suffix="", decimals=1):
+    """Format a possibly-None number for display, returning 'n/a' for None."""
+    if value is None:
+        return "n/a"
+    return f"{value:.{decimals}f}{suffix}"
+
+
+def _generate_qc_card(
+    *,
+    component_table,
+    cross_comp_metrics_dict,
+    decay_metrics_dict,
+    kappa_elbow,
+    rho_elbow,
+    n_vols,
+    n_comps,
+    tree_node_count,
+    version,
+):
+    """Assemble display rows of precomputed run-level QC values for the summary card.
+
+    Returns a list of ``{"label": str, "value": str}`` rows. Performs no scientific
+    computation; every value is read from precomputed inputs.
+    """
+    counts = component_table["classification"].value_counts().to_dict()
+    n_total = int(len(component_table))
+    n_accepted = int(counts.get("accepted", 0))
+    n_rejected = int(counts.get("rejected", 0))
+    n_ignored = int(counts.get("ignored", 0))
+
+    ccm = cross_comp_metrics_dict or {}
+
+    rows = [
+        {
+            "label": "Components",
+            "value": (
+                f"{n_total} total | {n_accepted} accepted | "
+                f"{n_rejected} rejected | {n_ignored} ignored"
+            ),
+        },
+        {"label": "Variance accepted", "value": _fmt_num(ccm.get("accepted_variance"), "%")},
+        {"label": "Variance rejected", "value": _fmt_num(ccm.get("rejected_variance"), "%")},
+        {"label": "Variance unmodeled", "value": _fmt_num(ccm.get("unmodeled_variance"), "%")},
+        {"label": "Variance retained", "value": _fmt_num(ccm.get("retained_variance"), "%")},
+        {"label": "Kappa elbow", "value": _fmt_num(kappa_elbow, decimals=2)},
+        {"label": "Rho elbow", "value": _fmt_num(rho_elbow, decimals=2)},
+        {
+            "label": "Dimensions",
+            "value": (
+                f"{n_vols} volumes | {ccm.get('n_echos', 'n/a')} echoes | " f"{n_comps} components"
+            ),
+        },
+        {
+            "label": "Decision-tree nodes",
+            "value": "n/a" if tree_node_count is None else str(tree_node_count),
+        },
+        {"label": "tedana version", "value": str(version)},
+    ]
+
+    if decay_metrics_dict:
+        rows.append(
+            {"label": "Mean T2*", "value": _fmt_num(decay_metrics_dict.get("t2star_mean"), " ms")}
+        )
+        rows.append(
+            {"label": "Median RMSE", "value": _fmt_num(decay_metrics_dict.get("rmse_median"))}
+        )
+        fit_vox = decay_metrics_dict.get("n_voxels_fit_mask")
+        base_vox = decay_metrics_dict.get("n_voxels_base_mask")
+        fit_vox_txt = "n/a" if fit_vox is None else str(fit_vox)
+        base_vox_txt = "n/a" if base_vox is None else str(base_vox)
+        rows.append(
+            {
+                "label": "Fit-mask voxels",
+                "value": f"{fit_vox_txt} of {base_vox_txt} base-mask voxels",
+            }
+        )
+        if decay_metrics_dict.get("n_fit_failures") is not None:
+            after = decay_metrics_dict.get("n_fit_failures_after_interpolation")
+            after_txt = "n/a" if after is None else str(after)
+            rows.append(
+                {
+                    "label": "Fit failures",
+                    "value": (
+                        f"{decay_metrics_dict['n_fit_failures']} first-pass | "
+                        f"{after_txt} after interpolation"
+                    ),
+                }
+            )
+
+    return rows
+
+
 def generate_report(io_generator: OutputGenerator, cluster_labels, similarity_t_sne) -> None:
     """Generate an HTML report.
 
@@ -545,6 +641,33 @@ def get_elbow_val(elbow_prefix):
     # Create the decision tree tables
     tree_table, status_table = _generate_tree_tables(io_generator)
 
+    # Read the precomputed decay metrics, if present.
+    decay_metrics_dict = None
+    decay_metrics_path = io_generator.get_name("decay metrics json")
+    if os.path.exists(decay_metrics_path):
+        decay_metrics_dict = load_json(decay_metrics_path)
+
+    # Number of decision-tree nodes, if the tree JSON exists.
+    tree_node_count = None
+    tree_path = io_generator.get_name("ICA decision tree json")
+    if os.path.exists(tree_path):
+        tree_data = load_json(tree_path)
+        nodes = tree_data.get("nodes")
+        if nodes is not None:
+            tree_node_count = len(nodes)
+
+    qc_card = _generate_qc_card(
+        component_table=component_table,
+        cross_comp_metrics_dict=cross_comp_metrics_dict,
+        decay_metrics_dict=decay_metrics_dict,
+        kappa_elbow=kappa_elbow,
+        rho_elbow=rho_elbow,
+        n_vols=n_vols,
+        n_comps=n_comps,
+        tree_node_count=tree_node_count,
+        version=__version__,
+    )
+
     body = _update_template_bokeh(
         bokeh_id=kr_div,
         info_table=info_table,
@@ -556,6 +679,7 @@ def get_elbow_val(elbow_prefix):
         tsne=tsne_html,
         tree_table=tree_table,
         status_table=status_table,
+        qc_card=qc_card,
     )
     html = _save_as_html(body)
     with open(opj(io_generator.out_dir, f"{io_generator.prefix}tedana_report.html"), "wb") as f: