diff --git a/src/vstarstack/library/data.py b/src/vstarstack/library/data.py index 7b225a12..674db4a2 100644 --- a/src/vstarstack/library/data.py +++ b/src/vstarstack/library/data.py @@ -126,8 +126,16 @@ def add_parameter(self, value, name : str): """Add parameter""" self.params[name] = value - def get_channel(self, channel) -> Tuple[dict, np.ndarray]: - """Get channel image""" + def get_channel(self, channel : str) -> Tuple[np.ndarray, dict]: + """ + Get channel image. + + Parameters: + channel (str) - channel name + + Returns: + tuple(image, options) - ndarray with pixel data and channel options + """ return self.channels[channel]["data"], self.channels[channel]["options"] def get_channels(self) -> List[str]: diff --git a/src/vstarstack/library/fwhm.py b/src/vstarstack/library/fwhm.py new file mode 100644 index 00000000..c581a545 --- /dev/null +++ b/src/vstarstack/library/fwhm.py @@ -0,0 +1,97 @@ +from typing import Tuple +from vstarstack.library.data import DataFrame +import numpy as np + +def get_profile(image : np.ndarray): + y = int(image.shape[0]/2) + x = int(image.shape[1]/2) + profile1 = image[y,:] + profile2 = image[:,x] + return (profile1 + profile2)/2 + +def get_star(image : np.ndarray, x : int, y : int, r : int): + w = image.shape[1] + h = image.shape[0] + x1 = x - r + x2 = x + r + 1 + y1 = y - r + y2 = y + r + 1 + if x < r or y < r: + return None + if x > w-r-1 or y > h-r-1: + return None + return image[y1:y2, x1:x2] + +def _interpolate(v1 : float, v2 : float, v : float): + if v1 > v: + return 0 + if v2 < v: + return 1 + return (v - v1)/(v2-v1) + +def get_width2(profile : np.ndarray): + size = profile.shape[0] + part = int(size/8) + left = profile[0:part] + right = profile[size-part-1:size-1] + background = np.median(list(left) + list(right)) + noback = profile - background + + maxv = np.amax(noback) + center = int(size/2) + for i in range(size): + if noback[i] > noback[center]: + center = i + x1 = center + while True: + val = noback[x1] + if val <= maxv/2: + break + x1 -= 1 + if x1 == 0: + break + + x2 = center + while True: + val = noback[x2] + if val <= maxv/2: + break + x2 += 1 + if x2 == size-1: + break + + vl1 = noback[x1] + vl2 = noback[x1+1] + + vr1 = noback[x2] + vr2 = noback[x2-1] + + d = _interpolate(vl1, vl2, maxv/2) + x1 += d + + d = _interpolate(vr1, vr2, maxv/2) + x2 -= d + + return x2 - x1 + +def find_fwhm(image : np.ndarray, x : int, y : int, radius : int) -> float: + maxv = np.amax(image) + si = get_star(image, x, y, radius*8+1) + if si is None: + return None + if np.amax(si) == maxv: + return None + profile = get_profile(si) + width = get_width2(profile) + return width + +def find_fwhm_df(image : DataFrame, x : int, y : int, radius : int) -> Tuple[dict,float]: + channels = image.get_channels() + widths = {} + for channel in channels: + layer, opts = image.get_channel(channel) + if not opts["brightness"]: + continue + widths[channel] = find_fwhm(layer, x, y, radius) + mean = np.median(list(widths.values())) + return widths, mean diff --git a/src/vstarstack/tool/analyzers/analyzers.py b/src/vstarstack/tool/analyzers/analyzers.py index 775dd493..12fd35ee 100644 --- a/src/vstarstack/tool/analyzers/analyzers.py +++ b/src/vstarstack/tool/analyzers/analyzers.py @@ -14,4 +14,5 @@ commands = { "measure-mag" : ("vstarstack.tool.analyzers.measure_mag", "measure star magnitude"), + "measure-fwhm" : ("vstarstack.tool.analyzers.measure_fwhm", "measure full width at half maximum"), } diff --git a/src/vstarstack/tool/analyzers/measure_fwhm.py b/src/vstarstack/tool/analyzers/measure_fwhm.py new file mode 100644 index 00000000..ffab3997 --- /dev/null +++ b/src/vstarstack/tool/analyzers/measure_fwhm.py @@ -0,0 +1,65 @@ +import vstarstack.library.data +import vstarstack.library.fwhm +import vstarstack.library.stars.detect +import vstarstack.tool.common +import numpy as np +import os + +def _measure_single_fwhm(project, argv): + fname = argv[0] + x = int(argv[1]) + y = int(argv[2]) + if len(argv) > 3: + r = int(argv[3]) + else: + r = 4 + df = vstarstack.library.data.DataFrame.load(fname) + channels, mean = vstarstack.library.fwhm.find_fwhm_df(df, x, y, r) + print("Mean: %f" % mean) + for channel in channels: + print("\t%s: %f" % (channel, channels[channel])) + +def _measure_mean_fwhm_file(fname : str): + df = vstarstack.library.data.DataFrame.load(fname) + vstarstack.library.stars.detect.configure_detector(thresh_coeff=1.2) + chvs = {} + for channel in df.get_channels(): + layer, opts = df.get_channel(channel) + if not opts["brightness"]: + continue + stars = vstarstack.library.stars.detect.detect_stars(layer) + values = [] + for star in stars: + x = int(star["x"]+0.5) + y = int(star["y"]+0.5) + r = int(star["radius"]) + fwhm = vstarstack.library.fwhm.find_fwhm(layer, x, y, r) + if fwhm is not None: + values.append(fwhm) + fwhm = np.median(values) + chvs[channel] = fwhm + print("Mean: %f" % np.median(list(chvs.values()))) + for channel in chvs: + print("\t%s: %f" % (channel, chvs[channel])) + +def _measure_mean_fwhm_dir(dirname : str): + files = vstarstack.tool.common.listfiles(dirname, ".zip") + for name, filename in files: + print(name) + _measure_mean_fwhm_file(filename) + +def _measure_mean_fwhm(project, argv): + path = argv[0] + if os.path.isdir(path): + _measure_mean_fwhm_dir(path) + else: + _measure_mean_fwhm_file(path) + +commands = { + "star" : (_measure_single_fwhm, + "calculate FWHM of specific star", + "image.zip x y"), + "mean" : (_measure_mean_fwhm, + "calculate mean FWHM of stars", + "(image.zip | path/)"), +}