Multiple documentation updates

docs/approach.rst

@@ -36,6 +36,11 @@ manner.
 
 .. image:: /_static/a02_echo_value_distributions.png
 
+.. note::
+    In this example, the non-steady state volumes at the beginning of the run are
+    excluded. Some pulse sequences save these initial volumes and some do not. If
+    they are saved, then the first few volume in a run will have much larger relative
+    magnitudes. These initial volumes should be removed before running ``tedana``
 
 ************************
 Adaptive mask generation
@@ -49,20 +54,25 @@ have good signal for some echoes.
 In order to avoid using bad signal from affected echoes in calculating
 :math:`T_{2}^*` and :math:`S_{0}` for a given voxel, ``tedana`` generates an
 adaptive mask, where the value for each voxel is the number of echoes with
-"good" signal.
-When :math:`T_{2}^*` and :math:`S_{0}` are calculated below, each voxel's values
-are only calculated from the first :math:`n` echoes, where :math:`n` is the
-value for that voxel in the adaptive mask.
+"good" signal. The voxel in the shortest echo with the 33rd percentile mean signal
+across time is identified. The threshold for each echo is the signal in the same voxel
+divided by 3. This is an arbitrary, but conservative threshold in that it only
+excludes voxels where the signal is much lower than other measured signals in
+each echo. When :math:`T_{2}^*` and :math:`S_{0}` are calculated below, each
+voxel's valusesare only calculated from the first :math:`n` echoes, where
+:math:`n` is the value for that voxel in the adaptive mask. By default, the
+optimally combined and denoised time series will include voxels where there
+is at least one good echo, but ICA and the fit maps require at least three
+good echoes.
 
 .. note::
     ``tedana`` allows users to provide their own mask.
     The adaptive mask will be computed on this explicit mask, and may reduce
     it further based on the data.
     If a mask is not provided, ``tedana`` runs :func:`nilearn.masking.compute_epi_mask`
     on the first echo's data to derive a mask prior to adaptive masking.
-    The workflow does this because the adaptive mask generation function
-    sometimes identifies almost the entire bounding box as "brain", and
-    ``compute_epi_mask`` restricts analysis to a more reasonable area.
+    Some brain masking is required because the percentile-based thresholding
+    in the adaptive mask will be flawed if it includes all out-of-brain voxels.
 
 .. image:: /_static/a03_adaptive_mask.png
   :width: 600 px

docs/faq.rst

@@ -63,6 +63,25 @@ The standard space template in this example is "MNI152NLin2009cAsym", but will d
 
     <script src="https://gist.github.com/tsalo/f9f38e9aba901e99ddb720465bb5222b.js"></script>
 
+************************************************
+[tedana] Why are there so few voxels in my mask?
+************************************************
+
+Several steps in ``tedana`` require excluding voxels outside of the brain. If the
+provided data are not already masked, ``tedana`` runs
+:func:`nilearn.masking.compute_epi_mask` to create a mask. This often works, but
+it is a fairly simple method and there are better masking tools in most fMRI
+processing pipelines, such as `AFNI`_. If your full preprocessing pipeline already
+generates a mask for fMRI data, we recommend inputting that mask to ``tedana`` using
+the ``--mask`` option. ``tedana`` additionally creates an adaptive mask that will
+exclude voxels from some steps that have large dropout in later echoes. The adaptive
+mask will flag some voxels in common dropout regions, but should not radically alter
+the inputted mask. Here is more information on the
+`creation and use of the adaptive mask`_.
+
+.. _AFNI: http://afni.nimh.nih.gov
+.. _creation and use of the adaptive mask: approach.html#adaptive-mask-generation
+
 ************************************
 [tedana] ICA has failed to converge.
 ************************************

docs/multi-echo.rst

@@ -440,8 +440,25 @@ For single-echo EPI data, that excitation time would be the same regardless of t
 and the same is true when one is collecting multiple echoes after a single excitation pulse.
 Therefore, we suggest using the same slice timing for all echoes in an ME-EPI series.
 
+3. Apply spatial normalization and susceptibility distortion correction consistently across echoes
+==================================================================================================
 
-3. Perform distortion correction, spatial normalization, smoothing, and any rescaling or filtering **after** denoising
+One key feature of susceptibility distortion is that it is primarily a factor of readout pattern and
+total readout time, rather than echo time. This means that, for most multi-echo sequences, even though
+dropout will increase with echo time, distortion will not (at least not to a noticeable/meaningful extent).
+
+For this reason, if you are applying TOPUP-style (blip-up/blip-down) "field maps",
+we recommend using your first echo time, as this will exhibit the least dropout.
+If your first echo time is very short, and exhibits poor gray/white contrast, then a later echo time may
+be preferable. In any case, you should calculate the spatial transform from just one of your echoes and
+apply it across all of them.
+
+Similarly, if spatial normalization to a template is done, the spatial transform should be calculated
+once and the same transformation (ideally on transformation for motion correction, distortion correction,
+and spatial normalization) should be applied to all echoes.
+
+
+4. Perform smoothing, and any rescaling or filtering **after** denoising
 ======================================================================================================================
 
 Any step that will alter the relationship of signal magnitudes between echoes should occur after denoising and combining
@@ -452,32 +469,19 @@ and the subsequent calculation of voxelwise T2* values will be distorted or inco
 time point.
 
 .. note::
-    We are assuming that spatial normalization and distortion correction, particularly non-linear normalization methods
-    with higher order interpolation functions, are likely to distort the relationship between echoes while rigid body
-    motion correction would linearly alter each echo in a similar manner. This assumption has not yet been empirically
-    tested and an affine normalzation with bilinear interpolation may not distort the relationship between echoes.
-    Additionally, there are benefits to applying only one spatial transform to data rather than applying one spatial
-    transform for motion correction and a later transform for normalization and distortion correction. Our advice
-    against doing normalization and distortion correction is a conservative choice and we encourage additional
-    research to better understand how these steps can be applied before denoising.
-
-
-4. Apply susceptibility distortion correction consistently across echoes
-========================================================================
-
-One key feature of susceptibility distortion is that it is primarily a factor of readout pattern and total readout time, rather than echo time.
-This means that, for most multi-echo sequences, even though dropout will increase with echo time,
-distortion will not (at least not to a noticeable/meaningful extent).
-
-For this reason, if you are applying TOPUP-style (blip-up/blip-down) "field maps",
-we recommend using your first echo time, as this will exhibit the least dropout.
-If your first echo time is very short, and exhibits poor gray/white contrast, then a later echo time may be preferable.
-In any case, you should calculate the spatial transform from just one of your echoes and apply it across all of them.
-
+    Spatial normalization and distortion correction, particularly non-linear normalization methods
+    with higher order interpolation functions to regrid voxels in a new space, may distort the relationship between
+    echoes more than bilinear interpolation. This has the potential to distort the relationship between echoes and
+    there have been anecdotal cases where this might be an issue. Still, since serial spatial transforms spatially
+    smooth the data and most modern pipeline combine all spatial transforms into a single step, we recommend doing
+    these steps before running denoising. Particularly for data with high intensity heterogeneity between the surface
+    and center of the brain, we recommend checking if distoration correction and normalization adversely affect the
+    relationship between echoes.
 
 .. _fMRIPrep: https://fmriprep.readthedocs.io
 .. _afni_proc.py: https://afni.nimh.nih.gov/pub/dist/doc/program_help/afni_proc.py.html
 
+
 *****************
 General Resources
 *****************
@@ -542,6 +546,11 @@ For more details, see the `fmriprep workflows page`_ and :ref:`collecting fMRIPr
 
 .. _fmriprep workflows page: https://fmriprep.readthedocs.io/en/stable/workflows.html
 
+`fmrwhy`_ runs BIDS-compatible fMRI analysis with SPM12 and supports multi-echo data, 
+but it is no longer being actively maintained.
+
+.. _fmrwhy: https://fmrwhy.readthedocs.io
+
 Currently SPM and FSL do not natively support multi-echo fmri data processing.
 
 

docs/outputs.rst

@@ -314,22 +314,22 @@ analyses. This page describes the plots forming the reports and well as
 information on how to take advantage of the interactive functionalities.
 You can also play around with `our demo`_.
 
-.. _our demo: https://me-ica.github.io/tedana-ohbm-2020/
+.. _our demo: https://me-ica.github.io/ohbm-2023-multiecho/tedana/tedana_results_minimal_five-echo/tedana_report.html
 
 
 Report Structure
 ----------------
 
-The image below shows a representative report, which has two sections: a) the summary view,
-and b) the individual component view.
+The image below shows a representative report. The left is a summary view
+which contains information on all components and the right presents additional
+information for an individual component. One can hover over any pie chart wedge
+or data point in the summary view to see additional information about a
+component. Clicking on a component will select the component and the additional
+information will appear to the right.
 
 .. image:: /_static/rep01_overallview.png
   :align: center
 
-.. note::
-  When a report is initially loaded, as no component is selected on the
-  summary view, the individual component view appears empty.
-
 
 Summary View
 ------------
@@ -339,37 +339,46 @@ selection results. It includes four different plots.
 
 * **Kappa/Rho Scatter:** This is a scatter plot of `Kappa` vs. `Rho` features for all components.
   In the plot, each dot represents a different component. The x-axis represents the kappa feature, and the
-  y-axis represents the rho feature. These are two of the most
-  informative features describing the likelihood of the component
-  being BOLD or non-BOLD. Additional information is provided via color
-  and size. In particular, color informs about its classification
-  status (e.g., accepted, rejected); while size relates to
-  the amount of variance explained by the component (larger dot,
-  larger variance).
+  y-axis represents the rho feature. `Kappa` is a summary metric on for how much
+  BOLD information might be in a component and `rho` is a summary metric for how
+  much non-BOLD information is in a component. Thus a component with a higher `kappa`
+  and lower `rho` value is more likely to be retained. The solid gray line is
+  :math:`\kappa=\rho`.  Color is used to label accepted (green) or rejected (red)
+  components. The size of the circle is the amount of variance explained by the
+  component so larger circle (higher variance) that seems misclassified is worth
+  closer inspection. The component classification process uses kappa and rho elbow
+  thresholds (black dashed lines) along with other criteria. Most accepted
+  components should be greater than the kappa elbow and less than the rho elbow.
+  Accepted or rejected components that don't cross those thresholds might be
+  worth additional inspection. Hovering over a component also shows a `Tag`
+  that explains why a component received its classification.
 
 .. image:: /_static/rep01_kapparhoScatter.png
   :align: center
   :height: 400px
 
 * **Kappa Scree Plot:** This scree plot provides a view of the components ranked by `Kappa`.
   As in the previous plot, each dot represents a component. The color of the dot informs us
-  about classification status. In this plot, size is not related to variance explained.
+  about classification status. The dashed line is the kappa elbow threshold.
+  In this plot, size is not related to variance explained, but you can see the variance
+  explained by hovering over any dot.
 
 .. image:: /_static/rep01_kappaScree.png
   :align: center
   :height: 400px
 
 * **Rho Scree Plot:** This scree plot provides a view of the components ranked by `Rho`.
   As in the previous plot, each dot represents a component. The color of the dot informs us
-  about classification status. In this plot, size is not related to variance explained.
+  about classification status. The dashed line is the rho elbow threshold.
+  Size is not related to variance explained.
 
 .. image:: /_static/rep01_rhoScree.png
   :align: center
   :height: 400px
 
 * **Variance Explained Plot:** This pie plot provides a summary of how much variance is explained
-  by each individual component, as well as the total variance explained by each of the three
-  classification categories (i.e., accepted, rejected, ignored). In this plot, each component is
+  by each individual component, as well as the total variance explained by each of the two
+  classification categories (i.e., accepted, rejected). In this plot, each component is
   represented as a wedge, whose size is directly related to the amount of variance explained. The
   color of the wedge inform us about the classification status of the component. For this view,
   components are sorted by classification first, and inside each classification group by variance
@@ -387,25 +396,41 @@ This view provides detailed information about an individual
 component (selected in the summary view, see below). It includes three different plots.
 
 * **Time series:** This plot shows the time series associated with a given component
-  (selected in the summary view). The x-axis represents time (in units of TR), and the
-  y-axis represents signal levels (in arbitrary units). Finally, the color of the trace
-  informs us about the component classification status.
+  (selected in the summary view). The x-axis represents time (in units of TR and seconds),
+  and the y-axis represents signal levels (in arbitrary units).
+  Finally, the color of the trace informs us about the component classification status.
+  Plausibly BOLD-weighted components might have reponses the follow the a task design, 
+  while components that are less likely to be BOLD-weighted might have large signal
+  spikes or slow drifts. If a high variance component time series initially has a few
+  very high magnitude volumes, that is a sign non-steady state volumes were not removed
+  before running ``tedana``. Keeping these volumes might results in a suboptimal ICA
+  results. ``tedana`` should be run without any initial non-steady state volumes.
 
 .. image:: /_static/rep01_tsPlot.png
   :align: center
   :height: 150px
 
 * **Component beta map:** This plot shows the map of the beta coefficients associated with
   a given component (selected in the summary view). The colorbar represents the amplitude
-  of the beta coefficients.
+  of the beta coefficients. The same weights could be flipped postive/negative so relative
+  values are more relevant that what is very positive vs negative.
+  Plausibly BOLD-weighted components should have larger hotspots in area that follow
+  cortical or cerebellar brain structure. Hotspots in ventricles, on the edges of the
+  brain or slice-specific or slice-alternating effects are signs of artifacts.
 
 .. image:: /_static/rep01_betaMap.png
   :align: center
   :height: 400px
 
 * **Spectrum:** This plot shows the spectrogram associated with a given component
   (selected in the summary view). The x-axis represents frequency (in Hz), and the
-  y-axis represents spectral amplitude.
+  y-axis represents spectral amplitude. BOLD-weighted signals will likely have most
+  power below 0.1Hz. Peaks at higher frequencies are signs of non-BOLD signals. A
+  respiration artifact might be around 0.25-0.33Hz and a cardiac artifact might be
+  around 1Hz. This plot shows the maximum resolvable frequency given the TR, so
+  those higher frequencies might fold over to different peaks that are still above
+  0.1Hz. Respirator and cardiac fluctuation artifacts are also sometimes visible
+  in the time series.
 
 .. image:: /_static/rep01_fftPlot.png
   :align: center
@@ -424,6 +449,7 @@ figures.
   :align: center
   :height: 25px
 
+
 The table below includes information about all available interactions
 
 .. |Reset| image:: /_static/rep01_tool_reset.png
@@ -490,9 +516,11 @@ Carpet plots
 ************
 
 In additional to the elements described above, ``tedana``'s interactive reports include carpet plots for the main outputs of the workflow:
-the optimally combined data, the denoised data, the high-Kappa (accepted) data, and the low-Kappa (rejected) data.
-
-These plots may be useful for visual quality control of the overall denoising run.
+the optimally combined data, the denoised data, the high-Kappa (accepted) data, and the low-Kappa (rejected) data. Each row is a voxel
+and the grayscale is the relative signal changes across time. After denoising, voxels that look very different from others across time
+or time points that are uniformly high or low across voxels are concerning. These carpet plots can be help as a quick quality check for
+data, but since some neural signals really are more global than others and there are voxelwise differences in responses, quality checks
+should not overly focus on carpet plots and should examine these results in context with other quality measures.
 
 .. image:: /_static/carpet_overview.png
   :align: center