利用平行弦法求日面圆心：单文件有效多文件失效问题排查

我通过平行弦法基于FITS文件计算太阳盘面（圆形）中心坐标，Python代码逻辑如下：

1. 读取FITS文件数据，通过阈值筛选日面区域；
2. 生成间隔200像素的横竖线获取平行弦，记录弦与圆的交点及弦长；
3. 利用x方向极值计算半径，结合勾股定理推导圆心距离；
4. 取坐标平均值得到圆心。

代码在模板FITS文件上能得到正确结果，但更换其他FITS文件时，计算出的圆心坐标出现偏移。我仅更换了输入数据，原以为代码已处理图像截断问题，却仍异常，请求帮忙排查原因。

代码实现 #

import fits
import numpy as np

# 读取FITS文件
hdul = fits.open(r'template_file.fits')
header = hdul[0].header
data = hdul[0].data
hdul.close()

# 阈值筛选日面区域
y_indices, x_indices = np.where(data >= 6e4) 
# this threshold selects the disk as the intensity is much higher than of the background

def calculate_indices(x_indices, y_indices):
    x_values = range(200, 1201, 200)
    y_values = range(200, 1001, 200)

    # Initialize dictionaries to store the results
    y_1 = {}
    y_2 = {}
    ylen = {}
    x_1 = {}
    x_2 = {}
    xlen = {}

    # Loop over the x values
    for x in x_values:
        # Get the y values corresponding to the current x value
        y_at_x = y_indices[x_indices == x]
        
        # Store the first and last y values and their difference
        y_1[x] = y_at_x[0]
        y_2[x] = y_at_x[-1]
        ylen[x] = y_2[x] - y_1[x]

    # Loop over the y values
    for y in y_values:
        # Get the x values corresponding to the current y value
        x_at_y = x_indices[y_indices == y]
        
        # Store the first and last x values and their difference
        x_1[y] = x_at_y[0]
        x_2[y] = x_at_y[-1]
        xlen[y] = x_2[y] - x_1[y]

    return y_1, y_2, ylen, x_1, x_2, xlen

def calculate_radius(x_indices):

    left = x_indices.min() #using only x_indices as some images are cut off in the y-direction
    right = x_indices.max()

    r = (right - left)/2

    return r

def calculate_results(r, ylen, xlen):
    y_results = {}
    y_vals = range(200, 1201, 200)

    for y in y_vals:
        y_results[y] = y + np.sqrt(r**2 - (ylen[y]/2)**2)

    x_results = {}
    x_vals = range(200, 1001, 200)

    for x in x_vals:
        x_results[x] = x + np.sqrt(r**2 - (xlen[x]/2)**2)
    
    return y_results, x_results

def calculate_average(x_results, y_results): 
    a_values = [200, 400, 600] #arbitrary name
    #only taking the values in the positive hemisphere hence why up to 600

    a = np.array([y_results[a] for a in a_values])
    b = np.array([x_results[a] for a in a_values])

    av_x = np.average(a)
    av_y = np.average(b)

    return av_x, av_y

问题排查方向 #

1. 固定阈值适配性差 #

代码中硬编码6e4作为日面筛选阈值，但不同FITS文件的日面亮度、背景噪声差异极大：

• 若新文件日面亮度低于该阈值，会导致筛选出的日面区域不完整；
• 若背景噪声过高，会误将背景像素纳入日面，干扰弦端点的计算。

验证方法：对新FITS文件，可视化(x_indices, y_indices)点集，检查是否完整覆盖太阳盘面且无多余背景点。

2. 半径计算逻辑存在缺陷 #

calculate_radius仅依赖x方向的极值差计算半径，默认x方向完整包含太阳圆的直径：

• 若新文件x方向存在截断，或日面在x方向并非完整直径，计算出的半径会偏差，直接导致后续勾股定理推导的圆心坐标错误。

优化思路：

• 结合所有有效弦长估算半径，比如对每个弦长ylen[x]，计算对应的半径候选值sqrt( (ylen[x]/2)^2 + dx^2 )（dx为弦到圆心的水平距离，可通过迭代优化）；
• 对多个半径候选值取平均，提升鲁棒性。

3. 弦端点提取未排序 #

np.where返回的索引是按行优先（先y后x）排列的，并非对每个x对应的y值排序。代码中直接取y_at_x[0]和y_at_x[-1]会得到错误的弦端点，导致弦长计算错误。

修复代码：对每个x对应的y值先排序：

y_at_x = np.sort(y_indices[x_indices == x])

同理处理y对应的x值：

x_at_y = np.sort(x_indices[y_indices == y])

4. 圆心平均样本量过少 #

calculate_average仅选取[200,400,600]三个位置的结果取平均，样本量不足且覆盖范围有限：

• 若这几个位置的弦因截断或噪声出现误差，会直接影响最终圆心；
• 若新文件在这些位置存在截断，计算出的结果会完全偏离真实圆心。

优化思路：

• 纳入所有有效弦的计算结果进行平均；
• 先判断弦是否完整（比如弦长是否符合半径预期），过滤截断的弦后再计算平均。

5. 坐标系混淆风险 #

需确认FITS文件的data维度是否为(y, x)，代码中y_indices, x_indices = np.where(...)是否匹配图像真实坐标；部分FITS文件包含WCS坐标信息，若代码未考虑，可能导致坐标映射错误。

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