diff --git a/CHANGES.rst b/CHANGES.rst
index 14e6bb0d47..efb1faea83 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -10,6 +10,7 @@ Cubeviz
 
 - Added functionality to Collapse and Spectral Extraction plugins to save results to FITS file. [#2586]
 
+- Moment map plugin now supports linear per-spaxel continuum subtraction. [#2587]
 
 Imviz
 ^^^^^
@@ -26,6 +27,9 @@ Specviz2d
 API Changes
 -----------
 
+- ``width`` argument in Line Analysis plugin is renamed to ``continuum_width`` and ``width``
+  will be removed in a future release. [#2587]
+
 Cubeviz
 ^^^^^^^
 
diff --git a/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py b/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py
index 3feb986787..c52381f93a 100644
--- a/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py
+++ b/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py
@@ -17,6 +17,8 @@
                                         SpectralSubsetSelectMixin,
                                         AddResultsMixin,
                                         SelectPluginComponent,
+                                        SpectralContinuumMixin,
+                                        skip_if_no_updates_since_last_active,
                                         with_spinner)
 from jdaviz.core.user_api import PluginUserApi
 
@@ -30,7 +32,7 @@
 @tray_registry('cubeviz-moment-maps', label="Moment Maps",
                viewer_requirements=['spectrum', 'image'])
 class MomentMap(PluginTemplateMixin, DatasetSelectMixin, SpectralSubsetSelectMixin,
-                AddResultsMixin):
+                SpectralContinuumMixin, AddResultsMixin):
     """
     See the :ref:`Moment Maps Plugin Documentation <moment-maps>` for more details.
 
@@ -44,6 +46,14 @@ class MomentMap(PluginTemplateMixin, DatasetSelectMixin, SpectralSubsetSelectMix
       Dataset to use for computing line statistics.
     * ``spectral_subset`` (:class:`~jdaviz.core.template_mixin.SubsetSelect`):
       Subset to use for the line, or ``Entire Spectrum``.
+    * ``continuum`` (:class:`~jdaviz.core.template_mixin.SubsetSelect`):
+      Subset to use for the continuum, or ``None`` to skip continuum subtraction,
+      or ``Surrounding`` to use a region surrounding the
+      subset set in ``spectral_subset``.
+    * ``continuum_width``:
+      Width, relative to the overall line spectral region, to fit the linear continuum
+      (excluding the region containing the line). If 1, will use endpoints within line region
+      only.
     * ``n_moment``
     * ``output_unit``
       Choice of "Wavelength" or "Velocity", applicable for n_moment >= 1.
@@ -53,6 +63,7 @@ class MomentMap(PluginTemplateMixin, DatasetSelectMixin, SpectralSubsetSelectMix
     * :meth:`calculate_moment`
     """
     template_file = __file__, "moment_maps.vue"
+    uses_active_status = Bool(True).tag(sync=True)
 
     n_moment = IntHandleEmpty(0).tag(sync=True)
     filename = Unicode().tag(sync=True)
@@ -88,20 +99,32 @@ def __init__(self, *args, **kwargs):
 
     @property
     def _default_image_viewer_reference_name(self):
-        return self.jdaviz_helper._default_image_viewer_reference_name
+        return getattr(
+            self.app._jdaviz_helper, '_default_spectrum_viewer_reference_name', 'flux-viewer'
+        )
 
     @property
     def _default_spectrum_viewer_reference_name(self):
-        return self.jdaviz_helper._default_spectrum_viewer_reference_name
+        return getattr(
+            self.app._jdaviz_helper, '_default_spectrum_viewer_reference_name', 'spectrum-viewer'
+        )
 
     @property
     def user_api(self):
         # NOTE: leaving save_as_fits out for now - we may want a more general API to do that
         # accross all plugins at some point
-        return PluginUserApi(self, expose=('dataset', 'spectral_subset', 'n_moment',
+        return PluginUserApi(self, expose=('dataset', 'spectral_subset',
+                                           'continuum', 'continuum_width',
+                                           'n_moment',
                                            'output_unit', 'reference_wavelength',
                                            'add_results', 'calculate_moment'))
 
+    @observe('is_active')
+    def _is_active_changed(self, msg):
+        for pos, mark in self.continuum_marks.items():
+            mark.visible = self.is_active
+        self._calculate_continuum(msg)
+
     @observe("dataset_selected", "dataset_items", "n_moment")
     def _set_default_results_label(self, event={}):
         label_comps = []
@@ -120,6 +143,21 @@ def _set_data_units(self, event={}):
             else:
                 self.dataset_spectral_unit = ""
 
+    @observe("dataset_selected", "spectral_subset_selected",
+             "continuum_subset_selected", "continuum_width")
+    @skip_if_no_updates_since_last_active()
+    def _calculate_continuum(self, msg=None):
+        if not hasattr(self, 'dataset') or self.app._jdaviz_helper is None:  # noqa
+            # during initial init, this can trigger before the component is initialized
+            return
+
+        # NOTE: there is no use in caching this, as the continuum will need to be re-computed
+        # per-spaxel to use within calculating the moment map
+        _ = self._get_continuum(self.dataset,
+                                None,
+                                self.spectral_subset,
+                                update_marks=True)
+
     @with_spinner()
     def calculate_moment(self, add_data=True):
         """
@@ -130,12 +168,19 @@ def calculate_moment(self, add_data=True):
         add_data : bool
             Whether to add the resulting data object to the app according to ``add_results``.
         """
-        # Retrieve the data cube and slice out desired region, if specified
-        if "_orig_spec" in self.dataset.selected_obj.meta:
-            cube = self.dataset.selected_obj.meta["_orig_spec"]
+        if self.continuum.selected == 'None':
+            if "_orig_spec" in self.dataset.selected_obj.meta:
+                cube = self.dataset.selected_obj.meta["_orig_spec"]
+            else:
+                cube = self.dataset.selected_obj
         else:
-            cube = self.dataset.selected_obj
+            _, _, cube = self._get_continuum(self.dataset,
+                                             'per-pixel',
+                                             self.spectral_subset,
+                                             update_marks=False)
 
+        # slice out desired region
+        # TODO: should we add a warning for a composite spectral subset?
         spec_min, spec_max = self.spectral_subset.selected_min_max(cube)
         slab = manipulation.spectral_slab(cube, spec_min, spec_max)
 
diff --git a/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.vue b/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.vue
index 22cf667b09..7c230ac47b 100644
--- a/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.vue
+++ b/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.vue
@@ -2,8 +2,16 @@
   <j-tray-plugin
     description='Create a 2D image from a data cube.'
     :link="docs_link || 'https://jdaviz.readthedocs.io/en/'+vdocs+'/'+config+'/plugins.html#moment-maps'"
+    :uses_active_status="uses_active_status"
+    @plugin-ping="plugin_ping($event)"
+    :keep_active.sync="keep_active"
     :popout_button="popout_button">
 
+    <j-plugin-section-header>Cube</j-plugin-section-header>
+    <v-row>
+      <j-docs-link>Choose the input cube and spectral subset.</j-docs-link>
+    </v-row>
+
     <plugin-dataset-select
       :items="dataset_items"
       :selected.sync="dataset_selected"
@@ -22,6 +30,50 @@
       hint="Spectral region to compute the moment map."
     />
 
+    <j-plugin-section-header>Continuum Subtraction</j-plugin-section-header>
+    <v-row>
+      <j-docs-link v-if="continuum_subset_selected==='None'">
+        Choose whether and how to compute the continuum for continuum subtraction.
+      </j-docs-link>
+      <j-docs-link v-else>
+        {{continuum_subset_selected=='Surrounding' && spectral_subset_selected=='Entire Spectrum' ? "Since using the entire spectrum, the end points will be used to fit a linear continuum." : "Choose a region to fit a linear line as the underlying continuum."}}  
+        {{continuum_subset_selected=='Surrounding' && spectral_subset_selected!='Entire Spectrum' ? "Choose a width in number of data points to consider on each side of the line region defined above." : null}}
+        When this plugin is opened, a visual indicator will show on the spectrum plot showing the continuum fitted as a thick line, and interpolated into the line region as a thin line.
+        When computing the moment map, these same input options will be used to compute and subtract a linear continuum for each spaxel, independently.
+      </j-docs-link>
+    </v-row>
+
+    <plugin-subset-select 
+      :items="continuum_subset_items"
+      :selected.sync="continuum_subset_selected"
+      :show_if_single_entry="true"
+      :rules="[() => continuum_subset_selected!==spectral_subset_selected || 'Must not match line selection.']"
+      label="Continuum"
+      hint="Select spectral region that defines the continuum."
+    />
+
+    <v-row v-if="continuum_subset_selected=='Surrounding' && spectral_subset_selected!='Entire Spectrum'">
+      <!-- DEV NOTE: if changing the validation rules below, also update the logic to clear the results
+           in line_analysis.py  -->
+      <v-text-field
+        label="Width"
+        type="number"
+        v-model.number="continuum_width"
+        step="0.1"
+        :rules="[() => continuum_width!=='' || 'This field is required.',
+                 () => continuum_width<=10 || 'Width must be <= 10.',
+                 () => continuum_width>=1 || 'Width must be >= 1.']"
+        hint="Width, relative to the overall line spectral region, to fit the linear continuum (excluding the region containing the line).  If 1, will use endpoints within line region only."
+        persistent-hint
+      >
+      </v-text-field>
+    </v-row>
+
+    <j-plugin-section-header>Moment</j-plugin-section-header>
+    <v-row>
+      <j-docs-link>Options for generating the moment map.</j-docs-link>
+    </v-row>
+
     <v-row>
       <v-text-field
         ref="n_moment"
diff --git a/jdaviz/configs/cubeviz/plugins/moment_maps/tests/test_moment_maps.py b/jdaviz/configs/cubeviz/plugins/moment_maps/tests/test_moment_maps.py
index 7b512b71e1..26374790ee 100644
--- a/jdaviz/configs/cubeviz/plugins/moment_maps/tests/test_moment_maps.py
+++ b/jdaviz/configs/cubeviz/plugins/moment_maps/tests/test_moment_maps.py
@@ -12,6 +12,30 @@
 from jdaviz.configs.cubeviz.plugins.moment_maps.moment_maps import MomentMap
 
 
+def test_user_api(cubeviz_helper, spectrum1d_cube):
+    with warnings.catch_warnings():
+        warnings.filterwarnings("ignore", message="No observer defined on WCS.*")
+        cubeviz_helper.load_data(spectrum1d_cube, data_label='test')
+
+    mm = cubeviz_helper.plugins['Moment Maps']
+    assert not mm._obj.continuum_marks['center'].visible
+    with mm.as_active():
+        assert mm._obj.continuum_marks['center'].visible
+        mm.n_moment = 0
+        # no continuum so marks should be empty
+        assert len(mm._obj.continuum_marks['center'].x) == 0
+
+        mom = mm.calculate_moment()
+
+        mm.continuum = 'Surrounding'
+        mm.continuum_width = 10
+        assert len(mm._obj.continuum_marks['center'].x) > 0
+
+        mom_sub = mm.calculate_moment()
+
+        assert mom != mom_sub
+
+
 def test_moment_calculation(cubeviz_helper, spectrum1d_cube, tmpdir):
     dc = cubeviz_helper.app.data_collection
     with warnings.catch_warnings():
diff --git a/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.py b/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.py
index 1c4d340abb..4d1d772214 100644
--- a/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.py
+++ b/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.py
@@ -1,5 +1,5 @@
 import os
-import warnings
+import logging
 
 import numpy as np
 from glue.core.message import (SubsetDeleteMessage,
@@ -8,9 +8,7 @@
 from traitlets import Bool, List, Float, Unicode, observe
 from astropy import units as u
 from specutils import analysis, Spectrum1D
-from specutils.manipulation import extract_region
 
-from jdaviz.core.custom_traitlets import FloatHandleEmpty
 from jdaviz.core.events import (AddDataMessage,
                                 RemoveDataMessage,
                                 SpectralMarksChangedMessage,
@@ -18,19 +16,16 @@
                                 RedshiftMessage,
                                 GlobalDisplayUnitChanged,
                                 SnackbarMessage)
-from jdaviz.core.marks import (LineAnalysisContinuum,
-                               LineAnalysisContinuumCenter,
-                               LineAnalysisContinuumLeft,
-                               LineAnalysisContinuumRight,
-                               ShadowLine)
 from jdaviz.core.registries import tray_registry
 from jdaviz.core.template_mixin import (PluginTemplateMixin,
                                         DatasetSelectMixin,
                                         SpectralSubsetSelectMixin,
                                         DatasetSpectralSubsetValidMixin,
                                         SubsetSelect,
+                                        SpectralContinuumMixin,
                                         SPATIAL_DEFAULT_TEXT,
-                                        skip_if_no_updates_since_last_active)
+                                        skip_if_no_updates_since_last_active,
+                                        with_spinner)
 from jdaviz.core.user_api import PluginUserApi
 from jdaviz.core.tools import ICON_DIR
 
@@ -70,7 +65,7 @@ def _coerce_unit(quantity):
 
 @tray_registry('specviz-line-analysis', label="Line Analysis", viewer_requirements='spectrum')
 class LineAnalysis(PluginTemplateMixin, DatasetSelectMixin, SpectralSubsetSelectMixin,
-                   DatasetSpectralSubsetValidMixin):
+                   DatasetSpectralSubsetValidMixin, SpectralContinuumMixin):
     """
     The Line Analysis plugin returns specutils analysis for a single spectral line.
     See the :ref:`Line Analysis Plugin Documentation <line-analysis>` for more details.
@@ -90,7 +85,7 @@ class LineAnalysis(PluginTemplateMixin, DatasetSelectMixin, SpectralSubsetSelect
     * ``continuum`` (:class:`~jdaviz.core.template_mixin.SubsetSelect`):
       Subset to use for the continuum, or ``Surrounding`` to use a region surrounding the
       subset set in ``spectral_subset``.
-    * :attr:`width`:
+    * ```continuum_width``:
       Width, relative to the overall line spectral region, to fit the linear continuum
       (excluding the region containing the line). If 1, will use endpoints within line region
       only.
@@ -104,10 +99,6 @@ class LineAnalysis(PluginTemplateMixin, DatasetSelectMixin, SpectralSubsetSelect
     spatial_subset_items = List().tag(sync=True)
     spatial_subset_selected = Unicode().tag(sync=True)
 
-    continuum_subset_items = List().tag(sync=True)
-    continuum_subset_selected = Unicode().tag(sync=True)
-
-    width = FloatHandleEmpty(3).tag(sync=True)
     results_computing = Bool(False).tag(sync=True)
     results = List().tag(sync=True)
     results_centroid = Float().tag(sync=True)  # stored in AA units
@@ -125,17 +116,14 @@ def __init__(self, *args, **kwargs):
 
         self.update_results(None)
 
-        self.continuum = SubsetSelect(self,
-                                      'continuum_subset_items',
-                                      'continuum_subset_selected',
-                                      default_text='Surrounding',
-                                      filters=['is_spectral'])
-
         # when accessing the selected data, access the spectrum-viewer version
         self.dataset._viewers = [self._default_spectrum_viewer_reference_name]
         # require entries to be in spectrum-viewer (not other cubeviz images, etc)
         self.dataset.add_filter('layer_in_spectrum_viewer')
 
+        # continuum selection is mandatory for line-analysis
+        self._continuum_remove_none_option()
+
         if self.app.state.settings.get("configuration") == "cubeviz":
             self.spatial_subset = SubsetSelect(self,
                                                'spatial_subset_items',
@@ -167,10 +155,23 @@ def _default_spectrum_viewer_reference_name(self):
             self.app._jdaviz_helper, '_default_spectrum_viewer_reference_name', 'spectrum-viewer'
         )
 
+    # backwards compatibility for width (replace with user API deprecation)
+    @property
+    def width(self):
+        logging.warning(f"DeprecationWarning: width was replaced by continuum_width in 3.9 and will be removed in a future release")  # noqa
+        return self.continuum_width
+
+    @width.setter
+    def width(self, width):
+        logging.warning("DeprecationWarning: width was replaced by continuum_width in 3.9 and will be removed in a future release")  # noqa
+        self.continuum_width = width
+
     @property
     def user_api(self):
+        # deprecated: width was replaced with continuum_width in 3.9 so should be removed from the
+        # user API and the property and setter above as soon as 3.11.
         return PluginUserApi(self, expose=('dataset', 'spatial_subset', 'spectral_subset',
-                                           'continuum', 'width', 'get_results'))
+                                           'continuum', 'width', 'continuum_width', 'get_results'))
 
     def _on_viewer_data_changed(self, msg):
         viewer_id = self.app._viewer_item_by_reference(
@@ -220,52 +221,17 @@ def _is_active_changed(self, msg):
         if self.disabled_msg:
             return
 
-        for pos, mark in self.marks.items():
+        for pos, mark in self.continuum_marks.items():
             mark.visible = self.is_active
         self._calculate_statistics(msg)
 
-    @property
-    def marks(self):
-        marks = {}
-        viewer = self.app.get_viewer(self._default_spectrum_viewer_reference_name)
-        if viewer is None:
-            return {}
-        for mark in viewer.figure.marks:
-            if isinstance(mark, LineAnalysisContinuum):
-                # NOTE: we don't use isinstance anymore because of nested inheritance
-                if mark.__class__.__name__ == 'LineAnalysisContinuumLeft':
-                    marks['left'] = mark
-                elif mark.__class__.__name__ == 'LineAnalysisContinuumCenter':
-                    marks['center'] = mark
-                elif mark.__class__.__name__ == 'LineAnalysisContinuumRight':
-                    marks['right'] = mark
-
-        if not len(marks):
-            if not viewer.state.reference_data:
-                # we don't have data yet for scales, defer initializing
-                return {}
-            # then haven't been initialized yet, so initialize with empty
-            # marks that will be populated once the first analysis is done.
-            marks = {'left': LineAnalysisContinuumLeft(viewer, visible=self.is_active),
-                     'center': LineAnalysisContinuumCenter(viewer, visible=self.is_active),
-                     'right': LineAnalysisContinuumRight(viewer, visible=self.is_active)}
-            shadows = [ShadowLine(mark, shadow_width=2) for mark in marks.values()]
-            # NOTE: += won't trigger the figure to notice new marks
-            viewer.figure.marks = viewer.figure.marks + shadows + list(marks.values())
-
-        return marks
-
-    def update_results(self, results=None, mark_x={}, mark_y={}):
-        for pos, mark in self.marks.items():
-            mark.update_xy(mark_x.get(pos, []), mark_y.get(pos, []))
-
+    def update_results(self, results=None):
         if results is None:
             self.results = [{'function': function, 'result': ''} for function in FUNCTIONS]
+            self._update_continuum_marks()
         else:
             self.results = results
 
-        self.results_computing = False
-
     def get_results(self):
         # user-facing API call to force updating and retrieving results, even if the plugin
         # is closed
@@ -293,17 +259,19 @@ def _on_identified_line_changed(self, msg):
             # in which case we'll default to the identified line
             self.selected_line = self.identified_line
 
-    @observe("spatial_subset_selected", "spectral_subset_selected", "dataset_selected",
-             "continuum_subset_selected", "width")
+    @observe("dataset_selected", "spatial_subset_selected", "spectral_subset_selected",
+             "continuum_subset_selected", "continuum_width")
     @skip_if_no_updates_since_last_active()
+    @with_spinner('results_computing')
     def _calculate_statistics(self, msg={}):
         """
         Run the line analysis functions on the selected data/subset and
         display the results.
         """
-        if not hasattr(self, 'dataset') or self.app._jdaviz_helper is None or self.dataset_selected == '':  # noqa
+        if not hasattr(self, 'dataset') or self.app._jdaviz_helper is None:  # noqa
             # during initial init, this can trigger before the component is initialized
             return
+
         if self.disabled_msg != '':
             self.update_results(None)
             return
@@ -315,114 +283,14 @@ def _calculate_statistics(self, msg={}):
             self.update_results(None)
             return
 
-        # show spinner with overlay
-        self.results_computing = True
-        full_spectrum = self.dataset.selected_spectrum_for_spatial_subset(self.spatial_subset_selected,  # noqa
-                                                                          use_display_units=True)
-
-        if full_spectrum is None or self.width == "":
-            # this can happen DURING a unit conversion change
+        spectrum, continuum, spec_subtracted = self._get_continuum(self.dataset,
+                                                                   self.spatial_subset,
+                                                                   self.spectral_subset,
+                                                                   update_marks=True)
+        if spectrum is None:
             self.update_results(None)
             return
 
-        spectral_axis = full_spectrum.spectral_axis
-        if spectral_axis.unit == u.pix:
-            # plugin should be disabled so not get this far, but can still get here
-            # before the disabled message is set
-            self.update_results(None)
-            return
-
-        if self.continuum_subset_selected == self.spectral_subset_selected:
-            # already raised a validation error in the UI
-            self.update_results(None)
-            return
-
-        if self.spectral_subset_selected == "Entire Spectrum":
-            spectrum = full_spectrum
-        else:
-            sr = self.app.get_subsets(self.spectral_subset_selected,
-                                      simplify_spectral=True, use_display_units=True)
-            spectrum = extract_region(full_spectrum, sr, return_single_spectrum=True)
-            sr_lower = np.nanmin(spectrum.spectral_axis[spectrum.spectral_axis.value >= sr.lower.value])  # noqa
-            sr_upper = np.nanmax(spectrum.spectral_axis[spectrum.spectral_axis.value <= sr.upper.value])  # noqa
-
-        # compute continuum
-        if self.continuum_subset_selected == "Surrounding" and self.spectral_subset_selected == "Entire Spectrum": # noqa
-            # we know we'll just use the endpoints, so let's be efficient and not even
-            # try extracting from the region
-            continuum_mask = np.array([0, len(spectral_axis)-1])
-            mark_x = {'left': np.array([]),
-                      'center': np.array([min(spectral_axis.value), max(spectral_axis.value)]),
-                      'right': np.array([])}
-
-        elif self.continuum_subset_selected == "Surrounding":
-            # self.spectral_subset_selected != "Entire Spectrum"
-            if self.width > 10 or self.width < 1:
-                # DEV NOTE: if changing the limits, make sure to also update the form validation
-                # rules in line_analysis.vue
-                self.update_results(None)
-                return
-
-            spectral_region_width = sr_upper - sr_lower
-            # convert width from total relative width, to width per "side"
-            width = (self.width - 1) / 2
-            left, = np.where((spectral_axis < sr_lower) &
-                             (spectral_axis > sr_lower - spectral_region_width*width))
-            if not len(left):
-                # then no points matching the width are available outside the line region,
-                # so we'll default to the left-most point of the line region.
-                left, = np.where(spectral_axis == min(spectrum.spectral_axis))
-
-            right, = np.where((spectral_axis > sr_upper) &
-                              (spectral_axis < sr_upper + spectral_region_width*width))
-            if not len(right):
-                # then no points matching the width are available outside the line region,
-                # so we'll default to the right-most point of the line region.
-                right, = np.where(spectral_axis == max(spectrum.spectral_axis))
-
-            continuum_mask = np.concatenate((left, right))
-            mark_x = {'left': np.array([min(spectral_axis.value[continuum_mask]), sr_lower.value]),
-                      'center': np.array([sr_lower.value, sr_upper.value]),
-                      'right': np.array([sr_upper.value, max(spectral_axis.value[continuum_mask])])}
-
-        else:
-            # we'll access the mask of the continuum and then apply that to the spectrum.  For a
-            # spatially-collapsed spectrum in cubeviz, this will access the mask from the full
-            # cube, but still apply that to the spatially-collapsed spectrum.
-            continuum_mask = ~self._specviz_helper.get_data(
-                self.dataset.selected,
-                spectral_subset=self.continuum_subset_selected,
-                use_display_units=False).mask
-            spectral_axis_nanmasked = spectral_axis.value.copy()
-            spectral_axis_nanmasked[~continuum_mask] = np.nan
-            if self.spectral_subset_selected == "Entire Spectrum":
-                mark_x = {'left': spectral_axis_nanmasked,
-                          'center': spectral_axis.value,
-                          'right': []}
-            else:
-                mark_x = {'left': spectral_axis_nanmasked[spectral_axis.value < sr_lower.value],
-                          'right': spectral_axis_nanmasked[spectral_axis.value > sr_upper.value]}
-                # Center should extend (at least) across the line region between the full
-                # range defined by the continuum subset(s).
-                # OK for mark_x to be all NaNs.
-                with warnings.catch_warnings():
-                    warnings.simplefilter('ignore', category=RuntimeWarning)
-                    mark_x_min = np.nanmin(mark_x['left'])
-                    mark_x_max = np.nanmax(mark_x['right'])
-                left_min = np.nanmin([mark_x_min, sr_lower.value])
-                right_max = np.nanmax([mark_x_max, sr_upper.value])
-                mark_x['center'] = np.array([left_min, right_max])
-
-        continuum_x = spectral_axis[continuum_mask].value
-        min_x = min(spectral_axis.value)
-        continuum_y = full_spectrum.flux[continuum_mask].value
-        # DEV NOTE: could replace this with internal calls to the model fitting infrastructure
-        # to enable other model-types and to give visual feedback (by labeling the model
-        # as line_analysis:continuum, for example)
-        slope, intercept = np.polyfit(continuum_x-min_x, continuum_y, deg=1)
-        continuum = slope * (spectrum.spectral_axis.value-min_x) + intercept
-        mark_y = {k: slope * (v-min_x) + intercept for k, v in mark_x.items()}
-
         def _uncertainty(result):
             if getattr(result, 'uncertainty', None) is not None:
                 # we'll keep the uncertainty and result in the same unit (so
@@ -434,7 +302,6 @@ def _uncertainty(result):
                 return ''
 
         temp_results = []
-        spec_subtracted = spectrum - continuum
 
         if spec_subtracted.mask is not None:
             # temporary fix while mask may contain None:
@@ -553,7 +420,7 @@ def _uncertainty(result):
             # default to the identified line
             self.selected_line = self.identified_line
 
-        self.update_results(temp_results, mark_x, mark_y)
+        self.update_results(temp_results)
 
     def _compute_redshift_for_selected_line(self):
         index = self.line_items.index(self.selected_line)
diff --git a/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.vue b/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.vue
index 087ad9171e..f40d84da40 100644
--- a/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.vue
+++ b/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.vue
@@ -49,7 +49,7 @@
     <j-plugin-section-header>Continuum</j-plugin-section-header>
     <v-row>
       <j-docs-link>
-        {{spectral_subset_selected=='Entire Spectrum' ? "Since using the entire spectrum, the end points will be used to fit a linear continuum." : "Choose a region to fit a linear line as the underlying continuum."}}  
+        {{continuum_subset_selected=='Surrounding' && spectral_subset_selected=='Entire Spectrum' ? "Since using the entire spectrum, the end points will be used to fit a linear continuum." : "Choose a region to fit a linear line as the underlying continuum."}}  
         {{continuum_subset_selected=='Surrounding' && spectral_subset_selected!='Entire Spectrum' ? "Choose a width in number of data points to consider on each side of the line region defined above." : null}}
         When this plugin is opened, a visual indicator will show on the spectrum plot showing the continuum fitted as a thick line, and interpolated into the line region as a thin line.
       </j-docs-link>
@@ -70,11 +70,11 @@
       <v-text-field
         label="Width"
         type="number"
-        v-model.number="width"
+        v-model.number="continuum_width"
         step="0.1"
-        :rules="[() => width!=='' || 'This field is required.',
-                 () => width<=10 || 'Width must be <= 10.',
-                 () => width>=1 || 'Width must be >= 1.']"
+        :rules="[() => continuum_width!=='' || 'This field is required.',
+                 () => continuum_width<=10 || 'Width must be <= 10.',
+                 () => continuum_width>=1 || 'Width must be >= 1.']"
         hint="Width, relative to the overall line spectral region, to fit the linear continuum (excluding the region containing the line).  If 1, will use endpoints within line region only."
         persistent-hint
       >
diff --git a/jdaviz/core/template_mixin.py b/jdaviz/core/template_mixin.py
index 303ccb2eb5..5271143b3e 100644
--- a/jdaviz/core/template_mixin.py
+++ b/jdaviz/core/template_mixin.py
@@ -9,6 +9,7 @@
 import os
 import threading
 import time
+import warnings
 
 from echo import delay_callback
 from functools import cached_property
@@ -28,6 +29,7 @@
 from glue_jupyter.registries import viewer_registry
 from glue_jupyter.widgets.linked_dropdown import get_choices as _get_glue_choices
 from specutils import Spectrum1D
+from specutils.manipulation import extract_region
 from traitlets import Any, Bool, HasTraits, List, Unicode, observe
 
 from ipywidgets import widget_serialization
@@ -35,10 +37,16 @@
 
 from jdaviz import __version__
 from jdaviz.components.toolbar_nested import NestedJupyterToolbar
+from jdaviz.core.custom_traitlets import FloatHandleEmpty
 from jdaviz.core.events import (AddDataMessage, RemoveDataMessage,
                                 ViewerAddedMessage, ViewerRemovedMessage,
                                 ViewerRenamedMessage, SnackbarMessage,
                                 AddDataToViewerMessage)
+from jdaviz.core.marks import (LineAnalysisContinuum,
+                               LineAnalysisContinuumCenter,
+                               LineAnalysisContinuumLeft,
+                               LineAnalysisContinuumRight,
+                               ShadowLine)
 from jdaviz.core.region_translators import _get_region_from_spatial_subset
 from jdaviz.core.user_api import UserApiWrapper, PluginUserApi
 from jdaviz.style_registry import PopoutStyleWrapper
@@ -52,7 +60,7 @@
            'SelectPluginComponent', 'UnitSelectPluginComponent', 'EditableSelectPluginComponent',
            'PluginSubcomponent',
            'SubsetSelect', 'SpatialSubsetSelectMixin', 'SpectralSubsetSelectMixin',
-           'DatasetSpectralSubsetValidMixin',
+           'DatasetSpectralSubsetValidMixin', 'SpectralContinuumMixin',
            'ViewerSelect', 'ViewerSelectMixin',
            'LayerSelect', 'LayerSelectMixin',
            'NonFiniteUncertaintyMismatchMixin',
@@ -1979,6 +1987,214 @@ def _check_non_finite_uncertainty_mismatch(self, event={}):
         self.non_finite_uncertainty_mismatch = bool(mismatch)
 
 
+class SpectralContinuumMixin(VuetifyTemplate, HubListener):
+    """
+    Plugin select to choose options for a linear spectral continuum.
+    """
+    continuum_subset_items = List().tag(sync=True)
+    continuum_subset_selected = Unicode().tag(sync=True)
+
+    continuum_width = FloatHandleEmpty(3).tag(sync=True)
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.continuum = SubsetSelect(self,
+                                      'continuum_subset_items',
+                                      'continuum_subset_selected',
+                                      manual_options=['None', 'Surrounding'],
+                                      default_mode='first',
+                                      filters=['is_spectral'])
+
+    def _continuum_remove_none_option(self):
+        self.continuum.items = [item for item in self.continuum.items
+                                if item['label'] != 'None']
+        self.continuum._apply_default_selection()
+
+    @property
+    def continuum_marks(self):
+        marks = {}
+        viewer = self.app.get_viewer(self._default_spectrum_viewer_reference_name)
+        if viewer is None:
+            return {}
+        for mark in viewer.figure.marks:
+            if isinstance(mark, LineAnalysisContinuum):
+                # NOTE: we don't use isinstance anymore because of nested inheritance
+                if mark.__class__.__name__ == 'LineAnalysisContinuumLeft':
+                    marks['left'] = mark
+                elif mark.__class__.__name__ == 'LineAnalysisContinuumCenter':
+                    marks['center'] = mark
+                elif mark.__class__.__name__ == 'LineAnalysisContinuumRight':
+                    marks['right'] = mark
+
+        if not len(marks):
+            if not viewer.state.reference_data:
+                # we don't have data yet for scales, defer initializing
+                return {}
+            # then haven't been initialized yet, so initialize with empty
+            # marks that will be populated once the first analysis is done.
+            marks = {'left': LineAnalysisContinuumLeft(viewer, visible=self.is_active),
+                     'center': LineAnalysisContinuumCenter(viewer, visible=self.is_active),
+                     'right': LineAnalysisContinuumRight(viewer, visible=self.is_active)}
+            shadows = [ShadowLine(mark, shadow_width=2) for mark in marks.values()]
+            # NOTE: += won't trigger the figure to notice new marks
+            viewer.figure.marks = viewer.figure.marks + shadows + list(marks.values())
+
+        return marks
+
+    def _update_continuum_marks(self, mark_x={}, mark_y={}):
+        for pos, mark in self.continuum_marks.items():
+            mark.update_xy(mark_x.get(pos, []), mark_y.get(pos, []))
+
+    def _get_continuum(self, dataset, spatial_subset, spectral_subset, update_marks=False):
+        if dataset.selected == '':
+            self._update_continuum_marks()
+            return None, None, None
+
+        if spatial_subset == 'per-pixel':
+            if self.app.config != 'cubeviz':
+                raise ValueError("per-pixel only supported for cubeviz")
+            full_spectrum = self.dataset.selected_obj
+        else:
+            full_spectrum = dataset.selected_spectrum_for_spatial_subset(spatial_subset.selected if spatial_subset is not None else None,  # noqa
+                                                                         use_display_units=True)
+
+        if full_spectrum is None or self.continuum_width == "":
+            self._update_continuum_marks()
+            return None, None, None
+
+        spectral_axis = full_spectrum.spectral_axis
+        if spectral_axis.unit == u.pix:
+            # plugin should be disabled so not get this far, but can still get here
+            # before the disabled message is set
+            self._update_continuum_marks()
+            return None, None, None
+
+        if self.continuum_subset_selected == spectral_subset.selected:
+            # already raised a validation error in the UI
+            self._update_continuum_marks()
+            return None, None, None
+
+        if spectral_subset.selected == "Entire Spectrum":
+            spectrum = full_spectrum
+        else:
+            sr = self.app.get_subsets(spectral_subset.selected,
+                                      simplify_spectral=True,
+                                      use_display_units=True)
+            spectrum = extract_region(full_spectrum, sr, return_single_spectrum=True)
+            sr_lower = np.nanmin(spectrum.spectral_axis[spectrum.spectral_axis.value >= sr.lower.value])  # noqa
+            sr_upper = np.nanmax(spectrum.spectral_axis[spectrum.spectral_axis.value <= sr.upper.value])  # noqa
+
+        if self.continuum_subset_selected == 'None':
+            self._update_continuum_marks()
+            return spectrum, np.zeros_like(spectrum.flux.value), spectrum
+
+        # compute continuum
+        if self.continuum_subset_selected == "Surrounding" and spectral_subset.selected == "Entire Spectrum": # noqa
+            # we know we'll just use the endpoints, so let's be efficient and not even
+            # try extracting from the region
+            continuum_mask = np.array([0, len(spectral_axis)-1])
+            if update_marks:
+                mark_x = {'left': np.array([]),
+                          'center': np.array([min(spectral_axis.value),
+                                              max(spectral_axis.value)]),
+                          'right': np.array([])}
+
+        elif self.continuum_subset_selected == "Surrounding":
+            # self.spectral_subset_selected != "Entire Spectrum"
+            if self.continuum_width > 10 or self.continuum_width < 1:
+                # DEV NOTE: if changing the limits, make sure to also update the form validation
+                # rules in line_analysis.vue
+                self._update_continuum_marks()
+                return None, None, None
+
+            spectral_region_width = sr_upper - sr_lower
+            # convert width from total relative width, to width per "side"
+            width = (self.continuum_width - 1) / 2
+            left, = np.where((spectral_axis < sr_lower) &
+                             (spectral_axis > sr_lower - spectral_region_width*width))
+            if not len(left):
+                # then no points matching the width are available outside the line region,
+                # so we'll default to the left-most point of the line region.
+                left, = np.where(spectral_axis == min(spectrum.spectral_axis))
+
+            right, = np.where((spectral_axis > sr_upper) &
+                              (spectral_axis < sr_upper + spectral_region_width*width))
+            if not len(right):
+                # then no points matching the width are available outside the line region,
+                # so we'll default to the right-most point of the line region.
+                right, = np.where(spectral_axis == max(spectrum.spectral_axis))
+
+            continuum_mask = np.concatenate((left, right))
+            if update_marks:
+                mark_x = {'left': np.array([min(spectral_axis.value[continuum_mask]),
+                                            sr_lower.value]),
+                          'center': np.array([sr_lower.value,
+                                              sr_upper.value]),
+                          'right': np.array([sr_upper.value,
+                                             max(spectral_axis.value[continuum_mask])])}
+
+        else:
+            # we'll access the mask of the continuum and then apply that to the spectrum.  For a
+            # spatially-collapsed spectrum in cubeviz, this will access the mask from the full
+            # cube, but still apply that to the spatially-collapsed spectrum.
+            continuum_mask = ~self._specviz_helper.get_data(
+                dataset.selected,
+                spectral_subset=self.continuum_subset_selected,
+                use_display_units=False).mask
+            spectral_axis_nanmasked = spectral_axis.value.copy()
+            spectral_axis_nanmasked[~continuum_mask] = np.nan
+            if not update_marks:
+                pass
+            elif spectral_subset.selected == "Entire Spectrum":
+                mark_x = {'left': spectral_axis_nanmasked,
+                          'center': spectral_axis.value,
+                          'right': []}
+            else:
+                mark_x = {'left': spectral_axis_nanmasked[spectral_axis.value < sr_lower.value],
+                          'right': spectral_axis_nanmasked[spectral_axis.value > sr_upper.value]}
+                # Center should extend (at least) across the line region between the full
+                # range defined by the continuum subset(s).
+                # OK for mark_x to be all NaNs.
+                with warnings.catch_warnings():
+                    warnings.simplefilter('ignore', category=RuntimeWarning)
+                    mark_x_min = np.nanmin(mark_x['left'])
+                    mark_x_max = np.nanmax(mark_x['right'])
+                left_min = np.nanmin([mark_x_min, sr_lower.value])
+                right_max = np.nanmax([mark_x_max, sr_upper.value])
+                mark_x['center'] = np.array([left_min, right_max])
+
+        continuum_x = spectral_axis[continuum_mask].value
+        min_x = min(spectral_axis.value)
+        if spatial_subset == 'per-pixel':
+            # full_spectrum.flux is a cube, so we want to act on all spaxels independently
+            continuum_y = full_spectrum.flux[:, :, continuum_mask].value
+
+            def fit_continuum(continuum_y_spaxel):
+                return np.polyfit(continuum_x-min_x, continuum_y_spaxel, deg=1)
+
+            # compute the linear fit for each spaxel independently, along the spectral axis
+            slopes_intercepts = np.apply_along_axis(fit_continuum, 2, continuum_y)
+            slopes = slopes_intercepts[:, :, 0]
+            intercepts = slopes_intercepts[:, :, 1]
+
+            # spectrum.spectral_axis is an array, but we need our continuum to have the same
+            # shape as the fluxes in the cube, so let's just duplicate to the correct shape
+            spectral_axis_cube = np.zeros(spectrum.flux.shape)
+            spectral_axis_cube[:, :] = spectrum.spectral_axis.value
+
+            continuum = slopes[:, :, np.newaxis] * (spectral_axis_cube-min_x) + intercepts[:, :, np.newaxis]  # noqa
+        else:
+            continuum_y = full_spectrum.flux[continuum_mask].value
+            slope, intercept = np.polyfit(continuum_x-min_x, continuum_y, deg=1)
+            continuum = slope * (spectrum.spectral_axis.value-min_x) + intercept
+
+        if update_marks:
+            mark_y = {k: slope * (v-min_x) + intercept for k, v in mark_x.items()}
+            self._update_continuum_marks(mark_x, mark_y)
+
+        return spectrum, continuum, spectrum - continuum
+
+
 class ViewerSelect(SelectPluginComponent):
     """
     Plugin select for viewers, with support for single or multi-selection.