Add Full Width at Half Maximum analyzer (#109)

Author: vladtcvs

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]:

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

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"),
 }

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/)"),
+}