天体画像で遊ぶ (3) 天体画像に星表を組み合わせる

　[初回](https://qiita.com/amutou/items/cc6d010c4df8a7637f8b) [前回](https://qiita.com/amutou/items/3a4a9a15490112050e54)

***

　スミソニアン天文台星表を使用する。英語版Wikipediaの外部リンクからダウンロードできるので、これを使用する。

　[Smithsonian Astrophysical Observatory Star Catalog - Wikipedia](https://en.wikipedia.org/wiki/Smithsonian_Astrophysical_Observatory_Star_Catalog)

　Referencesの`Catalog translated to J2000`のリンクからダウンロードする。ただしChromium系ブラウザからアクセスすると多数の警告を出した上でgzを展開して保存されるので、適当なツール(wgetとか)を使用してダウンロードする。圧縮されたデータは約13.1MBで、参考までに手元のファイルのMD5ハッシュ値は`6af2c48d41013bf7796cf6d55c97f305`だった。また、今回はWindowsのダウンロードフォルダに`heasarc_sao.tdat.gz`ファイルを保存している。

***

　SAOファイルの読み込みは以下のように行う。

```C#
static List<(int Id, double RA, double Dec, float VMag)> LoadSAO()
{
    using var sr =
        new StreamReader(
            new GZipStream(
                new FileStream(
                    Path.Combine(
                        Environment.GetFolderPath(Environment.SpecialFolder.Desktop),
                        "../downloads/heasarc_sao.tdat.gz"),
                    FileMode.Open, FileAccess.Read),
                CompressionMode.Decompress));

    if ("<HEADER>" != sr.ReadLine()) { throw new Exception(); }
    while ("<DATA>" != (sr.ReadLine() ?? throw new Exception())) { }

    List<(int, double, double, float)> dst = new(258997);
    string line;
    while ((line = sr.ReadLine() ?? throw new Exception()) != "<END>")
    {
        var tmp = line.Split('|');
        if (float.TryParse(tmp[13], out var vmag))
        {
            dst.Add((int.Parse(tmp[0][4..]), double.Parse(tmp[1]), double.Parse(tmp[5]), vmag));
        }
    }

    return dst;
}
```

　ファイルの場所を探したりgzの展開で面倒なことになっているが、あとはヘッダを確認して読み捨て、データを1行ずつ拾っているだけ。SAOカタログには視等級が設定されていないデータが含まれるが、これは読み捨てている。
　また、固有運動の情報は使用していない。取り扱うのが面倒という理由もあるし、このファイルに含まれている固有運動の情報が信用しきれないという理由もあるし、たかだか25年程度ではそれほど大きく動かないだろう、そもそもfloatの精度では固有運動を表現しきれないだろう、という願望などによる。

***

　天体のR.A./Dec.はまず単位ベクトルへ変換を行い(後述)、そのベクトルを平面に投影するには以下のように行う。

```C#
static List<Vector2> Projection(ReadOnlySpan<Vector3> src, Quaternion atti, float scale, RectangleF rect)
{
    List<Vector2> dst = new();
    var inv = Quaternion.Inverse(atti);
    foreach (var a in src)
    {
        var b = Vector3.Transform(a, inv);
        if (0 < b.X)
        {
            var c = new Vector2(-b.Y, -b.Z) * (scale / b.X);
            if (rect.Contains(c.X, c.Y))
            {
                dst.Add(c);
            }
        }
    }
    return dst;
}
```

　四元数で視野の姿勢を与え、視野の広さはscaleとrectで与える。

***

　画像化は以下のように行う。

```C#
static void Main()
{
    var fLen = 238.9e-3f;
    var pixelPitch = 16.88e-6f;
    var scale = fLen / pixelPitch;

    var RA = 23.694077f;
    var Dec = 70.93123f;
    var Roll = 152.3225f;
    var atti =
        Quaternion.CreateFromAxisAngle(Vector3.UnitZ, RA / 12.0f * float.Pi) *
        Quaternion.CreateFromAxisAngle(-Vector3.UnitY, Dec / 180.0f * float.Pi) *
        Quaternion.CreateFromAxisAngle(Vector3.UnitX, Roll / 180.0f * float.Pi);

    var vec = LoadSAO().OrderBy(a => a.VMag).Select(a =>
    {
        var RA = (float)(a.RA / 180.0 * double.Pi);
        var Dec = (float)(a.Dec / 180.0 * double.Pi);
        return new Vector3(
            float.Cos(Dec) * float.Cos(RA),
            float.Cos(Dec) * float.Sin(RA),
            float.Sin(Dec));
    }).ToArray();

    var bmp = LoadImage("DSC01337.JPG");
    var g = Graphics.FromImage(bmp);
    g.TranslateTransform(bmp.Width / 2, bmp.Height / 2);

    var rect = new RectangleF(-bmp.Width / 2, -bmp.Height / 2, bmp.Width, bmp.Height);

    foreach (var a in Projection(vec, atti, scale, rect))
    { g.DrawEllipse(Pens.Lime, a.X - 7, a.Y - 7, 15, 15); }

    bmp.Save("log.png");
}
```

　最初に視野のパラメータを与える。焦点距離`fLen`および画素ピッチ`pixelPitch`はメートル単位で与え、座標変換の際はその比`scale`を使用する。
　視野の向きはR.A.(時)、Dec.(度)、Roll(度)の3種類を使用し、順にZ軸、-Y軸、X軸で回転させて四元数を得る。
　SAOカタログはまず等級でソートを行った上でR.A./Dec.から単位ベクトルを得る。
　その後、座標変換を行い、画像に書き込む。

***

　結果、以下のような画像が得られる。

![](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264452/37dc2efa-49ed-0b32-12f1-5f1f41ccf1f8.png)

　周辺部はレンズの歪みの影響でズレるが、中心部は概ね輝点と同じ位置にマーカーが現れている。

　SAOカタログにはおよそ25万個の恒星が含まれているが、35mm換算250mm程度の画角でも意外とスカスカな印象になる。

　なお、今回のシリーズで紹介する予定はないが、例えばGDR3(Gaia Data Release 3)というカタログを使った場合、以下のようになる。

![](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264452/fa0916d0-7b11-cfdb-5bd3-a96bf8964ba5.png)

　オレンジの点がザラッと並んでいるが、この点1個1個が1つの恒星の位置を示しており、1ピクセルに複数の恒星が含まれている場所もある。
　GDR3カタログはデータ量が非常に巨大で(圧縮時に約700GB)、複数の領域に分割して配布されているため、1つのファイル(今回は253016-253253領域)を使用するとこのような複雑な形状になる。このファイル1個にはおよそ100万個の恒星が含まれており、その内で画角内には32万個の恒星が含まれている(つまりこの画像1枚のごく一部の領域だけでSAOカタログの全恒星を超える星が含まれている)。

　12.5等級までの恒星(およそ1200個)を表示すると以下のようになる。

![log.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264452/497d3b74-1dae-e7b8-8ab3-33591d153f20.png)

　画像で見えている恒星はほぼすべてマーカーがついており、この写真はおよそ12.5等級までの星が写っていると考えられる。
　なお、SAOカタログには最大でおよそ12等級までの恒星が含まれているが、実際には9等級あたりから星の数が減り始めるため、SAOカタログは9等級あたりまでを使用するのが無難。