星を写した画像を用いてカメラの歪みの検出し、
それを用いて画像加工をするのためのプログラム群です。

=====================================
解凍後の内容物
◎プログラム本体
noisec.exe:固定ノイズを引き算して除去する。
phtlnr.exe:写真上の星の位置と、指定した図法での星図上の位置から、
           その相互関係をアフィン変換で近似し、その結果と誤差を
           テキストと画像で表示。
           さらにそれらの計算結果を用いて、元の画像に合成するための
           赤経赤緯線や星座線を描画した画像を作成する。
phtply.exe:同じカメラ（レンズ、倍率）で撮影した複数の写真に対して
           それぞれsphtlnr.exeを適用して得られたデータを元に、
           そのカメラの歪み具合を多項式近似する。
suming.exe:ほぼ同じ領域をちょっとずつずらして撮影した写真を、
           phtlnr.exeとphtply.exeで得られたデータを元にして
           １枚の画像に重ね合わせて、それによってランダムノイズを
           消去したり暗い星を浮かび上がらせたりする。
           または複数の画像を貼り合わせて広域の画像を作る。
compos.exe:phtlnr.exeで作成した赤経赤緯線や星座線の画像を
           元の画像に上乗せする。
F1A0RI63.DLL,F1A0RN63.DLL,F1A0RW63.DLL:DLLファイル(CLRCON.exeと共通です)
◎関連データ
ET-UTC.txt:世界時(UTC)と暦表時(ET)・力学時(TD)の差に関するデータ
zenten5.txt:星表データ(5.5等級より明るい星,星座所属等の情報付)
senzanet.net:星座線データ(赤経赤緯表示)
◎サンプルデータ
rdn(番号).jpg:ブダペスト市内で撮影した北斗七星の写真。
              都心部で月齢16の悪条件下です。
rdn(番号).txt:phtply.exe用の各写真用データ。撮影地、日時などの基本データ、
            上記写真から読み取った星の写真上の座標、
            その星に対応するzenten5.txt中での星表番号のリスト。
blk4s.BMP:上記写真の撮影に使ったデジカメの固有ノイズデータ（夜中に布団の中で撮影したもの）。
time.txt:上記写真の撮影日時(デジカメによるタイムスタンプを、ls -l からリダイレクトで取得したもの)
         及び関連するデータ、メモ
noisec-t.txt:noisec.exe用のテストスクリプト
phtlnr-t.txt:phtlnr.exe用のテストスクリプト
phtply-t.txt:phtply.exe用のテストスクリプト
sumimg-t.txt:sumimg.exe用のテストスクリプト
=====================================
◎基本的な操作の流れ
1)写真に固定ノイズが強く残っている場合はそれを除去する(noisec.exe)
2a)写真に写っている星をhayami.exeなどで同定し、
テキストファイルで「写真での座標」と「赤経赤緯(または星表番号)」を
対応付ける。手作業で:-)。
2b)写真と指定した図法(通常は正距方位図法)の座標の間のアフィン変換を
最小自乗法で求める(phtlnr.exe)
3)複数の写真で2)の作業を行った後で、そこで得られたデータを元に、
正距方位図法と写真での写り方のずれを最小自乗法で多項式近似する(phtply.exe)
この時点でカメラの歪みの近似が出来上がる。

4)ユーザー定義正距図法の定義式として4)の結果を用いて、
元の写真に合わせて星座線や赤経赤緯線を引いた画像を作成し(phtlnr.exe)、
元の写真と合成する(compos.exe)。
5)複数の画像で1)-3)の作業を行い、
そのデータを用いてそれらの画像を重ね合わせる(sumimg.exe)

=====================================
各ソフトの使用方法
1)noisec.exe
このソフトは高感度撮影の際に出る、カメラ固有のノイズを除去するものです。
したがって、カメラ固有ノイズのデータが必要になります。
ただしこのデータの入手方法は簡単です。星空を撮影する時と同じ高感度撮影設定にして、
夜中に電灯を消して布団の中で撮影するだけです。多少条件を変えて数枚撮影しましょう。
これをMS-Windows付属のpaintなどでフルカラーBMPに変換し、適当な名前を付けます。
(サンプルデータではblk4s.BMP)。

多数の写真をまとめて処理するためにスクリプトで動作させる事も出来ます。
スクリプトの形式は、サンプルのnoisec-t.txtで分かるように単純で、
１行につき

(オリジナル画像名[jpg,bmp]),(処理後の保存画像名[bmp])

として、並べていくだけです。
「処理後の保存画像名」は省略できます。
「処理後の保存画像名」を省略した場合の画像名は
a)jpgであればオリジナル画像の拡張子を.bmpに変えた名前
b)bmpの場合はオリジナル画像名の.bmpの直前に"c"を挿入名前
です。
なおjpgの場合はオリジナル画像を一旦bmpに変換する必要があるために、
(オリジナル名).bmpが生成されます。処理後の画像を同じ名前にするのは構いませんが、
同じ名前の画像が別にある場合に断り無く上書きします。
オリジナル画像がbmpの場合は、同じ名前で上書きする事も出来ます。

（スクリプトを使う場合はそれを作成した後）noisec.exeを起動します。
まずノイズデータ画像の名前を入力し、続けてにノイズ除去率、除去閾値を入力します。
この２つの数値の設定は画像の状態次第や目的で変わってきます。
都会など背景が明るい写真の場合は、除去率が高いと黒い穴が目立ってしまいます。
また閾値が低いと、暗い星がノイズと一緒に消されてしまう恐れがあります。
一方で除去率が低くかったり閾値が高かったりすると、星とノイズの判別が
難しい事もあります。判別の際には除去率を高めて確実な星だけを判別し、
多くの枚数を重ね合わせる時は除去率を下げて、重ね合わせによるノイズキャンセル効果に
任せる、など使い分ける方がよいでしょう。

次にスクリプトを使うかどうかを入力します。
スクリプト使用の場合はスクリプトファイル名を入力すれば終わりです。
スクリプトを使用しない場合は、ひとつずつオリジナル画像名と出力画像名を入力します。
何も入力しないでEnterを押したら終了です。

-----------------------------------
2)phtlnr.exe(基本データ収集編)
このソフトは、画像と選択した図法との間の関係を最小自乗法で線形（アフィン変換）
近似する行列を求め、画像中心点の仰角方位角・赤経赤緯・カメラの傾き、図角の大きさ、
近似式での各星における誤差などを計算し、テキストと画像で吐き出します。
誤差情報画像の代わりに赤経赤緯線と星座線を描いた画像を生成し、
compos.exeでオリジナルと合成する事も出来ますが、これは後述します。

a)まずは計算をするために必要な情報をテキストファイルに列挙していきます。
サンプルのrdn(画像番号).txt(同じ名前のjpgに対応)をコピー＆ペイストして
必要な部分だけ書き換えればいいのですが、一応述べておくと、定幅形式で次の形です。
定幅形式なので桁数などは変えないでください。

date :2007/03/05,07:07:14  <-年月日時分秒
zone :+09h00m              <-タイムゾーン
place:+019deg05am25as,+47deg30am22as,10m <-撮影地の経度緯度(60進表示)と標高(m,km)
size :( 640, 480) <-画像のサイズ（横,縦）
( 496, 306)UMa0006              <-星の画像上の座標と対応する星表番号
( 424, 263)[200.9809,+54.9254]  <-星表番号の代わりに赤経赤緯を度単位(10進小数)表示
…(以下星の数だけデータを並べる)…
#------------------------       <-終端

時刻はデジカメのタイムスタンプを見ればいいのですが、その場合のタイムゾーンは
デジカメの時計のタイムゾーンにする必要がある点に注意します。
日本時間に合わせたデジカメを、設定を変えずに外国に持っていった場合に注意が必要です。
標高は、普通は0mで問題ありません。よほど高い山から月でも撮るんでなければ関係ないです。

あとは、paint辺りで画像を拡大表示して、せっせと画像での座標を読み取ります:-)
CLRCON.exe ならば、拡大して星の上でクリックすればクリップボードから
すぐに貼り付け出来ます（書式も合わせてあります）。
星がpixelの境界上に来たり、分解能ギリギリの二重星が写っている場合に
同定が難しい事も少なくありませんが、その時は「適当に読み取ってください」。
適当に読み取っても誤差は大きくて1pixel程度です。
星の数を多くすれば、最小自乗法処理で実質的に誤差がキャンセルされます。
そして、それぞれの星に対応する星表番号か赤経赤緯を記述します。
hayami.exeのクリック検索機能で「星表番号表示」に設定し、クリックすれば星表番号が
出ますし、暗い星で星表番号が出ない場合でも赤経赤緯が出ます。
一応３点とればアフィン変換行列は求まりますが、カメラの歪みまで検出するために
写真１枚につき１０点程度は取っておきましょう。


b)写真と星の対応を記述したテキストファイルが出来上がった後、phtlnr.exeを起動します。
まずは比較対象とする図法を選択します。
といっても、1番と10番以外は技術的な興味がなければほとんど意味がありません。
また、10番はカメラの歪みの多項式近似が出来た後でないと意味がないので、
実質1番(正距方位図法)だけ意味があります。

次に、スクリプトを使用するかどうか尋ねてきます。
使用する場合は、以降に述べる全ての誤差表示機能がオンになりますが、
その機能を全て述べてからスクリプト作成について述べます。

スクリプトを使用しない場合、作成したテキストファイルの名前を入力し、
計算結果を出力するテキストファイル名を入力します。
次にビットマップ名を入力するのですが、これは誤差の程度を一目で確認するためのものです。
テキストファイルだけでも分かるので、必要ない場合はEnterだけを押して終了します。
ビットマップ名を入力すると、次に「1:check picture,2:overdraw picture」の選択肢が出ます。
ここでは誤差情報確認ですから、1を選択します。以下、描画する内容を尋ねてきます。

b1)「line for checking」は、画像の中心を中心とする等間隔の円と、縦と横の大円線を描画します。
経緯度線ではなく大円線です。心射図法で直線になる線であり、地表物ならば直線からなる図形が
写真上でどのくらい歪むか見るためのものです。
経線は大円なので同じようなものですが、緯線は大円ではないので違ってきます。

b2)「Do you draw stars?」で「y」を選択すると、写真上の星の位置に小さな円を描き、
アフィン近似での計算上の位置との誤差を強調して直線で表示します。
その強調倍率を次に尋ねてきます。１倍にすれば、計算上の位置と写真上の実際の位置を
結ぶ線になりますが、それでは小さ過ぎて見えないので、10倍くらいを選びます。

b3)「Data of Errors?」で「y」を選択すると、誤差を動径方向成分と方位角方向成分に分解して、
出力画像左下に、中心からの距離を横軸にしてグラフを描きます。
アフィン変換（線形近似）からどのくらいずれているか見るためのものです。
こちらも強調倍率は10倍くらいが適当です。

多数の画像情報を処理する場合のスクリプト(サンプルではphtlnr-t.txt)は、１行ごとに

(画像と星の対応情報ファイル),(出力テキストファイル),(画像).bmp,誤差拡大倍率(２種類の倍率共通)

の４つを並べたものを、画像の数だけ並べます。画像での選択肢は３つとも描画で、
２つの拡大倍率は共通値とします。スクリプトを使用する方を選んだ後にスクリプト名を入力すると、
後は自動的に計算結果テキストと誤差情報画像を１画像につきひとつずつ出力します。


c)これで１画像につき、テキストファイルと画像(bmp)ファイルの２つを出力しますが、
まずは画像ファイルを簡易ビュアで見ましょう。
星の位置を示す直線が、全て画像の大きさの１割程度かそれ以下ならば、特に問題はありませんが、
極端に長いものがあれば、座標読取か星の同定が大きく間違っている可能性があります。
極端に大きな誤差がひとつあると、正しく入力した星にも逆向きに数pixel〜十数pixelの
誤差が発生します。ですから全体的に長い線が多数出て、「極端に大きな誤差」が見えにくくなります。
そういう場合は、テキストファイルの方を見て、誤差が極端に大きいもの（普通はひとつかふたつ）
を見つけて、もう一度座標と星の同定を確認修正して、もう一度phtlnr.exeを通します。

出力テキストの方には、画像中心点の方位角仰角、赤経赤緯、垂直や赤経線に対するカメラの傾き、
写真座標と図法投影平面との近似アフィン変換行列とその逆行列、Ｘ軸やＹ軸のずれ、
縦軸と横軸の長さのずれ、対角線の図角などが記述された後に、
各点での写真上でのずれの長さ、図法投影平面上でのずれの動径方向、方位角方向の各成分、
それらの平均値が表示されます。
ずれの平均値は、線形近似であるために、通常1pixelを上回ります。
正距方位図法の場合に2pixelを越える事は少ないと思いますが。

なお、写真中心点の方位角・仰角・傾きは、記述された星の星表位置を、写真上の距離の逆数で重みを
付けた加重平均を計算して仮の中心とし、最小自乗法でアフィン変換をとりあえず求めます。
仮の中心を元に計算しているのでアフィン変換式で中心がずれるはずですから、
そこで中心を補正する。これを６回繰り返して求めたものです。再帰的な近似計算なので、
若干誤差が残りますが、１万分の１〜１億分の１pixel程度なのでまあいいかと。
Ｘ軸の傾きは同様の計算をするので十分小さくなりますが、Ｙ軸は単純なアフィン近似の結果なので、
画像そのものの歪みを表しています。同じく縦横比も画像の歪みを表しています。
-----------------------------
3)phtply.exe
このソフトは、2)で作成した各写真のデータを基にして、カメラの歪みを多項式近似します。
カメラのレンズは、スキャン型など特殊なものを除いて中心軸に関して対称なので、
方位図法になると想定して、中心軸からの角度と写真中心からの距離の関係を最小自乗法で
多項式近似します。

このソフトでもスクリプトで動作させる事が出来ます。phtply-t.txtを見ればお分かりの通り、
phtlnr.exeで出力したテキストファイルの名前を、１行に１個ずつ書いて行くだけです。
ただしここで列挙するファイルは、同じ条件で撮影した写真に関するファイルだけを並べます。
「同じ条件」というのは、同じカメラ、同じ画像サイズ、同じレンズ、同じ倍率くらいです。
写真の歪み方は、レンズや倍率が違えば当然違ってきますから、違うものを並べても意味がありません。
同じズームレンズの倍率を変えた場合に関しては、検討の余地があるかも知れませんが、
それを知るためにも個々の倍率のデータを比較する必要があります。
もしかしたら、何度も倍率変更を繰り返していた場合に、その都度微妙に変化するかも知れませんが、
どの程度の精度で歪み検出したいか、どの程度条件を揃えられるか、などを考慮して
適当に決めてください。画像サイズに関しては、うちのデジカメでは解像度を変えて
画像サイズが変わっても、同じ光学レンズ倍率ならば同じような近似式が出ます。
しかし全てのデジカメでそうなるかは分かりません。
なおスクリプト動作でなければ、１つの写真からのデータだけで計算します。
その場合は2)の時点で十分な点の数を確保しておいてください。

phtply.exeを起動すると、まず
「多項式の係数を奇数次だけ非零にして偶数次は０にするかどうか」
の選択肢が出てきます。明確な根拠があるわけではないのですが、奇数次だけの方がいいような
感じがしないでもないです。ま、何通りか試してください。
次に多項式の最高次数を尋ねてきます。正直あまり大きくしても意味がありません。
５次くらいでいいと思います。３次で十分かも知れません。
そしてスクリプトを使用するかどうかを選択し、使用する場合はスクリプト名、
使用しない場合はphtlnr.exeの出力ファイル名を入力します。
そして出力テキスト名と出力画像名を入力したらおしまいです。

まずは出力画像を見てみましょう。画像の対角線の半分の長さを一辺とする正方形の画像です。
2)で採取した点での「正距方位図法からのずれ」が、中心点からの角距離をＸ軸として
プロットされ、最小自乗法で求めた近似式の曲線が合わせて表示されます。
出力テキストの方には、方位図法での角距離と図法平面での中心からの距離との関係を
近似する多項式係数、その多項式を用いた方位図法と写真との間のアフィン変換行列があり、
他にphtlnr.exeの時と同じようなデータが並びます。ただし各点の誤差ベクトルは、
通常の直交座標ではなく、動径方向と方位角方向に分解して表示してあります。
細かいレンズの歪み具合を知りたいのでなければ細かい数値はどうでもいいんですが、
このテキストファイルを後述する各種ソフトで使います。
------------------------------------------
4)phtlnr.exeを使って、オリジナル画像に合わせて方位角仰角、赤経赤緯線、黄道線、星座線を
描画し、compos.exeで重ね合わせる事が出来ます。

a)まずphtlnr.exeを起動します。2)の時と同様に図法選択が出てきますが、
既にレンズの歪みは検出済みなので、「10:ユーザー定義方位図法」を使えます（他のでもいいですが）。
10を入力すると、多項式係数ファイルの入力を求めてきますが、これは3)のphtply.exeで
出力したファイルです。その名前を入力します。
次にスクリプト使用について尋ねてきますが、これは誤差検出専用なので「n」を選択します。

次に2)において手作業で作成した、星の写真位置と星表位置を対応させたファイル名を入力します。
この情報を元に描画を行います。
次に出力テキストファイル名を入力します。これは2)と同じく誤差情報などを記述したもので、
本当はあってもなくてもいいものですが、とりあえず出力します。

次に出力画像名を入力し、その次に「2:overdraw picture」を選択します。
以下、天頂座標線、赤道座標線、黄道、星座線、それぞれを描画するかどうか尋ねてきます。
天頂座標線と赤道座標線の間隔は、対角線の画角と画像そのものの大きさで自動的に決めます。
黄道線は黄道(黄緯0度線)だけです。
星座線を描く場合は、星に重なって描画しないように空白域(というか黒なんですけど)を設けるので、
その空白域の半径を入力します。
以上で重ね合わせ用の画像が完成します。

b)composit.exeは、同じ機能を持った画像処理ソフトが多数あるでしょうが、
自分が持ってなかったから作っただけです:-)
起動して、オリジナル画像名（星の写真,jpgでもbmpでも可）を入力し、
4a)で作成した重ね合わせ用画像名を入力し、出力画像名を入力して終わりです。
重ね合わせ側の黒((R,B,G)=(0,0,0))を透明として重ね描きします。
オリジナル画像がjpgの場合は作業用に同名.bmpを作成します。
------------------------------------------
5)sumimg.exeは、複数の画像を重ね合わせてランダムノイズを除去し暗い星を浮かび上がらせる、
または複数の画像を貼り合わせて大きな画像を作るためのソフトです。
画像と関連情報を読み込み、画像の歪みを除去して赤経赤緯を計算し、
同じ赤経赤緯になる点をRGBそれぞれ足し合わせていきます。当然255を越える値になるので、
適当に割り算をして、画像を出力します。

このソフトはスクリプト使用だけで使います。起動後はスクリプト名を入力するだけです。
そのスクリプトは次のような構成になります。

projection factor:phtply-u.txt <-多項式係数ファイル,3)での出力テキスト
brightness factor:1.1          <-画像全体で最も明るい点の明るさの調整係数
output image name:sumimg-t.bmp <-出力画像名
output image size:640          <-出力画像の大きさ（このpixelの正方形）
output projection:1            <-出力画像での図法(番号)
angle of scope :24.56          <-画像の横幅の半分に相当する視野角
right ascension:178            <-画像中心の赤経(度単位,10進表示)
declination    :+60            <-画像中心の赤緯(度単位,10進表示)
rotation,clockw:245.5          <-真上を基準として天の北極が向く方向(回転角)
function:log                   <-linear,log,expのいずれか
factor1 :2                     <-係数1
factor2 :-80                   <-係数2
factor3 :266                   <-係数3
input images                   <-ここ以降が入力画像に関するデータ
rdn01543.jpg,rdn11543.txt      <-画像名(jpg or bmp),2)でphtlnr.exeが出力したファイル
rdn01551.jpg,rdn11551.txt
…(以下画像の数だけ続ける)…

output projection の番号は
 1.正距方位図法
 2.正積方位図法
 3.正角方位図法
 4.メルカトル図法
 5.縦メルカトル図法
 6.平射図法
 7.モルワイデ図法
 8.正射影図法
 9.サンソン図法
10.(欠番)
11.正距円筒図法
12.ミラー図法
13.ハンメル図法
14.ランバート正積円筒図法
15.アイトフ図法

function,factor1,factor2,factor3 は CLRCON.exeと共通です。
function:linear の時は、Y=(factor1)X+(factor2)で変換します。
function:log,exp の時は,明るさ 0 が factor2 に、明るさ最大値が factor3 になります。
factor2<factor3(単調増加)の時は,log が上に凸,exp が下に凸になります。
factor1が大きければ凸の度合いが大きくなります。

省略した場合はfunction=linear,factor1=brightness factor,factor2=0 とします。

functionを省略した時に有効になる brightness factor は、
1とした時に一番明るい点のRGB値が255になるように調整します。
1より小さくすれば暗くなります。1より大きくした時は、明るめの星・色がすべて「255」になり、
その分暗い星が強調されます。ただし背景も明るくなる可能性があります。


画像の大きさは、オリジナルの画像から見て等倍程度が良いでしょう。
拡大すると隙間が出来ますし、縮小すると大きな円盤に広がっている星が
小さく縮小される分明るさが積算されて、暗い星との差が大きくなります。

元の画像に近い形で表示したい場合、視野角、中心点の赤経赤緯、回転角などは
写真によってそれぞれ違いますが、2)においてphtlnr.exeが出力した
ファイルの中にそれぞれの画像についての値が記載されているので、
それを参考に適当に決めてください。

入力画像をjpgにした場合は、作業用に(同名).bmpを生成します。
入力画像は、枚数が少ない時はnoisec.exeで固定ノイズを少し低めの除去率で取り除き、
そのbmpを張り合わせた方が良いでしょう。背景が明るい時は、オリジナルのjpgを
そのまま重ね合わせてもいい感じに出来上がります。

===============================================
Reference

hayami.exeのサブルーチンや補助ファイルつかいまくりなのでそっち見てください。
あとは最小自乗法使ってるだけです。

===============================================
2008/01/27:公開
2008/02/24:sumimg.exe の明るさ調整機能拡大
2008/02/28:DLL添付し忘れてたので添付
2008/05/16:正距方位図法以外に対応
===============================================
May-16th,2008
ftcenter@mth.biglobe.ne.jp
http://www2s.biglobe.ne.jp/~ftcenter
