Add T2* and S0 figures (#1040)

* Add T2* and S0 figures.

* Add seaborn as a dependency.

* Update static_figures.py

* Fix.

* Update expected outputs.

* Reclassify doesn't need T2*/S0 plots.

* Add png_cmap to T2*/S0 plots.

* Add type hints and fix test outputs.

* Update html_report.py

* Update report_body_template.html

* Try improving the figures.

* Update report_body_template.html

* Drop seaborn.

* Update static_figures.py

* Update static_figures.py

* Update static_figures.py

* Update static_figures.py

* Document the figures.

* Update docs/outputs.rst

Co-authored-by: Dan Handwerker <7406227+handwerkerd@users.noreply.github.com>

* Switch to grayscale colormap.

* Update tedana/reporting/static_figures.py

Co-authored-by: Dan Handwerker <7406227+handwerkerd@users.noreply.github.com>

* Update tedana/reporting/static_figures.py

Co-authored-by: Dan Handwerker <7406227+handwerkerd@users.noreply.github.com>

* Update static_figures.py

* Update figure.

---------

Co-authored-by: Dan Handwerker <7406227+handwerkerd@users.noreply.github.com>

docs/outputs.rst

@@ -527,6 +527,22 @@ should not overly focus on carpet plots and should examine these results in cont
   :height: 400px
 
 
+************************
+T2* and S0 Summary Plots
+************************
+
+Below the carpet plots are summary plots for the T2* and S0 maps.
+Each map has two figures: a spatial map of the values and a histogram of the voxelwise values.
+The T2* map should look similar to T2 maps and be brightest in the ventricles and darkest in areas of largest susceptibility.
+The S0 map should roughly follow the signal-to-noise ratio and will be brightest near the surface near RF coils.
+
+It is important to note that the histogram is limited from 0 to the 98th percentile of the data to improve readability.
+
+.. image:: /_static/t2star_plots.png
+  :align: center
+  :height: 400px
+
+
 **************************
 Citable workflow summaries
 **************************

tedana/reporting/data/html/report_body_template.html

@@ -1,4 +1,9 @@
 <style type="text/css">
+  .brainplot {
+    width: 50%;
+    height: auto;
+  }
+
   .report {
     display: flex;
     flex-direction: column;
@@ -174,21 +179,39 @@ <h1>Carpet plots</h1>
       <img id="imgCarpetPlot" src="$initialCarpet" />
     </div>
   </div>
-  <div class="info">
-    <h1>Info</h1>
-    $info
+</div>
+<div class="carpet-plots">
+  <h1>T2* and S0</h1>
+  <h2>T2*</h2>
+  <div class="carpet-plots-image">
+    <img id="t2starBrainPlot" src="$t2starBrainPlot" class="brainplot" />
+  </div>
+  <div class="carpet-plots-image">
+    <img id="t2starHistogram" src="$t2starHistogram" />
   </div>
-  <div class="about">
-    <h1>About tedana</h1>
-    $about
+  <h2>S0</h2>
+  <div class="carpet-plots-image">
+    <img id="s0BrainPlot" src="$s0BrainPlot" class="brainplot" />
   </div>
-  <div class="content references">
-    <h1>References</h1>
-    <ul>
-      $references
-    </ul>
+  <div class="carpet-plots-image">
+    <img id="s0Histogram" src="$s0Histogram" />
   </div>
 </div>
+</div>
+<div class="info">
+  <h1>Info</h1>
+  $info
+</div>
+<div class="about">
+  <h1>About tedana</h1>
+  $about
+</div>
+<div class="content references">
+  <h1>References</h1>
+  <ul>
+    $references
+  </ul>
+</div>
 
 <script type="text/javascript">
   function updateCarpetPlot(name) {

tedana/reporting/html_report.py

@@ -14,7 +14,7 @@
 from pybtex.plugin import find_plugin
 
 from tedana import __version__
-from tedana.io import load_json
+from tedana.io import OutputGenerator, load_json
 from tedana.reporting import dynamic_figures as df
 
 LGR = logging.getLogger("GENERAL")
@@ -143,6 +143,12 @@ def _update_template_bokeh(bokeh_id, info_table, about, prefix, references, boke
     # Initial carpet plot (default one)
     initial_carpet = f"./figures/{prefix}carpet_optcom.svg"
 
+    # T2* and S0 images
+    t2star_brain = f"./figures/{prefix}t2star_brain.svg"
+    t2star_histogram = f"./figures/{prefix}t2star_histogram.svg"
+    s0_brain = f"./figures/{prefix}s0_brain.svg"
+    s0_histogram = f"./figures/{prefix}s0_histogram.svg"
+
     # Convert bibtex to html
     references, bibliography = _bib2html(references)
 
@@ -159,6 +165,10 @@ def _update_template_bokeh(bokeh_id, info_table, about, prefix, references, boke
         about=about,
         prefix=prefix,
         initialCarpet=initial_carpet,
+        t2starBrainPlot=t2star_brain,
+        t2starHistogram=t2star_histogram,
+        s0BrainPlot=s0_brain,
+        s0Histogram=s0_histogram,
         references=references,
         javascript=bokeh_js,
         buttons=buttons,
@@ -212,18 +222,17 @@ def _generate_info_table(info_dict):
     return info_html
 
 
-def generate_report(io_generator):
+def generate_report(io_generator: OutputGenerator) -> None:
     """Generate an HTML report.
 
     Parameters
     ----------
-    io_generator : tedana.io.OutputGenerator
+    io_generator : :obj:`tedana.io.OutputGenerator`
         io_generator object for this workflow's output
 
-    Returns
-    -------
-    HTML : file
-        A generated HTML report
+    Notes
+    -----
+    This writes out an HTML report to a file.
     """
     # Load the component time series
     comp_ts_path = io_generator.get_name("ICA mixing tsv")

tedana/reporting/static_figures.py

@@ -8,7 +8,7 @@
 
 matplotlib.use("AGG")
 import matplotlib.pyplot as plt
-from nilearn import plotting
+from nilearn import masking, plotting
 
 from tedana import io, stats, utils
 
@@ -463,3 +463,80 @@ def pca_results(criteria, n_components, all_varex, io_generator):
     pca_variance_explained_name = os.path.join(io_generator.out_dir, "figures", plot_name)
     plt.savefig(pca_variance_explained_name)
     plt.close()
+
+
+def plot_t2star_and_s0(
+    *,
+    io_generator: io.OutputGenerator,
+    mask: np.ndarray,
+) -> None:
+    """Create T2* and S0 maps and histograms.
+
+    Parameters
+    ----------
+    io_generator : :obj:`~tedana.io.OutputGenerator`
+        The output generator for this workflow
+    mask : (S,) :obj:`numpy.ndarray`
+        Binary mask used to apply to the data.
+    """
+    t2star_img = io_generator.get_name("t2star img")
+    s0_img = io_generator.get_name("s0 img")
+    mask_img = io.new_nii_like(io_generator.reference_img, mask.astype(int))
+    assert os.path.isfile(t2star_img), f"File {t2star_img} does not exist"
+    assert os.path.isfile(s0_img), f"File {s0_img} does not exist"
+
+    # Plot histograms
+    t2star_data = masking.apply_mask(t2star_img, mask_img)
+    t2s_p02, t2s_p98 = np.percentile(t2star_data, [2, 98])
+    t2star_histogram = f"{io_generator.prefix}t2star_histogram.svg"
+
+    fig, ax = plt.subplots(figsize=(10, 6))
+    ax.hist(t2star_data[t2star_data <= t2s_p98], bins=100)
+    ax.set_xlim(0, t2s_p98)
+    ax.set_title("T2*", fontsize=20)
+    ax.set_ylabel("Count", fontsize=16)
+    ax.set_xlabel("Seconds\n(limited to 98th percentile)", fontsize=16)
+    fig.tight_layout()
+    fig.savefig(os.path.join(io_generator.out_dir, "figures", t2star_histogram))
+
+    s0_data = masking.apply_mask(s0_img, mask_img)
+    s0_p02, s0_p98 = np.percentile(s0_data, [2, 98])
+    s0_histogram = f"{io_generator.prefix}s0_histogram.svg"
+
+    fig, ax = plt.subplots(figsize=(10, 6))
+    ax.hist(s0_data[s0_data <= s0_p98], bins=100)
+    ax.set_xlim(0, s0_p98)
+    ax.set_title("S0", fontsize=20)
+    ax.set_ylabel("Count", fontsize=16)
+    ax.set_xlabel("Arbitrary Units\n(limited to 98th percentile)", fontsize=16)
+    fig.tight_layout()
+    fig.savefig(os.path.join(io_generator.out_dir, "figures", s0_histogram))
+
+    # Plot T2* and S0 maps
+    t2star_plot = f"{io_generator.prefix}t2star_brain.svg"
+    plotting.plot_stat_map(
+        t2star_img,
+        bg_img=None,
+        display_mode="mosaic",
+        symmetric_cbar=False,
+        black_bg=True,
+        cmap="gray",
+        vmin=t2s_p02,
+        vmax=t2s_p98,
+        annotate=False,
+        output_file=os.path.join(io_generator.out_dir, "figures", t2star_plot),
+    )
+
+    s0_plot = f"{io_generator.prefix}s0_brain.svg"
+    plotting.plot_stat_map(
+        s0_img,
+        bg_img=None,
+        display_mode="mosaic",
+        symmetric_cbar=False,
+        black_bg=True,
+        cmap="gray",
+        vmin=s0_p02,
+        vmax=s0_p98,
+        annotate=False,
+        output_file=os.path.join(io_generator.out_dir, "figures", s0_plot),
+    )

tedana/tests/data/cornell_three_echo_outputs.txt

@@ -98,6 +98,10 @@ figures/comp_067.png
 figures/comp_068.png
 figures/pca_criteria.png
 figures/pca_variance_explained.png
+figures/s0_brain.svg
+figures/s0_histogram.svg
+figures/t2star_brain.svg
+figures/t2star_histogram.svg
 references.bib
 report.txt
 tedana_report.html

tedana/tests/data/fiu_four_echo_outputs.txt

@@ -89,3 +89,7 @@ figures/sub-01_comp_017.png
 figures/sub-01_comp_018.png
 figures/sub-01_comp_019.png
 figures/sub-01_comp_020.png
+figures/sub-01_s0_brain.svg
+figures/sub-01_s0_histogram.svg
+figures/sub-01_t2star_brain.svg
+figures/sub-01_t2star_histogram.svg

tedana/tests/data/nih_five_echo_outputs_verbose.txt

@@ -138,3 +138,7 @@ figures/sub-01_comp_048.png
 figures/sub-01_comp_049.png
 figures/sub-01_comp_050.png
 figures/sub-01_comp_051.png
+figures/sub-01_s0_brain.svg
+figures/sub-01_s0_histogram.svg
+figures/sub-01_t2star_brain.svg
+figures/sub-01_t2star_histogram.svg

tedana/workflows/tedana.py

@@ -902,6 +902,7 @@ def tedana_workflow(
             io_generator=io_generator,
             png_cmap=png_cmap,
         )
+        reporting.static_figures.plot_t2star_and_s0(io_generator=io_generator, mask=mask_denoise)
 
         LGR.info("Generating dynamic report")
         reporting.generate_report(io_generator)