星跟踪算法中星点像素坐标平均方法及中心提取咨询

针对二值化后单星对应多像素的问题，常用以下几种方法提取星点中心：

1. 亮度加权质心法 #

这是星跟踪中最常用的方法，恒星亮度分布近似高斯，加权质心能精准反映真实星点中心：

• 提取所有亮像素的坐标和原始亮度值
• 以亮度为权重，计算坐标的加权平均值作为中心
• 若需整数像素中心，对结果四舍五入即可

2. 连通区域标记+区域中心法 #

先区分二值图像中的独立星点区域，再计算每个区域的中心：

• 用连通区域标记算法拆分不同星点
• 几何中心：取区域最小/最大坐标的平均值，计算简单但精度略低
• 区域加权质心：基于区域内亮度分布计算，精度接近全局加权质心

3. 高斯拟合法 #

对星点亮度分布做二维高斯函数拟合，峰值位置即为中心，精度最高但计算复杂度稍高：

• 截取单个星点的局部区域
• 定义二维高斯函数，通过曲线拟合求解参数
• 适合高精度场景，但星点过小或噪声大时可能拟合失败

优化后的代码示例 #

替换原循环二值化逻辑，用numpy提升效率并实现上述方法：

import pandas as pd 
import math
from PIL import Image, ImageDraw
from pillow_heif import register_heif_opener
import numpy as np
from scipy.ndimage import label, center_of_mass
from scipy.optimize import curve_fit

# 初始化与图像读取
register_heif_opener()
img = Image.open("src/sample_star.heic")
img_np = np.array(img)
data = pd.read_csv("src/data.csv")

# 快速二值化（替代原循环）
threshold = 200
brightness = np.mean(img_np, axis=2)  # 计算每个像素的亮度
binary = (brightness >= threshold).astype(np.uint8)  # 生成二值数组

# ------------------- 方法1：亮度加权质心 -------------------
y_coords, x_coords = np.where(binary == 1)
intensities = brightness[y_coords, x_coords]

if len(x_coords) > 0:
    # 计算加权质心
    centroid_x = np.sum(x_coords * intensities) / np.sum(intensities)
    centroid_y = np.sum(y_coords * intensities) / np.sum(intensities)
    # 转为整数像素中心
    centroid_x_int = round(centroid_x)
    centroid_y_int = round(centroid_y)
    print(f"全局加权质心中心：({centroid_x_int}, {centroid_y_int})")

# ------------------- 方法2：连通区域中心 -------------------
labeled_regions, star_count = label(binary)
print(f"共识别到{star_count}个星点区域")

for idx in range(1, star_count + 1):
    region_mask = labeled_regions == idx
    y, x = np.where(region_mask)
    
    # 几何中心
    geo_center_x = (x.min() + x.max()) / 2
    geo_center_y = (y.min() + y.max()) / 2
    # 区域加权质心
    weighted_center = center_of_mass(brightness, labeled_regions, idx)
    
    print(f"\n星点{idx}几何中心：({round(geo_center_x)}, {round(geo_center_y)})")
    print(f"星点{idx}区域加权质心：({round(weighted_center[1])}, {round(weighted_center[0])})")

# ------------------- 方法3：高斯拟合 -------------------
if star_count > 0:
    # 取第一个星点区域进行拟合
    first_region_mask = labeled_regions == 1
    y, x = np.where(first_region_mask)
    x_min, x_max = x.min(), x.max()
    y_min, y_max = y.min(), y.max()
    star_region = brightness[y_min:y_max+1, x_min:x_max+1]
    
    # 定义二维高斯函数
    def gaussian_2d(x, y, x0, y0, sigma_x, sigma_y, A, offset):
        return A * np.exp(-((x-x0)**2/(2*sigma_x**2) + (y-y0)**2/(2*sigma_y**2))) + offset
    
    # 拟合初始参数
    initial_guess = [(x_min+x_max)/2, (y_min+y_max)/2, 1, 1, star_region.max(), star_region.min()]
    xx, yy = np.meshgrid(np.arange(x_min, x_max+1), np.arange(y_min, y_max+1))
    
    try:
        popt, _ = curve_fit(gaussian_2d, (xx.flatten(), yy.flatten()), star_region.flatten(), p0=initial_guess)
        fit_x0, fit_y0 = popt[0], popt[1]
        print(f"\n第一个星点高斯拟合中心：({round(fit_x0)}, {round(fit_y0)})")
    except RuntimeError:
        print("\n高斯拟合失败：星点区域过小或噪声干扰")

# 可选：在图像上标记中心
img_marked = img.copy()
draw = ImageDraw.Draw(img_marked)
draw.ellipse((centroid_x_int-2, centroid_y_int-2, centroid_x_int+2, centroid_y_int+2), fill=(255,0,0), width=1)
img_marked.show()

内容的提问来源于stack exchange，提问作者theunnoticable