Calculating day length based on angular displacement

day_length.py

from sklearn.cluster import KMeans
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt


time_diff_mins = 1
shape = [3468, 4624]
cutoff = 95
box = (0, shape[1]*0.1, shape[0]*0.5, shape[1]*0.4)
file_names = ["1.jpg", "2.jpg"]
star_count = 27

img_1, img_2 = map(
        lambda file_name: np.array(Image.open(file_name).crop(box).convert("L")),
        file_names
    )

# Maximising the contrast
img_1_1D_array = np.reshape(img_1, (img_1.shape[0]*img_1.shape[1]))
img_1_1D_array[img_1_1D_array <= cutoff] = 0
img_1_1D_array[img_1_1D_array > cutoff] = 1
img_1 = np.reshape(img_1_1D_array, img_1.shape)

img_2_1D_array = np.reshape(img_2, (img_2.shape[0]*img_2.shape[1]))
img_2_1D_array[img_2_1D_array <= cutoff] = 0
img_2_1D_array[img_2_1D_array > cutoff] = 1
img_2 = np.reshape(img_2_1D_array, img_2.shape)

# seeing images
img_1_img = Image.fromarray(img_1*255)
img_1_img.save("1_contrast_100.png")

img_2_img = Image.fromarray(img_2*255)
img_2_img.save("2_contrast_100.png")

##############
# star centers
##############
# creating the feature vectors
x_1, x_2 = np.nonzero(img_1)
X_1 = np.array([x_1, x_2]).T
# clustering
kmeans_1 = KMeans(star_count, init="random").fit(X_1)
cluster_centers_1 = kmeans_1.cluster_centers_

# creating the feature vectors
x_1, x_2 = np.nonzero(img_2)
X_2 = np.array([x_1, x_2]).T
# clustering
kmeans_2 = KMeans(star_count, init="random").fit(X_2)
cluster_centers_2 = kmeans_2.cluster_centers_

# seeing them
y_1, x_1 = cluster_centers_1.T
y_2, x_2 = cluster_centers_2.T
plt.clf()
fig, ax = plt.subplots(1, 1)
ax.scatter(x_1, y_1, label="1", alpha=0.5)
ax.scatter(x_2, y_2, label="2", alpha=0.5)
ax.invert_yaxis()
plt.legend()
plt.grid()
plt.savefig("cluster_centers.png")

# pairing them
def distance(a, b):
    c = a - b
    return np.sqrt(c[0]**2 + c[1]**2)

pairs = []
for center_1 in cluster_centers_1:
    min_dist = 100
    for center_2 in cluster_centers_2:
        if distance(center_1, center_2) < min_dist:
            pair = center_2
            min_dist = distance(center_1, pair)
    pairs.append([center_1, pair, distance(center_1, pair)])

# filtering outliers
centers_1, centers_2, distances = np.array(list(filter(lambda x: x[2] < 30, pairs))).T

# average displacement vector
diff_of_centers = centers_1-centers_2
adv = sum(diff_of_centers)/len(diff_of_centers)

# results
pixel_displacement = np.sqrt(adv[0]**2 + adv[1]**2)
angular_displacement = pixel_displacement/416.535  # by comparing distance in orion's legs
angular_displacement/=(np.cos(np.deg2rad(20.25)))
minutes = 360/(angular_displacement/time_diff_mins)
hours = minutes/60
print(f"minutes: {minutes}")
print(f"hours: {hours}")
print(f"days: {hours/24}")