Using chords to find centre of a circle - getting correct coordinates in one file but wrong coordinates for others

I have been attempting to use parallel chords to find the coordinates of the centre of a circle (the sun's disk, in this case).

I have FITS files with data that I have been manipulating, and have written the below code:

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

To 'make' parallel chords, I have decided to put vertical and horizontal lines across the image at intervals of ~200 pixels.

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

This function returns the points where the line crosses the circle (x1, x2, y1, y2) and the length of each chord (xlen, ylen) for each value in the range y_values, x_values.

I find the radius using this function:

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

Then, using Pythagoras, it is possible to determine the distance to the centre of the circle using this function:

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

I have decided it would be best to take an average of the coordinates to get the central coordinates I require. This is done here:

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

This returns the coordinates I plot.

This appears to work for my template image, but as soon as I apply the same functions to a different file the centre is off.

The only values that are changed is the input data from the FITS file - I don't know what could be causing this issue, as I thought what I had written solves the issue faced by part of the image cutting off.

Can anyone help point out what's going wrong?

Thanks

I think you have over-complicated the maths and the code a little bit 😉

If you find the first and last white pixel in any row of the image, you can average them and that will give you the x-coordinate of the image centre - as long as the horizontal chord you have drawn is complete, i.e. it doesn't touch the edges of the image.

Likewise, if you find the top-most and bottom-most white pixel in any column of the image, you can average them and that will give you the y-coordinate of the image centre - with same proviso as above.

The proviso about being a complete chord is expressed in this line in the code:

if L[x]>0 and R[x]<w

I also simplified and hopefully speeded up the code using Numpy to do the processing. You should consider trying to avoid iterating over pixels with for loops - they are slow and error-prone.

#!/usr/bin/env python3

import numpy as np
import cv2 as cv

# Load the image, greyscale and threshold at 128 to Boolean mask of sun's pixels
BGR  = cv.imread('image.png')
grey = cv.cvtColor(BGR, cv.COLOR_BGR2GRAY)
mask = grey > 128
h, w = grey.shape

# Get left-most and right-most True pixel in each row
# Get top-most and bottom-most True pixel in each column
L = np.argmax(mask, axis=1)
R = w - np.argmax(mask[..., ::-1], axis=1)
T = np.argmax(mask, axis=0)
B = h - np.argmax(mask[::-1, ...], axis=0)

# Get midpoint between left and right only for rows that are complete chords
xcentres = [ (L[x]+R[x])/2 for x in range(h) if L[x]>0 and R[x]<w ]
x = int(np.mean(xcentres))
# Get midpoint between top and bottom only for columns that are complete chords
ycentres = [ (T[x]+B[x])/2 for x in range(h) if T[x]>0 and B[x]<h ]
y = int(np.mean(ycentres))

print(f'{x=}, {y=}')

# Draw centre in red
BGR[y-3:y+3, x-3:x+3] = [0,0,255]
cv.imwrite('result.png', BGR)