From 813270d6f8d0be25cd3c0f5bbad28ee12b81123c Mon Sep 17 00:00:00 2001 From: Ivan Usov Date: Tue, 19 Oct 2021 10:59:06 +0200 Subject: [PATCH] Refactor fit_event --- pyzebra/app/panel_hdf_param_study.py | 95 ++++++---------------------- pyzebra/ccl_process.py | 26 +++++++- 2 files changed, 43 insertions(+), 78 deletions(-) diff --git a/pyzebra/app/panel_hdf_param_study.py b/pyzebra/app/panel_hdf_param_study.py index 652e9ea..502a680 100644 --- a/pyzebra/app/panel_hdf_param_study.py +++ b/pyzebra/app/panel_hdf_param_study.py @@ -1,6 +1,5 @@ import base64 import io -import math import os import numpy as np @@ -38,7 +37,6 @@ from bokeh.models import ( WheelZoomTool, ) from bokeh.palettes import Cividis256, Greys256, Plasma256 # pylint: disable=E0611 -from scipy.optimize import curve_fit import pyzebra @@ -447,68 +445,6 @@ def create(): ) proj_display_min_spinner.on_change("value", proj_display_min_spinner_callback) - def fit_event(scan): - p0 = [1.0, 0.0, 1.0] - maxfev = 100000 - - # wave = scan["wave"] - # ddist = scan["ddist"] - # cell = scan["cell"] - - # gamma = scan["gamma"][0] - # omega = scan["omega"][0] - # nu = scan["nu"][0] - # chi = scan["chi"][0] - # phi = scan["phi"][0] - - scan_motor = scan["scan_motor"] - var_angle = scan[scan_motor] - - x0 = int(np.floor(det_x_range.start)) - xN = int(np.ceil(det_x_range.end)) - y0 = int(np.floor(det_y_range.start)) - yN = int(np.ceil(det_y_range.end)) - fr0 = int(np.floor(frame_range.start)) - frN = int(np.ceil(frame_range.end)) - data_roi = scan["data"][fr0:frN, y0:yN, x0:xN] - - cnts = np.sum(data_roi, axis=(1, 2)) - coeff, _ = curve_fit(gauss, range(len(cnts)), cnts, p0=p0, maxfev=maxfev) - - # m = cnts.mean() - # sd = cnts.std() - # snr_cnts = np.where(sd == 0, 0, m / sd) - - frC = fr0 + coeff[1] - var_F = var_angle[math.floor(frC)] - var_C = var_angle[math.ceil(frC)] - # frStep = frC - math.floor(frC) - var_step = var_C - var_F - # var_p = var_F + var_step * frStep - - # if scan_motor == "gamma": - # gamma = var_p - # elif scan_motor == "omega": - # omega = var_p - # elif scan_motor == "nu": - # nu = var_p - # elif scan_motor == "chi": - # chi = var_p - # elif scan_motor == "phi": - # phi = var_p - - intensity = coeff[1] * abs(coeff[2] * var_step) * math.sqrt(2) * math.sqrt(np.pi) - - projX = np.sum(data_roi, axis=(0, 1)) - coeff, _ = curve_fit(gauss, range(len(projX)), projX, p0=p0, maxfev=maxfev) - x_pos = x0 + coeff[1] - - projY = np.sum(data_roi, axis=(0, 2)) - coeff, _ = curve_fit(gauss, range(len(projY)), projY, p0=p0, maxfev=maxfev) - y_pos = y0 + coeff[1] - - scan["fit"] = {"frame": frC, "x_pos": x_pos, "y_pos": y_pos, "intensity": intensity} - metadata_table_source = ColumnDataSource(dict(geom=[""], temp=[None], mf=[None])) metadata_table = DataTable( source=metadata_table_source, @@ -556,7 +492,15 @@ def create(): def proc_all_button_callback(): for scan in zebra_data: - fit_event(scan) + pyzebra.fit_event( + scan, + int(np.floor(frame_range.start)), + int(np.ceil(frame_range.end)), + int(np.floor(det_y_range.start)), + int(np.ceil(det_y_range.end)), + int(np.floor(det_x_range.start)), + int(np.ceil(det_x_range.end)), + ) _update_table() @@ -573,7 +517,15 @@ def create(): proc_all_button.on_click(proc_all_button_callback) def proc_button_callback(): - fit_event(det_data) + pyzebra.fit_event( + det_data, + int(np.floor(frame_range.start)), + int(np.ceil(frame_range.end)), + int(np.floor(det_y_range.start)), + int(np.ceil(det_y_range.end)), + int(np.floor(det_x_range.start)), + int(np.ceil(det_x_range.end)), + ) _update_table() @@ -628,14 +580,3 @@ def create(): tab_layout = column(row(import_layout, scan_layout, plots)) return Panel(child=tab_layout, title="hdf param study") - - -def gauss(x, *p): - """Defines Gaussian function - Args: - A - amplitude, mu - position of the center, sigma - width - Returns: - Gaussian function - """ - A, mu, sigma = p - return A * np.exp(-((x - mu) ** 2) / (2.0 * sigma ** 2)) diff --git a/pyzebra/ccl_process.py b/pyzebra/ccl_process.py index 650722e..b0c3fe5 100644 --- a/pyzebra/ccl_process.py +++ b/pyzebra/ccl_process.py @@ -1,7 +1,7 @@ import os import numpy as np -from lmfit.models import GaussianModel, LinearModel, PseudoVoigtModel, VoigtModel +from lmfit.models import Gaussian2dModel, GaussianModel, LinearModel, PseudoVoigtModel, VoigtModel from scipy.integrate import simpson, trapezoid from .ccl_io import CCL_ANGLES @@ -244,3 +244,27 @@ def get_area(scan, area_method, lorentz): area_s = np.abs(area_s * corr_factor) scan["area"] = (area_v, area_s) + + +def fit_event(scan, fr_from, fr_to, y_from, y_to, x_from, x_to): + data_roi = scan["data"][fr_from:fr_to, y_from:y_to, x_from:x_to] + + model = GaussianModel() + fr = np.arange(fr_from, fr_to) + counts_per_fr = np.sum(data_roi, axis=(1, 2)) + params = model.guess(counts_per_fr, fr) + result = model.fit(counts_per_fr, x=fr, params=params) + frC = result.params["center"].value + intensity = result.params["height"].value + + model = Gaussian2dModel() + xs, ys = np.meshgrid(np.arange(x_from, x_to), np.arange(y_from, y_to)) + xs = xs.flatten() + ys = ys.flatten() + counts = np.sum(data_roi, axis=0).flatten() + params = model.guess(counts, xs, ys) + result = model.fit(counts, x=xs, y=ys, params=params) + xC = result.params["centerx"].value + yC = result.params["centery"].value + + scan["fit"] = {"frame": frC, "x_pos": xC, "y_pos": yC, "intensity": intensity}