[ENH] Update Figure Generation Code

tedana/utils.py

@@ -231,3 +231,23 @@ def andb(arrs):
     result = np.sum(arrs, axis=0)
 
     return result
+
+
+def get_spectrum(data: np.array, tr: float = 1):
+    """
+    Returns the power spectrum and corresponding frequencies when provided
+    with a component time course and repitition time.
+
+    Parameters
+    ----------
+    data : (S, ) array_like
+            A timeseries S, on which you would like to perform an fft.
+    tr : :obj:`float`
+            Reptition time (TR) of the data
+    """
+
+    # adapted from @dangom
+    power_spectrum = np.abs(np.fft.rfft(data)) ** 2
+    freqs = np.fft.rfftfreq(power_spectrum.size * 2 - 1, tr)
+    idx = np.argsort(freqs)
+    return power_spectrum[idx], freqs[idx]

tedana/viz.py

@@ -5,15 +5,43 @@
 import os
 
 import numpy as np
+import matplotlib
+matplotlib.use('AGG')
 import matplotlib.pyplot as plt
 
 from tedana import model
+from tedana.utils import get_spectrum
 
 LGR = logging.getLogger(__name__)
+MPL_LGR = logging.getLogger('matplotlib')
+MPL_LGR.setLevel(logging.WARNING)
 
 
-def write_comp_figs(ts, mask, comptable, mmix, n_vols,
-                    acc, rej, midk, empty, ref_img):
+def trim_edge_zeros(arr):
+    mask = arr != 0
+    bounding_box = tuple(
+                         slice(np.min(indexes), np.max(indexes) + 1)
+                         for indexes in np.where(mask))
+    return arr[bounding_box]
+    """
+    Trims away the zero-filled slices that surround many 3/4D arrays
+
+    Parameters
+    ----------
+    ndarray: (S x T) array_like
+        an array with signal, surrounded by slices that contain only zeros
+        that should be removed.
+
+    Returnes
+    ---------
+    ndarray: (S x T) array_like
+        an array with reduced dimensions, such that the array contains only
+        non_zero values from edge to edge.
+    """
+
+
+def write_comp_figs(ts, mask, comptable, mmix, ref_img, out_dir,
+                    png_cmap):
     """
     Creates static figures that highlight certain aspects of tedana processing
     This includes a figure for each component showing the component time course,
@@ -32,55 +60,60 @@ def write_comp_figs(ts, mask, comptable, mmix, n_vols,
     mmix : (C x T) array_like
         Mixing matrix for converting input data to component space, where `C`
         is components and `T` is the same as in `data`
-    n_vols : :obj:`int`
-        Number of volumes in original time series
-    acc : :obj:`list`
-        Indices of accepted (BOLD) components in `mmix`
-    rej : :obj:`list`
-        Indices of rejected (non-BOLD) components in `mmix`
-    midk : :obj:`list`
-        Indices of mid-K (questionable) components in `mmix`
-    empty : :obj:`list`
-        Indices of ignored components in `mmix`
     ref_img : :obj:`str` or img_like
         Reference image to dictate how outputs are saved to disk
 
     """
+    # Get the lenght of the timeseries
+    n_vols = len(mmix)
 
+    # Check that colormap provided exists
+    if png_cmap not in plt.colormaps():
+        LGR.warning('Provided colormap is not valid, proceeding with coolwarm')
+        png_cmap = 'coolwarm'
     # regenerate the beta images
     ts_B = model.get_coeffs(ts, mmix, mask)
     ts_B = ts_B.reshape(ref_img.shape[:3] + ts_B.shape[1:])
-    # Mask out zeros
+    # trim edges from ts_B array
+    ts_B = trim_edge_zeros(ts_B)
+
+    # Mask out remaining zeros
     ts_B = np.ma.masked_where(ts_B == 0, ts_B)
 
     # Get repetition time from ref_img
     tr = ref_img.header.get_zooms()[-1]
 
-    # Start making plots
-    if not os.path.isdir('figures'):
-        os.mkdir('figures')
-
     # Create indices for 6 cuts, based on dimensions
     cuts = [ts_B.shape[dim] // 6 for dim in range(3)]
+    expl_text = ''
 
+    # Remove trailing ';' from rationale column
+    comptable['rationale'] = comptable['rationale'].str.rstrip(';')
     for compnum in range(0, mmix.shape[1], 1):
 
-        if compnum in acc:
+        if comptable.iloc[compnum]["classification"] == 'accepted':
             line_color = 'g'
-        elif compnum in rej:
+            expl_text = 'accepted'
+        elif comptable.iloc[compnum]["classification"] == 'rejected':
             line_color = 'r'
-        elif compnum in midk:
-            line_color = 'm'
-        else:
+            expl_text = 'rejection reason(s): ' + comptable.iloc[compnum]["rationale"]
+        elif comptable.iloc[compnum]["classification"] == 'ignored':
             line_color = 'k'
+            expl_text = 'ignored reason(s): ' + comptable.iloc[compnum]["rationale"]
+        else:
+            # Classification not added
+            # If new, this will keep code running
+            line_color = '0.75'
+            expl_text = 'unknown classification'
 
         allplot = plt.figure(figsize=(10, 9))
         ax_ts = plt.subplot2grid((5, 6), (0, 0),
                                  rowspan=1, colspan=6,
                                  fig=allplot)
 
         ax_ts.set_xlabel('TRs')
-        ax_ts.set_xbound(0, n_vols)
+        ax_ts.set_xlim(0, n_vols)
+        plt.yticks([])
         # Make a second axis with units of time (s)
         max_xticks = 10
         xloc = plt.MaxNLocator(max_xticks)
@@ -95,18 +128,19 @@ def write_comp_figs(ts, mask, comptable, mmix, n_vols,
             seconds_val = round(X * tr, 2)
             ax2Xs.append(seconds_val)
         ax_ts2.set_xticks(ax1Xs)
-        ax_ts2.set_xbound(ax_ts.get_xbound())
+        ax_ts2.set_xlim(ax_ts.get_xbound())
         ax_ts2.set_xticklabels(ax2Xs)
         ax_ts2.set_xlabel('seconds')
 
         ax_ts.plot(mmix[:, compnum], color=line_color)
 
         # Title will include variance from comptable
-        comp_var = "{0:.2f}".format(comptable.iloc[compnum][3])
-        comp_kappa = "{0:.2f}".format(comptable.iloc[compnum][1])
-        comp_rho = "{0:.2f}".format(comptable.iloc[compnum][2])
-        plt_title = 'Comp. {}: variance: {}%, kappa: {}, rho: {}'.format(compnum, comp_var,
-                                                                         comp_kappa, comp_rho)
+        comp_var = "{0:.2f}".format(comptable.iloc[compnum]["variance explained"])
+        comp_kappa = "{0:.2f}".format(comptable.iloc[compnum]["kappa"])
+        comp_rho = "{0:.2f}".format(comptable.iloc[compnum]["rho"])
+        plt_title = 'Comp. {}: variance: {}%, kappa: {}, rho: {}, {}'.format(compnum, comp_var,
+                                                                             comp_kappa, comp_rho,
+                                                                             expl_text)
         title = ax_ts.set_title(plt_title)
         title.set_y(1.5)
 
@@ -127,7 +161,7 @@ def write_comp_figs(ts, mask, comptable, mmix, n_vols,
                     to_plot = ts_B[:, :, imgslice * cuts[idx], compnum]
 
                 ax_im = ax.imshow(to_plot, vmin=imgmin, vmax=imgmax, aspect='equal',
-                                  cmap='coolwarm')
+                                  cmap=png_cmap)
 
         # Add a color bar to the plot.
         ax_cbar = allplot.add_axes([0.8, 0.3, 0.03, 0.37])
@@ -144,17 +178,18 @@ def write_comp_figs(ts, mask, comptable, mmix, n_vols,
         ax_fft.plot(freqs, spectrum)
         ax_fft.set_title('One Sided fft')
         ax_fft.set_xlabel('Hz')
-        ax_fft.set_xbound(freqs[0], freqs[-1])
+        ax_fft.set_xlim(freqs[0], freqs[-1])
+        plt.yticks([])
 
         # Fix spacing so TR label does overlap with other plots
         allplot.subplots_adjust(hspace=0.4)
         plot_name = 'comp_{}.png'.format(str(compnum).zfill(3))
-        compplot_name = os.path.join('figures', plot_name)
+        compplot_name = os.path.join(out_dir, plot_name)
         plt.savefig(compplot_name)
         plt.close()
 
 
-def write_kappa_scatter(comptable):
+def write_kappa_scatter(comptable, out_dir):
     """
     Creates a scatter plot of Kappa vs Rho values. The shape and size of the
     points is based on classification and variance explained, respectively.
@@ -196,16 +231,18 @@ def write_kappa_scatter(comptable):
     ax_scatter.xaxis.label.set_fontsize(20)
     ax_scatter.yaxis.label.set_fontsize(20)
     ax_scatter.title.set_fontsize(25)
-    scatter_title = os.path.join('figures', 'Kappa_vs_Rho_Scatter.png')
+    scatter_title = os.path.join(out_dir, 'Kappa_vs_Rho_Scatter.png')
     plt.savefig(scatter_title)
 
     plt.close()
 
 
-def write_summary_fig(comptable):
+def write_summary_fig(comptable, out_dir):
     """
-    Creates a bar graph showing total variance explained by each component
-    as well as the number of components identified for each category.
+    Creates a pie chart showing 1) The total variance explained by each
+    component in the outer ring, 2) the variance explained by each
+    individual component in the inner ring, 3) counts of each classification
+    and 4) the amount of unexplained variance.
 
     Parameters
     ----------
@@ -215,38 +252,61 @@ def write_summary_fig(comptable):
         component, and (5) normalized variance explained by component
     """
 
-    # Get the variance and count of each classification
     var_expl = []
+    ind_var_expl = {}
     counts = {}
+    # Get overall variance explained, each components variance and counts of comps
     for clf in ['accepted', 'rejected', 'ignored']:
         var_expl.append(np.sum(comptable[comptable.classification == clf]['variance explained']))
+        ind_var_expl[clf] = comptable[comptable.classification == clf]['variance explained'].values
         counts[clf] = '{0} {1}'.format(comptable[comptable.classification == clf].count()[0], clf)
 
-    fig, ax = plt.subplots(figsize=(10, 7))
-    plt.bar([1, 2, 3], var_expl, color=['g', 'r', 'k'])
-    plt.xticks([1, 2, 3], counts.values(), fontsize=20)
-    plt.yticks(fontsize=15)
-    plt.ylabel('Variance Explained', fontsize=20)
-    plt.title('Component Overview', fontsize=25)
-    sumfig_title = os.path.join('figures', 'Component_Overview.png')
+    # Genereate Colormaps for individual components
+    acc_colors = plt.cm.Greens(np.linspace(0.2, .6, len(ind_var_expl['accepted'].tolist())))
+    rej_colors = plt.cm.Reds(np.linspace(0.2, .6, len(ind_var_expl['rejected'].tolist())))
+    ign_colors = plt.cm.Greys(np.linspace(0.2, .8, len(ind_var_expl['ignored'].tolist())))
+    unxp_colors = np.atleast_2d(np.array(plt.cm.Greys(0)))
+
+    # Shuffle the colors so that neighboring wedges are (perhaps) visually seperable
+    np.random.shuffle(rej_colors)
+    np.random.shuffle(acc_colors)
+    np.random.shuffle(ign_colors)
+
+    # Decision on whether to include the unexplained variance in figure
+    unexpl_var = [100 - np.sum(var_expl)]
+    all_var_expl = []
+    if unexpl_var >= [0.001]:
+        var_expl += unexpl_var
+        counts['unexplained'] = 'unexplained variance'
+        # Combine individual variances from giant list
+        for value in ind_var_expl.values():
+            all_var_expl += value.tolist()
+        # Add in unexplained variance
+        all_var_expl += unexpl_var
+        outer_colors = np.stack((plt.cm.Greens(0.7), plt.cm.Reds(0.7),
+                                 plt.cm.Greys(0.7), plt.cm.Greys(0)))
+        inner_colors = np.concatenate((acc_colors, rej_colors, ign_colors, unxp_colors), axis=0)
+    else:
+        for value in ind_var_expl.values():
+            all_var_expl += value.tolist()
+        outer_colors = np.stack((plt.cm.Greens(0.7), plt.cm.Reds(0.7), plt.cm.Greys(0.7)))
+        inner_colors = np.concatenate((acc_colors, rej_colors, ign_colors), axis=0)
+
+    labels = counts.values()
+
+    fig, ax = plt.subplots(figsize=(16, 10))
+    size = 0.3
+    # Build outer, overall pie chart, and then inner individual comp pie
+    ax.pie(var_expl, radius=1, colors=outer_colors, labels=labels,
+           autopct='%1.1f%%', pctdistance=0.85, textprops={'fontsize': 20},
+           wedgeprops=dict(width=size, edgecolor='w'))
+
+    ax.pie(all_var_expl, radius=1 - size, colors=inner_colors,
+           wedgeprops=dict(width=size))
+
+    ax.set(aspect="equal")
+    ax.set_title('Variance Explained By Classification', fontdict={'fontsize': 28})
+    if unexpl_var < [0.001]:
+        plt.text(1, -1, '*Unexplained Variance less than 0.001', fontdict={'fontsize': 12})
+    sumfig_title = os.path.join(out_dir, 'Component_Overview.png')
     plt.savefig(sumfig_title)
-
-
-def get_spectrum(data: np.array, tr: float = 1):
-    """
-    Returns the power spectrum and corresponding frequencies when provided
-    with a component time course and repitition time.
-
-    Parameters
-    ----------
-    data : (S, ) array_like
-            A timeseries S, on which you would like to perform an fft.
-    tr : :obj:`float`
-            Reptition time (TR) of the data
-    """
-
-    # adapted from @dangom
-    power_spectrum = np.abs(np.fft.rfft(data)) ** 2
-    freqs = np.fft.rfftfreq(power_spectrum.size * 2 - 1, tr)
-    idx = np.argsort(freqs)
-    return power_spectrum[idx], freqs[idx]

tedana/workflows/tedana.py

@@ -147,6 +147,11 @@ def _get_parser():
                                 'maps, timecourse plots and other diagnostic '
                                 'images'),
                           default=False)
+    optional.add_argument('--png-cmap',
+                          dest='png_cmap',
+                          type=str,
+                          help=('Colormap for figures'),
+                          default='coolwarm')
     optional.add_argument('--maxit',
                           dest='maxit',
                           type=int,
@@ -179,7 +184,7 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None,
                     tedort=False, gscontrol=None, tedpca='mle',
                     ste=-1, combmode='t2s', verbose=False, stabilize=False,
                     out_dir='.', fixed_seed=42, maxit=500, maxrestart=10,
-                    debug=False, quiet=False, png=False):
+                    debug=False, quiet=False, png=False, png_cmap='coolwarm'):
     """
     Run the "canonical" TE-Dependent ANAlysis workflow.
 
@@ -221,6 +226,9 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None,
         Generate intermediate and additional files. Default is False.
     png : obj:'bool', optional
         Generate simple plots and figures. Default is false.
+    png_cmap : obj:'str', optional
+            Name of a matplotlib colormap to be used when generating figures.
+            --png must still be used to request figures. Default is 'coolwarm'
     out_dir : :obj:`str`, optional
         Output directory.
 
@@ -422,15 +430,22 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None,
     if png:
         LGR.info('Making figures folder with static component maps and '
                  'timecourse plots.')
+        # make figure folder first
+        if not op.isdir(op.join(out_dir, 'figures')):
+            os.mkdir(op.join(out_dir, 'figures'))
+
         viz.write_comp_figs(data_oc, mask=mask, comptable=comptable, mmix=mmix,
-                            n_vols=n_vols, acc=acc, rej=rej, midk=midk,
-                            empty=ign, ref_img=ref_img)
+                            ref_img=ref_img,
+                            out_dir=op.join(out_dir, 'figures'),
+                            png_cmap=png_cmap)
 
         LGR.info('Making Kappa vs Rho scatter plot')
-        viz.write_kappa_scatter(comptable=comptable)
+        viz.write_kappa_scatter(comptable=comptable,
+                                out_dir=op.join(out_dir, 'figures'))
 
         LGR.info('Making overall summary figure')
-        viz.write_summary_fig(comptable=comptable)
+        viz.write_summary_fig(comptable=comptable,
+                              out_dir=op.join(out_dir, 'figures'))
 
     LGR.info('Workflow completed')
     for handler in logging.root.handlers[:]: