WIP: Add ability to do a global subtract (fit based on picked spectrum).

nion/eels_analysis/BackgroundModel.py

@@ -49,23 +49,26 @@ def fit_background(self, *, spectrum_xdata: DataAndMetadata.DataAndMetadata,
                        fit_intervals: typing.Sequence[BackgroundInterval],
                        background_interval: BackgroundInterval, **kwargs) -> typing.Dict:
         return {
-            "background_model": self.__fit_background(spectrum_xdata, fit_intervals, background_interval),
+            "background_model": self.__fit_background(spectrum_xdata, None, fit_intervals, background_interval),
         }
 
     def subtract_background(self, *, spectrum_xdata: DataAndMetadata.DataAndMetadata,
                             fit_intervals: typing.Sequence[BackgroundInterval], **kwargs) -> typing.Dict:
         # set up initial values
         fit_minimum = min([fit_interval[0] for fit_interval in fit_intervals])
         signal_interval = fit_minimum, 1.0
-        subtracted_xdata = Core.calibrated_subtract_spectrum(spectrum_xdata, self.__fit_background(spectrum_xdata, fit_intervals, signal_interval))
+        subtracted_xdata = Core.calibrated_subtract_spectrum(spectrum_xdata, self.__fit_background(spectrum_xdata, None, fit_intervals, signal_interval))
         assert subtracted_xdata
         return {"subtracted": subtracted_xdata}
 
     def integrate_signal(self, *, spectrum_xdata: DataAndMetadata.DataAndMetadata,
+                         eels_spectrum_xdata: typing.Optional[DataAndMetadata.DataAndMetadata] = None,
                          fit_intervals: typing.Sequence[BackgroundInterval],
                          signal_interval: BackgroundInterval, **kwargs) -> typing.Dict:
         # set up initial values
-        subtracted_xdata = Core.calibrated_subtract_spectrum(spectrum_xdata, self.__fit_background(spectrum_xdata, fit_intervals, signal_interval))
+        subtracted_xdata = Core.calibrated_subtract_spectrum(spectrum_xdata,
+                                                             self.__fit_background(spectrum_xdata, eels_spectrum_xdata,
+                                                                                   fit_intervals, signal_interval))
         assert subtracted_xdata
         subtracted_data = subtracted_xdata.data
         assert subtracted_data is not None
@@ -81,6 +84,7 @@ def integrate_signal(self, *, spectrum_xdata: DataAndMetadata.DataAndMetadata,
             }
 
     def __fit_background(self, spectrum_xdata: DataAndMetadata.DataAndMetadata,
+                         eels_spectrum_xdata: typing.Optional[DataAndMetadata.DataAndMetadata],
                          fit_intervals: typing.Sequence[BackgroundInterval],
                          background_interval: BackgroundInterval) -> DataAndMetadata.DataAndMetadata:
         # fit polynomial to the data
@@ -93,6 +97,16 @@ def __fit_background(self, spectrum_xdata: DataAndMetadata.DataAndMetadata,
                  fit_intervals])
         else:
             ys = get_calibrated_interval_slice(spectrum_xdata, fit_intervals[0]).data
+        es: typing.Optional[numpy.ndarray]
+        if eels_spectrum_xdata:
+            if len(fit_intervals) > 1:
+                es = numpy.concatenate(
+                    [get_calibrated_interval_slice(eels_spectrum_xdata, fit_interval).data for fit_interval in
+                     fit_intervals])
+            else:
+                es = get_calibrated_interval_slice(eels_spectrum_xdata, fit_intervals[0]).data
+        else:
+            es = None
         # generate background model data from the series
         background_interval_start_pixel = round(spectrum_xdata.data_shape[-1] * background_interval[0])
         background_interval_end_pixel = round(spectrum_xdata.data_shape[-1] * background_interval[1])
@@ -106,7 +120,7 @@ def __fit_background(self, spectrum_xdata: DataAndMetadata.DataAndMetadata,
         if spectrum_xdata.is_navigable:
             calibrations = list(copy.deepcopy(spectrum_xdata.navigation_dimensional_calibrations)) + [calibration]
             yss = numpy.reshape(ys, (numpy.product(ys.shape[:-1]),) + (ys.shape[-1],))
-            fit_data = self._perform_fits(xs, yss, fs)
+            fit_data = self._perform_fits(xs, yss, fs, es)
             data_descriptor = DataAndMetadata.DataDescriptor(False, spectrum_xdata.navigation_dimension_count,
                                                              spectrum_xdata.datum_dimension_count)
             background_xdata = DataAndMetadata.new_data_and_metadata(numpy.reshape(fit_data, ys.shape[:-1] + (n,)),
@@ -119,10 +133,11 @@ def __fit_background(self, spectrum_xdata: DataAndMetadata.DataAndMetadata,
                                                                      intensity_calibration=spectrum_xdata.intensity_calibration)
         return background_xdata
 
-    def _perform_fits(self, xs: numpy.ndarray, yss: numpy.ndarray, fs: numpy.ndarray) -> numpy.ndarray:
+    def _perform_fits(self, xs: numpy.ndarray, yss: numpy.ndarray, fs: numpy.ndarray, es: typing.Optional[numpy.ndarray]) -> numpy.ndarray:
         # xs will be a set of x-values with shape (L) representing the energies at which to fit
         # ys will be an array of y-values with shape (m,L)
         # fs will be an array of x-values with shape (n) representing energies at which to generate fitted data
+        # es will be an optional array of y-values with shape (L) representing the spectrum on which the background fit may be based
         # return an ndarray of the fit with shape (m,n)
         # implement at least one of _perform_fits and _perform_fit
         yss_shape = yss.shape[:-1]
@@ -137,7 +152,7 @@ def _perform_fit(self, xs: numpy.ndarray, ys: numpy.ndarray, fs: numpy.ndarray)
         # fs will be an array of x-values with shape (n) representing energies at which to generate fitted data
         # return an ndarray of the fit with shape (n)
         # implement at least one of _perform_fits and _perform_fit
-        return numpy.reshape(self._perform_fits(xs, numpy.reshape(ys, (1,) + ys.shape), fs), fs.shape)
+        return numpy.reshape(self._perform_fits(xs, numpy.reshape(ys, (1,) + ys.shape), fs, None), fs.shape)
 
 
 class PolynomialBackgroundModel(AbstractBackgroundModel):
@@ -148,7 +163,7 @@ def __init__(self, background_model_id: str, deg: int, transform=None, untransfo
         self.transform = transform
         self.untransform = untransform
 
-    def _perform_fits(self, xs: numpy.ndarray, yss: numpy.ndarray, fs: numpy.ndarray) -> numpy.ndarray:
+    def _perform_fits(self, xs: numpy.ndarray, yss: numpy.ndarray, fs: numpy.ndarray, es: typing.Optional[numpy.ndarray]) -> numpy.ndarray:
         transform_data = self.transform or (lambda x: x)
         untransform_data = self.untransform or (lambda x: x)
         coefficients = numpy.polynomial.polynomial.polyfit(transform_data(xs), transform_data(yss.transpose()), self.deg)
@@ -163,6 +178,40 @@ def __unused_perform_fit(self, xs: numpy.ndarray, ys: numpy.ndarray, fs: numpy.n
         return untransform_data(series(fs))
 
 
+def power_law(x: numpy.ndarray, a: float, b: float) -> numpy.ndarray:
+    return a * (x) ** -b
+
+
+class SingleParamPowerLawBackgroundModel(AbstractBackgroundModel):
+
+    def __init__(self, background_model_id: str, title: str = None) -> None:
+        super().__init__(background_model_id, title)
+
+    def _perform_fit(self, xs: numpy.ndarray, ys: numpy.ndarray, fs: numpy.ndarray) -> numpy.ndarray:
+        yss = numpy.reshape(ys, (1,) + ys.shape)
+        coefficients = numpy.polynomial.polynomial.polyfit(xs, yss.transpose(), 1)
+        fit = numpy.polynomial.polynomial.polyval(fs, coefficients)
+        return numpy.reshape(numpy.where(numpy.isfinite(fit), fit, 0), fs.shape)
+
+        # params = scipy.optimize.curve_fit(power_law, xs, ys)
+        # global_model = numpy.array([power_law(i, *params[0]) for i in fs])
+        # global_model_fit = numpy.array([power_law(i, *params[0]) for i in xs])
+        # global_model_norm_factor = numpy.sqrt(numpy.sum(global_model_fit ** 2))
+        # global_model_fit_normed = global_model_fit / global_model_norm_factor
+        # global_model_normed = global_model / global_model_norm_factor
+        # amplitude = numpy.sum(global_model_fit_normed * ys)  # this is the "fit"
+        # return amplitude * global_model_normed
+
+    def _perform_fits(self, xs: numpy.ndarray, yss: numpy.ndarray, fs: numpy.ndarray, es: typing.Optional[numpy.ndarray]) -> numpy.ndarray:
+        # use the es array if available or else just pixel by pixel
+        zs = numpy.reshape(es, (1,) + es.shape).transpose() if es is not None else yss.transpose()
+        coefficients1 = numpy.polynomial.polynomial.polyfit(xs, zs, 1)
+        coefficients = numpy.empty((coefficients1.shape[0], yss.shape[0]), dtype=coefficients1.dtype)
+        coefficients[:] = coefficients1
+        fit = numpy.polynomial.polynomial.polyval(fs, coefficients)
+        return numpy.where(numpy.isfinite(fit), fit, 0)
+
+
 class TwoAreaBackgroundModel(AbstractBackgroundModel):
     # Fit power law or exponential background model using the two-area method described in Egerton chapter 4.
     # This approximation is slightly faster than the polynomial fit for mapping large SI, and may perform better for high-noise spectra.
@@ -177,7 +226,7 @@ def __init__(self,
         self.model_func = model_func
         self.params_func = params_func
 
-    def _perform_fits(self, xs: numpy.ndarray, yss: numpy.ndarray, fs: numpy.ndarray) -> numpy.ndarray:
+    def _perform_fits(self, xs: numpy.ndarray, yss: numpy.ndarray, fs: numpy.ndarray, es: typing.Optional[numpy.ndarray]) -> numpy.ndarray:
         half_interval = len(xs) // 2
         x_interval_1 = xs[:half_interval]
         x_interval_2 = xs[half_interval:2 * half_interval]
@@ -254,3 +303,5 @@ def exponential_func(x: numpy.ndarray, A: numpy.ndarray, tau: numpy.ndarray) ->
 
 Registry.register_component(TwoAreaBackgroundModel("exponential_two_area_background_model", params_func=exponential_params, model_func=exponential_func,
                                                    title=_("Exponential Two Area")), {"background-model"})
+
+Registry.register_component(SingleParamPowerLawBackgroundModel("power_law_global_background_model", title=_("Power Law Global")), {"background-model"})

nionswift_plugin/nion_eels_analysis/BackgroundSubtraction.py

@@ -170,6 +170,61 @@ def commit(self):
         self.computation.set_referenced_xdata("map", self.__mapped_xdata)
 
 
+class EELSMapBackgroundSubtractedSignal1:
+    label = _("EELS Map Background Subtracted Signal 1")
+    inputs = {
+        "spectrum_image_data_item": {"label": _("EELS Image")},
+        "eels_spectrum_data_item": {"label": _("EELS Spectrum")},
+        "background_model": {"label": _("Background Model"), "entity_id": "background_model"},
+        "fit_interval_graphics": {"label": _("Fit")},
+        "signal_interval_graphic": {"label": _("Signal")},
+        }
+    outputs = {
+        "map": {"label": _("EELS Signal")},
+    }
+
+    def __init__(self, computation, **kwargs):
+        self.computation = computation
+
+    def execute(self, spectrum_image_data_item: Facade.DataItem, eels_spectrum_data_item: Facade.DataItem, background_model, fit_interval_graphics, signal_interval_graphic):
+        try:
+            assert spectrum_image_data_item.xdata.is_datum_1d
+            assert spectrum_image_data_item.xdata.is_navigable
+            assert spectrum_image_data_item.xdata.datum_dimensional_calibrations[0].units == "eV"
+            spectrum_image_xdata = spectrum_image_data_item.xdata
+            spectrum_xdata = eels_spectrum_data_item.xdata
+            assert spectrum_xdata.is_datum_1d
+            assert spectrum_xdata.datum_dimensional_calibrations[0].units == "eV"
+            eels_spectrum_xdata = spectrum_xdata
+            # fit_interval_graphics.interval returns normalized coordinates. create calibrated intervals.
+            fit_intervals: typing.List[BackgroundModel.BackgroundInterval] = list()
+            for fit_interval_graphic in fit_interval_graphics:
+                fit_intervals.append(fit_interval_graphic.interval)
+            signal_interval = signal_interval_graphic.interval
+            mapped_xdata = None
+            if background_model._data_structure.entity:
+                entity_id = background_model._data_structure.entity.entity_type.entity_id
+                for component in Registry.get_components_by_type("background-model"):
+                    if entity_id == component.background_model_id:
+                        # import time
+                        # t0 = time.perf_counter()
+                        integrate_result = component.integrate_signal(spectrum_xdata=spectrum_image_xdata, eels_spectrum_xdata=eels_spectrum_xdata, fit_intervals=fit_intervals, signal_interval=signal_interval)
+                        # t1 = time.perf_counter()
+                        # print(f"{component.background_model_id} {((t1 - t0) * 1000)}ms")
+                        mapped_xdata = integrate_result["integrated"]
+            if mapped_xdata is None:
+                mapped_xdata = DataAndMetadata.new_data_and_metadata(numpy.zeros(spectrum_image_xdata.navigation_dimension_shape), dimensional_calibrations=spectrum_image_xdata.navigation_dimensional_calibrations)
+            self.__mapped_xdata = mapped_xdata
+        except Exception as e:
+            import traceback
+            print(traceback.format_exc())
+            print(e)
+            raise
+
+    def commit(self):
+        self.computation.set_referenced_xdata("map", self.__mapped_xdata)
+
+
 def add_background_subtraction_computation(api: Facade.API_1, library: Facade.Library, display_item: Facade.Display, data_item: Facade.DataItem, intervals: typing.Sequence[Facade.Graphic]) -> None:
     background = api.library.create_data_item(title="{} Background".format(data_item.title))
     signal = api.library.create_data_item(title="{} Subtracted".format(data_item.title))
@@ -278,8 +333,46 @@ def use_signal_for_map(api, window):
                     break
 
 
+def use_signal_for_global_map(api, window):
+    target_display = window.target_display
+    target_graphic = target_display.selected_graphics[0] if target_display and len(target_display.selected_graphics) == 1 else None
+    target_interval = target_graphic if target_graphic and target_graphic.graphic_type == "interval-graphic" else None
+    if target_display and target_interval:
+        target_display_item_data_items = target_display._display_item.data_items
+        for computation in api.library._document_model.computations:
+            if computation.processing_id == "eels.background_subtraction3":
+                if computation.get_input("eels_spectrum_data_item") in target_display_item_data_items and computation.get_output("subtracted") in target_display_item_data_items:
+                    eels_spectrum_data_item = computation.get_input("eels_spectrum_data_item")
+                    eels_spectrum_data_item = api._new_api_object(eels_spectrum_data_item)
+                    fit_interval_graphics = computation.get_input("fit_interval_graphics")
+                    fit_interval_graphics = [api._new_api_object(g) for g in fit_interval_graphics]
+                    background_model = computation.get_input("background_model")
+                    background_model = api._new_api_object(background_model)
+                    source_data_items = api.library._document_model.get_source_data_items(eels_spectrum_data_item._data_item)
+                    if len(source_data_items) == 1 and source_data_items[0].xdata.is_navigable and source_data_items[0].datum_dimension_count == 1:
+                        spectrum_image = api._new_api_object(source_data_items[0])
+                        map = api.library.create_data_item_from_data(numpy.zeros(spectrum_image._data_item.xdata.navigation_dimension_shape), title="{} Map".format(spectrum_image.title))
+                        signal_interval_graphic = target_interval
+                        api.library.create_computation(
+                            "eels.global_mapping",
+                            inputs={
+                                "spectrum_image_data_item": spectrum_image,
+                                "eels_spectrum_data_item": eels_spectrum_data_item,
+                                "fit_interval_graphics": fit_interval_graphics,
+                                "signal_interval_graphic": signal_interval_graphic,
+                                "background_model": background_model,
+                            },
+                            outputs={
+                                "map": map
+                            }
+                        )
+                        window.display_data_item(map)
+                    break
+
+
 Symbolic.register_computation_type("eels.background_subtraction3", EELSFitBackground)
 Symbolic.register_computation_type("eels.mapping3", EELSMapBackgroundSubtractedSignal)
+Symbolic.register_computation_type("eels.global_mapping", EELSMapBackgroundSubtractedSignal1)
 Symbolic.register_computation_type("eels.subtract_background", EELSSubtractBackground)
 
 BackgroundModelEntity = Schema.entity("background_model", None, None, {})

nionswift_plugin/nion_eels_analysis/__init__.py

@@ -48,6 +48,7 @@ def __build_menus(self, document_window):
         eels_menu.add_menu_item(_("Subtract Background"), functools.partial(BackgroundSubtraction.subtract_background, api, window))
         eels_menu.add_separator()
         eels_menu.add_menu_item(_("Map Signal"), functools.partial(BackgroundSubtraction.use_signal_for_map, api, window))
+        eels_menu.add_menu_item(_("Map Signal Global"), functools.partial(BackgroundSubtraction.use_signal_for_global_map, api, window))
         eels_menu.add_menu_item(_("Map Thickness"), functools.partial(ThicknessMap.map_thickness, api, window))
         eels_menu.add_separator()
         eels_menu.add_menu_item(_("Align ZLP (max method)"), functools.partial(AlignZLP.align_zlp, api, window))