From 2954d20dec19224abc5fd88ace7ef41ded29aeac Mon Sep 17 00:00:00 2001 From: Chris Meyer Date: Fri, 5 Nov 2021 08:31:23 -0700 Subject: [PATCH] WIP: Add ability to do a global subtract (fit based on picked spectrum). --- nion/eels_analysis/BackgroundModel.py | 67 +++++++++++-- .../BackgroundSubtraction.py | 93 +++++++++++++++++++ .../nion_eels_analysis/__init__.py | 1 + 3 files changed, 153 insertions(+), 8 deletions(-) diff --git a/nion/eels_analysis/BackgroundModel.py b/nion/eels_analysis/BackgroundModel.py index 0f5c8f9..671391d 100644 --- a/nion/eels_analysis/BackgroundModel.py +++ b/nion/eels_analysis/BackgroundModel.py @@ -49,7 +49,7 @@ 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, @@ -57,15 +57,18 @@ def subtract_background(self, *, spectrum_xdata: DataAndMetadata.DataAndMetadata # 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"}) diff --git a/nionswift_plugin/nion_eels_analysis/BackgroundSubtraction.py b/nionswift_plugin/nion_eels_analysis/BackgroundSubtraction.py index fee04be..e96d471 100644 --- a/nionswift_plugin/nion_eels_analysis/BackgroundSubtraction.py +++ b/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, {}) diff --git a/nionswift_plugin/nion_eels_analysis/__init__.py b/nionswift_plugin/nion_eels_analysis/__init__.py index 9fb22d4..d1d0052 100755 --- a/nionswift_plugin/nion_eels_analysis/__init__.py +++ b/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))