2022年12月18日

皆既月食2022(その8): Hough変換による月位置合わせ

一連の皆既月食記事の8本目。いい加減にそろそろ終わりにするつもりです。前回の記事はこちらから。

Hough変換

今回はHough変換を使ったサークル抽出で、月食中の月の位置合わせに挑戦です。FS-60CBで撮影した4時間分の広角の画像から月が画面中心に来るように位置変換をします。その後、タイムラプス映像にしてみます。

Hough変換は画像の中から特徴的な形を抽出するアルゴリズムで、その中に円を抽出する関数があります。Hough変換は色々な環境で使えますが、今回はOpenCVで用意されている関数をPythonで使うことにしました。検索するとサンプルプログラムなどたくさん出てくるのと、今回は対処療法で組んでいったのの積み重ねなので、あまりに汚いコードで人様に見せれるようなものではないです。なので、どうパラメータを取ったかだけの説明にとどめることにします。

基本的にはあるフォルダにある、月が映っているたくさんのファイルを順次読み込み、カラー画像をグレースケールに変換して、HoughCircles関数を呼び出すだけです。Hough関数は以下のようにしました。

HoughCircles(gray, cv2.HOUGH_GRADIENT, dp=1, minDist=1500, param1=30, param2=2, minRadius=350, maxRadius=400)

パラメータがいくつかありますが、少しだけコツを書いておきます。

HOUGH_GRADIENTとHOUGH_GRADIENT_ALTは両方試しましたが、ALT付きの方が誤検出が多かったのでALT無しの方にしました。
dpは1以下も試しましたが、位置精度に違いはあまりなかったです。大きな値にすると精度が悪くなりました。
円を複数検出するわけではないので、minDistは1000とかの相当大きな値にしておきます。
param1は最終的に50程度にしました。大きすぎたり小さすぎたり値では検出できなかったりしますが、位置精度にはあまり影響ないようです。
param2は位置精度に影響があるようです。大きすぎると精度が悪くなりますが、5以下くらいだと精度はこれ以上変わらないようです。
minRadiusとmaxRadiusは、月の大きさは一定なのでその半径を挟むような値を取ると効率が良いでしょう。

とりあえず結果を示します。

結局出た精度はこれくらいです。もうちょっとビシーっと止まってくれればいいのですが、まだちょこちょこズレています。

実は精度を出そうとして色々試してみてます。Hough変換のパラメータをいじるのはもちろんなのですが、大きなものは

グレー化後に、あらためて画像の2値化
グレー化後に、あらためて輝度やコントラストをいじる
Adobe Premiereのブレ補正
一旦Hough変換でラフに位置合わせして、少し大きめに切り抜いて、再度Hough変換

くらいでしょうか。でも結局どの方法も精度を劇的に上げることはできませんでした。Hough変換で目で見てもおかしいズレがあるので、そこだけでもずれが直ったらと思ったのですが、やはりダメなものはダメで、どうも苦手な画像があるようです。

今回はFS-60CBで広角で撮ったものから抜き出しているので、解像度が良くないこと。さらに動画の中から最初の1枚だけを抜き出しているので、大気揺らぎなどで多少のブレがあります。もしかしたらスタックをすると平均化されるので、細かいブレは少なくなるのかもしれません。スタックしてから、Hough変換ではなく、特徴点を抜き出して合わせるとかいう処理をした方が精度が出るかも知れませんが、ちょっと力尽きたので、今回はここまでにします。

ついでに同じ位置合わせの手法で、いくつか画像を抜き出して天王星の潜入画像をつくってみました。位置合わせ後、比較明合成しただけです。これくらいの精度ではそこそこ合って見えるんですけど、やはり長時間のタイムラプスだときついかもです。

まとめ

今回の皆既月食に関する記事はこれで一旦終わりにします。

撮影後、1ヶ月以上にわたって楽しむことができてかなり満足しました。その一方、そろそろ月も飽きてきました。ただ、まだ未処理ファイルが大量にあるので、気が向いたらもう少し追加で書くかもしれません。

先の記事でも書きましたが、今回は月食に関してはかなり満足して撮影できました。大きな課題は地球本影を位置補正することなく撮影することですが、次回月食ではこれ一本に絞ることにするかもしれません。

ここ数日、BlurXTerminatorがすごいことになっています。次は少しこちらの方に時間を費やそうかと思います。

タグ：: #月; #月食; #Hough変換; #位置合わせ; #皆既月食; #天王星食; #タイムラプス; #python; #OpenCV

2020年03月21日

星像切り出しスターベース仕様、Matlabバージョン、多画像形式対応

以前、星像切り出し用のコードをPythonで書いたのですが、

Matlab用に少し書き直したので公開します。結果はすでにここ最近の記事で使っているので、これまでに見ている方も多いかと思います。デザインはピンとくる方もいるかと思いますが、スターベース東京のブログの作例の切り出し画像のフォーマットに合わせてあります。

ファイル選択ですが、以前選んだフォルダを覚えるようにしました。以下のページを参考にしました。ありがとうございました。

GUIでのファイル選択について

前回の対応フォーマットはJPEGだけでしたが、今回はMatlabのimreadコマンドでサポートする画像ファイルに対応しています。なので、かなりのフォーマットに対応するはずです。実際に全部試したわけではないので分かりませんが、

BMP、JPEG、PNG、CUR、JPEG 2000、PPM、GIF、PBM、RAS、HDF4、PCX、TIFF、ICO、PGM、XWD

に対応するらしいです。天文で関連するのは主にJPG, PNG, TIFFくらいでしょうか。RAWファイルはあまり対応できないのですが、個別にFITS形式だけ対応させておきました。

カラーの8bit、16bitに対応しています。32bitには対応していません。あと、グレー画像は多分ダメです。結果はファイル名に”_cut25”というのが足されて、元のファイル形式と同じ形式（例えばjpgならjpg、tifならtif）に書き出されます。

あまり綺麗でないですが、ソースコードです。コピペして、Matlab上で走らせてみてください。上のような画像が出てくるはずです。簡単なコードなので、各自で希望に応じて適当に書き換えてみてください。3x3マスとかも簡単にできるはずです。

clear;

%%% Paramters
CS = 300;
BW = 5;

% File select
if ispref('MyPreferences','LastUigetfileFolder')
    folder = getpref('MyPreferences','LastUigetfileFolder');
    if ~ischar(folder)
        folder = '/Users/'; % for mac. if windows use 'C:\'.
    end
else
    folder = '/Users/'; % for mac. if windows use 'C:\'.
end
[file,path] = uigetfile({  '*.*',  'All files(*.*)'}, 'Pick a file','MultiSelect', 'on',folder);

if ~isnumeric(path)
    setpref('MyPreferences','LastUigetfileFolder',path)    
end

% Reading image with size and class (full size)
[filepath,name,ext] = fileparts(file);

if extractBefore(ext,5) == '.fit'
    Img = fitsread(fullfile(path, file));
else
    Img(:,:,:) = imread(fullfile(path, file));
end
[y, x, l] = size(Img);
if isa(Img,'uint16')
    cv = 256;
else
    cv = 1;
end

% Size of ASP-C
ay = round(size(Img,1)/1.6);
ax = round(size(Img,2)/1.6);

% empty cut image
CutImg = zeros(5*CS+6*BW, 5*CS+6*BW, 3, class(Img));

% Orange area
CutImg(:,:,1) = 235 *cv;
CutImg(:,:,2) = 170 *cv;
CutImg(:,:,3) = 80 *cv;

% Blue Area
BA = CS+BW+1:4*CS+5*BW;
CutImg(BA,BA,1) = 50 *cv;
CutImg(BA,BA,2) = 50 *cv;
CutImg(BA,BA,3) = 80 *cv;

% Define areas in cut image
C = zeros(CS,5);
for i = 1:5
    C(:,i) = (i-1)*CS+i*BW+1:i*(CS+BW);    
end

% Area in original image
IX(:,1) = 1:CS;
IX(:,2) = round((x-ax)/2)+1:round((x-ax)/2)+CS;
IX(:,3) = round((x-CS)/2)+1:round((x+CS)/2);
IX(:,4) = round((x+ax)/2)+1:round((x+ax)/2)+CS;
IX(:,5) = x-CS+1:x;

IY(:,1) = 1:CS;
IY(:,2) = round((y-ay)/2)+1:round((y-ay)/2)+CS;
IY(:,3) = round((y-CS)/2)+1:round((y+CS)/2);
IY(:,4) = round((y+ay)/2)+1:round((y+ay)/2)+CS;
IY(:,5) = y-CS+1:y;


% Filling cut image by original cut
for i = 1:5
    for j = 1:5
        if ( ((j==1)||(j==5)) && ((i==2)||(i==4)) ) || ( ((j==2)||(j==4)) && ((i==1)||(i==5)) )
            CutImg(C(:,i),C(:,j),:) = 256 *cv; % fill with white
        else
            CutImg(C(:,i),C(:,j),:) = Img(IY(:,j),IX(:,i),:);
        end
    end
end


%image(CutImg);
if extractBefore(ext,5) == '.fit'
    fitswrite (CutImg, fullfile(path, [name,'_cut25',ext]));
else
    imwrite (CutImg, fullfile(path, [name,'_cut25',ext]));
end

途中からいろいろ簡略化しているのでわかりにくいと思います。% Define areas in cut image以降を、以下のコードに置き換えるとまだわかりやすいかと思います。

% Define areas in cut image
C1 = 0*CS+1*BW+1:1*CS+1*BW;
C2 = 1*CS+2*BW+1:2*CS+2*BW;
C3 = 2*CS+3*BW+1:3*CS+3*BW;
C4 = 3*CS+4*BW+1:4*CS+4*BW;
C5 = 4*CS+5*BW+1:5*CS+5*BW;

% Area in original image
IX1 = 1:CS;
IX2 = round((x-ax)/2)+1:round((x-ax)/2)+CS;
IX3 = round((x-CS)/2)+1:round((x+CS)/2);
IX4 = round((x+ax)/2)+1:round((x+ax)/2)+CS;
IX5 = x-CS+1:x;

IY1 = 1:CS;
IY2 = round((y-ay)/2)+1:round((y-ay)/2)+CS;
IY3 = round((y-CS)/2)+1:round((y+CS)/2);
IY4 = round((y+ay)/2)+1:round((y+ay)/2)+CS;
IY5 = y-CS+1:y;

% Filling cut image by original cut
CutImg(C1,C1,:) = Img(IY1,IX1,:);
CutImg(C2,C1,:) = 256 *cv;
CutImg(C3,C1,:) = Img(IY3,IX1,:);
CutImg(C4,C1,:) = 256 *cv;
CutImg(C5,C1,:) = Img(IY5,IX1,:);

CutImg(C1,C2,:) = 256 *cv;
CutImg(C2,C2,:) = Img(IY2,IX2,:);
CutImg(C3,C2,:) = Img(IY3,IX2,:);
CutImg(C4,C2,:) = Img(IY4,IX2,:);
CutImg(C5,C2,:) = 256 *cv;

CutImg(C1,C3,:) = Img(IY1,IX3,:);
CutImg(C2,C3,:) = Img(IY2,IX3,:);
CutImg(C3,C3,:) = Img(IY3,IX3,:);
CutImg(C4,C3,:) = Img(IY4,IX3,:);
CutImg(C5,C3,:) = Img(IY5,IX3,:);

CutImg(C1,C4,:) = 256 *cv;
CutImg(C2,C4,:) = Img(IY2,IX4,:);
CutImg(C3,C4,:) = Img(IY3,IX4,:);
CutImg(C4,C4,:) = Img(IY4,IX4,:);
CutImg(C5,C4,:) = 256 *cv;

CutImg(C1,C5,:) = Img(IY1,IX5,:);
CutImg(C2,C5,:) = 256 *cv;
CutImg(C3,C5,:) = Img(IY3,IX5,:);
CutImg(C4,C5,:) = 256 *cv;
CutImg(C5,C5,:) = Img(IY5,IX5,:);

タグ：: #Matlab; #画像; #切り出し; #星像; #比較

2019年02月07日

お手軽プログラミング: 四隅の切り抜き図を作る

画像処理で出来上がったファイルの四隅の星像を比べたくなることがあると思います。これまでPhotoshopで手作業で切り取り貼り付けを駆使して作っていたのですが、何度もやるとなるとめんどくさいので、簡単にできるようにpythonとOpenCVを使って作ってみました。

超お手軽プロブラミングなので、内容はソースを見たらすぐにわかると思います。唯一工夫したところはファイル名をGUIで選択するところくらいでしょうか。Macで試しましたが、多分Windows他でも動くと思います。（すみません、チェックもしてません。）

必要な環境

環境を構築するのもいたって簡単です。今回はMacで試したので、Macでやる場合を書きます。WindowsなどでもPythonとOpenCVが動く環境なら走るはずなので、難しくもなんともないと思います。

1. Macの場合、Phython3はHomebrewからインストールするといいでしょう。Homebreはターミナルを立ち上げて、

/usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"

でインストールできます。

2. 次にPython3です。これもコマンド一発で、

brew install python3

3. 最後はOpenCVのインストールです。

pip3 install opencv-python

でインストール完了です。xcodeをインストールしていない場合は、最初にそれも必要かも。うまくいかない場合は、私なんかの説明よりも、適当に検索すればOpenCVが使えるようになるまでの説明はすぐに出てくるので、探してみてください。

ソースコード

さて、python3とOpenCVの環境が整ったら次はプログラミングです。今回はホントにチャカチャカっと書いたコードなので、全然綺麗ではありませんのでご容赦ください。ファイル選択の際文句が出たりもしますが、とりあえず動きます。

エディタなどで以下のコードをテキストファイルにコピぺして、適当なファイル名をつけて、拡張子を.pyにすれば準備完了です。

import cv2

import numpy as np

import tkinter

from tkinter import messagebox as tkMessageBox #python3

from tkinter import filedialog as tkFileDialog #python3

def pathname(fullpath):

n = fullpath.rfind('/') + 1

return fullpath[:n]

def select():

root=tkinter.Tk()

root.withdraw()

fTyp = [('JPEG file','*.jpg')]

iDir = '/Users/ユーザー名/' #自分のユーザー名を入れてください

filename = tkFileDialog.askopenfilename(filetypes=fTyp,initialdir=iDir,title = "ファイル選択")

iDir = pathname(filename)

return filename

filenameselect = select()

img = cv2.imread(filenameselect)

mlen = 250

size = img.shape

roi_UL = img[0:mlen, 0:mlen]

roi_UR = img[0:mlen, size[1]-mlen:size[1]]

roi_LL = img[size[0]-mlen:size[0], 0:mlen]

roi_LR = img[size[0]-mlen:size[0], size[1]-mlen:size[1]]

roi_CC = img[round((size[0]-mlen)/2):round((size[0]+mlen)/2), round((size[1]-mlen)/2):round((size[1]+mlen)/2)]

width = mlen*3

height = mlen*3

imageArray = np.zeros((height, width, 3), np.uint8)

count = np.array([[0,0],[0, mlen*3],[mlen*3, mlen*3],[mlen*3, 0]])

cv2.fillPoly(imageArray, pts=[count], color=(255,255,255))

imageArray[0:mlen, 0:mlen] = roi_UL

imageArray[0:mlen, mlen*2:mlen*3] = roi_UR

imageArray[mlen*2:mlen*3, 0:mlen] = roi_LL

imageArray[mlen*2:mlen*3, mlen*2:mlen*3] = roi_LR

imageArray[mlen:mlen*2, mlen:mlen*2] = roi_CC

savefile = filenameselect.replace('.jpg', '_4cut.jpg')

cv2.imwrite(savefile,imageArray)

実行結果

ターミナル上で上のプログラムがあるフォルダに行き、

python3.py ファイル名.py

とかで実行すると、画像ファイルを選択するダイアログが現れます。ダイアログが前面に出てこないことがあるので、下のDocからPythonというアイコンを選んでみてください。JPEGファイルしか選べませんが、プログラムの中の15行目を*.pngとかすれば多少他のファイル形式も選べます。

切り取りたいファイルを選択すると、同じフォルダの中に

オリジナルファイル名_4cut.jpg

というファイルが出来上がります。

出来上がりは下のようになります。元の画像から四隅250x250ドット分と真ん中を250x250ドット分を切り取って750x750ドットのファイルを作ってくれます。ホントにこれだけです。

ちなみにこれは、タカハシの新フラットナーで撮ったFS-60CBの星像です。カモメ星雲ですが現在画像処理中のものです。ほとんど流れていないことがわかると思います。元々の動機が、新旧フラットナー、レデューサー、エクステンダーでの四隅の星像を比較したくて、いちいちやるのがめんどくさいのでプログラミングに走ったというわけです。

今回のプログラムは色々制限もありますが、簡単なプログラムなので中身もすぐにわかると思います。各自で適当に改良してみてください。

例えば、26行目の250を他の値に変えると切り出してくれる画像の大きさが変わります。
ちょっと頑張れば上下左右も合わせて8方向の切り取りとかもできると思います。

画像処理プログラミングの取っ掛かりとしては、いい練習になるかと思います。

(追記: Matlabバージョンの切り出しプログラムも書いています。スターベースのブログと同じ形式で5x5マスでフルサイズとAPS-Cサイズを見るものです。よかったらこちらもどうぞ。)

参考ページ

以下のページを参考にしました。

Honebrew: MacにHomebrewを導入する方法
ブランク画像: OpenCV3とNumpyでブランク画像を作る
OpenCV: 画像上の基本的な処理
ファイルダイアログ: Python3でファイルを選択するダイアログボックス（初心者向け）

どのページもわかりやすくて、とても役に立ちました。どうもありがとうございました。

ほしぞloveログ

天体観測始めました。

カテゴリ：計算機 > プログラミング

皆既月食2022(その8): Hough変換による月位置合わせ

星像切り出しスターベース仕様、Matlabバージョン、多画像形式対応

お手軽プログラミング: 四隅の切り抜き図を作る