별 사진을 찍은 경험이 누구나 있을텐데요. 밤하늘에 아름답게 빛나는 별에 대해서 자세히 알거나 구분할 수 있는 사람은 의외로 드물 것 같습니다. 저도 어릴적에는 막연하게 낭만적으로 동경했던 정도인데, 근래에 시간과 여유가 생겨 우주에 관심을 가지고 취미생활(?)로 공부를 조금씩 해 보니 알아가는 과정이 상당히 재밌네요. 여러분도 이런 즐거움을 함께 느낄 수 있으면 좋겠습니다.^^;
사실 별은 하늘에 빽빽, 촘촘히 무수히 많을 테지만, 현재 주변환경에서 인간의 눈에 잘 관찰되는 북극성 등 247개의 대표적인 별을 대상으로 촬영한 별사진에 이름표를 붙여보려고 합니다! 이미 인터넷 사이트(nova.astrometry.net)에 접속하면 알아서 별자리를 그려주기도 하고 실시간으로 촬영화면에서 플레이트 솔버 기능을 하는 훌륭한 상용프로그램이 많은 줄 알지만,,,
제 블로그의 차별화되는 장점으로 말할 것 같으면, 별을 사랑하는 아마추어 일반인 시선으로 원리를 쉽게 접근하고 흥미를 가질 수 있게 순수히 제 노력과 고민으로 만든 것이라는 점입니다. ㅎㅎㅎ
사실 첨부터 의도해서 만든 것은 아니고 이전에 별자리판을 만들면서 힌트를 얻고, 재미없는 TV드라마 볼 바엔 그 시간에 만들어 보자는 생각이었습니다. 머릿속 아이디어가 대략 맞아 들어가고 코딩하면서 구체화되고 결국 구현이 된다는 사실에 저도 놀라운 경험이었습니다.
서론이 길었는데, 이제 본론으로 가겠습니다.
먼저 아이디어 설명입니다.
1. 특정 화각의 카메라를 하늘에 대고 어떻게 찍든, 우리가 천구의 중심에서 천구에 접하는 사각형을 볼 때 (적절히 조정만 하면) 이를 통해 보게 되는 천구화면과 사진이 동일한 시점이 있다.(천구에 접하는 사각형을 구하는 것이 관건이고 거기서 다 계산하는 것이다)
2. 아스라이 먼 별들이 천구의 벽 (반지름 r=1)에 박혀있고, 별 사이의 거리는 별 사이의 각과 같다.
3. 별 3개의 세변의 길이가 정해지면 기하적인 모양과 별이름을 특정할 수 있다.
4. 고등학교의 내적, 삼각함수와 대학교 구좌표>직각좌표 변환공식 정도를 이용했다.
그럼 이제 구체적으로 파이썬 소스코드를 살펴보겠습니다.
from numpy import *
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import sqlite3
def ang2rad(x) : return x*pi/180
def rad2ang(x) : return x*180/pi
def h2rad(x) : return x*pi/12
def rad2h(x) : return (x*12/pi+24)%24
def s2r(r, th, ph) : return r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)
def r2s(x, y, z) : return sqrt(x**2+ y**2+ z**2), arctan2(sqrt(x**2+ y**2), z), arctan2(y, x)
def rot(a, b, th) :
r = array([[cos(th), -sin(th)], [sin(th), cos(th)]])
ab= array([a.flatten(), b.flatten()])
c = dot(r, ab)
return c[0].reshape(a.shape), c[1].reshape(b.shape)
def sph(om, th, ph, x, y, z) :
y, z = rot(y, z, om) # x-axis rotation
z, x = rot(z, x, th) # y-axis rotation
x, y = rot(x, y, ph) # z-axis rotation
return x, y, z
def sph_r(ph, th, om, x, y, z) :
x, y = rot(x, y, ph)
z, x = rot(z, x, th)
y, z = rot(y, z, om)
return x, y, z
def cth(ax, ay, az, bx, by, bz) :
return arccos(dot([ax,ay,az], [bx,by,bz])/(sqrt(ax**2 + ay**2 + az**2)*sqrt(bx**2 + by**2 + bz**2)))
def read_db() :
con = sqlite3.connect("NexStarObjects.db")
cur = con.cursor()
query = cur.execute("SELECT * FROM stars;")
cols = [column[0] for column in query.description]
df = pd.DataFrame.from_records(data=query.fetchall(), columns=cols, index='Name')
con.close()
return df
def read_pic(fn='pic0.jpg', cam=[ang2rad(48.8), ang2rad(62.2)], adj=0.97, brt=194) :
fig, ax = plt.subplots()
img = plt.imread(fn)
gray = dot(img, [0.2989, 0.587, 0.144])
gray = gray/max(gray.flatten())*255
g = array(gray, dtype=uint8)
plt.imshow(img)
i = where(g>brt) # bright
iz = list(zip(i[0], i[1]))
xs=[iz[0]]
plt.text(iz[0][1], iz[0][0], str(iz[0]), color='g')
for x in iz :
if all([(abs(x[0]-y[0])>8)|(abs(x[1]-y[1])>8) for y in xs]) :
xs.append(x)
plt.text(x[1], x[0], str(x), color='g')
plt.show(block=False)
df = pd.DataFrame(array(xs), index=xs, columns=['m', 'n'])
df['brt'] = g[df['m'], df['n']]
cam = [cam[0]*adj, cam[1]*adj]
cen = (int(g.shape[0]/2), int(g.shape[1]/2))
w = sin(max(cam)/2) / max(cen)
x,_,_ = s2r(1, pi/2-(cam[0]/2), cam[1]/2)
print('w = ',w, 'x = ', x)
df['y'] = df['n']-cen[1]
df['z'] = -(df['m']-cen[0])
df['d'] = w * sqrt(df['y']**2 + df['z']**2)
df = df.sort_values(by='d')
df['arc'] = arctan2(df['d'],x)
df['d_th'] = arctan2(df['z'],df['y'])
for idx, Ay, Az in zip(df.index, w*df['y'].values, w*df['z'].values) :
l = array([])
for By, Bz in zip(w*df['y'].values, w*df['z'].values) :
if (Ay, Az) == (By, Bz) :
l = append(l, 0)
continue
l = append(l, cth(x, Ay, Az, x, By, Bz))
df[idx] = l
return g, df
def calc_pic(df, pdf) :
m, mn = [], []
for n in range(len(pdf)-2) :
m.append(list([pdf.index[n], pdf.index[n+1], pdf.index[n+2]]))
mn.append(list([pdf[pdf.index[n+1]][pdf.index[n]], pdf[pdf.index[n+2]][pdf.index[n+1]], pdf[pdf.index[n]][pdf.index[n+2]]]))
cm, es = [], []
for n in range(len(m)) :
sol, sol_es = tri_find(df.loc[:,'Acamar':'Zubeneshamali'], mn[n])
cm.append(sol)
es.append(sol_es)
#print(cm)
a = []
b = cm[0]
for i in range(len(cm)-1) :
for j in b :
for k in cm[i+1] :
if (j[-2] == k[0])&(j[-1] == k[1]) :
a.append(j + [k[2]])
ats = a
b = a
a = []
at = [(df.loc[x,'Mag'].sum()) for x in ats]
at = ats[at.index(min(at))]
if not at : at = cm[0]
av=df.loc[at[:3],['x','y','z']].values
bv = cos(pdf['arc'][:3].values)
sols = r2s(*linalg.solve(av, bv))
ptb, tb = [], []
for i in range(len(m)) : ptb.append(m[i] + mn[i])
for i in range(len(cm)) : tb.append(append(cm[i], es[i], axis=1))
return ptb, tb, at, sols
def tri_find(df, tri, e=0.01) :
ind0 = where(abs(df.values-tri[0])<e)
ind0 = array(ind0).transpose()
ind1 = where(abs(df.values-tri[1])<e)
ind1 = array(ind1).transpose()
ind2 = where(abs(df.values-tri[2])<e)
ind2 = array(ind2).transpose()
arr = []; arr2 = []
for a in ind0[:] :
if a[1] in ind1[:,0] :
temp = ind1[ind1[:,0]==a[1]]
for b in temp :
if b[1] in ind2[:,0] :
temp2 = ind2[ind2[:,0]==b[1]]
for c in temp2 :
if c[1]==a[0] :
arr.append([a[0],b[0],c[0]])
arr2.append([df.iloc[a[0],b[0]]-tri[0], df.iloc[b[0],c[0]]-tri[1], df.iloc[c[0],a[0]]-tri[2]])
y = [list(df.index[x]) for x in arr]
return y, arr2
def annot_pic(fn, rdf) :
fig, ax = plt.subplots()
img = plt.imread(fn)
plt.imshow(img)
for mn in rdf[['m','n', 'name']].values :
plt.text(mn[1], mn[0], mn[2], color='g')
plt.show(block=False)
df = read_db()
cam = (ang2rad(48.8), ang2rad(62.2)) # picamera noir v2.1 : 8MP (3280 * 2464)
adj = 0.97
#fn = 'pic15.jpg' # 12/12 stars(sirius)
#fn = 'pic11.jpg' # 6/8 stars(antares)
fn = 'pic19.jpg' # 14/18 stars(betelgeuse)
brt = 160
g, pdf = read_pic(fn, cam, adj, brt)
ptb, tb, at, sols = calc_pic(df, pdf)
pdf['name'] = ''
pdf.iloc[:len(at),-1] = at
print(pdf)
annot_pic(fn, pdf)좀 긴것 같지만 코드 170 라인 안쪽으로 기본적인 기능을 구현했습니다. 물론 상용기술과는 정확도와 알고리즘이 좀 떨어지겠지만 원리를 이해하는데 도움은 될 것입니다.
소스코드에 대한 부연설명을 하자면,,,
먼저 천체 별DB를 읽어 데이터프레임(df)으로 저장하고 촬영한 컬러사진을 명암사진으로 변환(g)하여 별의 좌표를 추출하여 저장(pdf)합니다. 구해진 좌표들을 3개씩 묶고 세변의 길이를 계산하여 ptb에 저장하고, 하나씩 길이로 오차범위(0.57º) 내의 3개의 별 짝들을 df에서 찾아 tb에 저장합니다. 각 단계에서 저장한 별들 짝이 중복되면 계속 연속해서 붙여가면서 별이름표 달기가 완성됩니다.
그러나 보완할 점이 많습니다. 초반 3개의 별 짝을 못 찾으면 시작도 못하고 에러납니다. 그래서 전처리가 중요한데, 여기서 cam, brt, 별 size, adj를 적절히 조정해서 어쨌든 초기 3개의 별 짝을 찾아야합니다. 그리고 정상적인 별 외에 노이즈들도 혼재하기때문에 초기에 pdf에서 삭제처리를 수작업으로 해야합니다.
제가 찍은 사진은 '라즈베리파이 noir 카메라' 와 '갤럭시A52S'로 야간촬영모드로 촬영을 했고, 원본사진과 별DB는 아래에 첨부하겠습니다.
그럼 최종적으로 파이썬으로 자작한 별이름붙이기(별주석) 결과물을 감상하시겠습니다.
다음은 스마트폰(A52S)으로 촬영한 사진입니다.
별들이 참 여전히 아름답습니다.
몇 만년 동안 인간은 변함없는 별들을 보며 신비해하고 우주를 탐구하려고 노력해 왔겠지요.^^;
'파이썬' 카테고리의 다른 글
| 지구에서 최단 비행경로, 위성궤도 시각화 (0) | 2023.04.10 |
|---|---|
| 로또번호 분석 & 제조기 (1) | 2023.01.14 |
| GPIB를 이용한 네트워크 분석기 활용 (0) | 2022.05.08 |
| 파이썬 GUI 실행파일 만들기(전력변환기) (0) | 2021.12.03 |
| 파이썬으로 별자리판 만들기(1) (0) | 2021.01.15 |
