天体画像で遊ぶ (4) 天体画像と星表の位置合わせ

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

***

　前々回は視野が重なっている状況での位置合わせを行い、前回は星表を平面に投影した。今回は天体画像と星表の位置合わせを行い、天体画像の姿勢を検出する。

```C#
static void Main()
{
    var labelingThreshold = 200;
    var dotMax = 1.1f;
    var ratioMax = 1.01f;
    var pairingTolerance = 15;
    var fLen = 238.9e-3f;
    var pixelPitch = 16.88e-6f;
    var scale = fLen / pixelPitch;

    var RA = 23.5f;
    var Dec = 71f;
    var atti =
        Quaternion.CreateFromAxisAngle(Vector3.UnitZ, RA / 12.0f * float.Pi) *
        Quaternion.CreateFromAxisAngle(-Vector3.UnitY, Dec / 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);

    var labels =
        Labeling(bmp.Colors().Select(a => (int)a.B), labelingThreshold, out _)
        .Skip(1).
        Where(a => 4 <= a.Area)
        .OrderByDescending(a => a.Weight)
        .Select(a => a.CenterOfGravity - new Vector2(bmp.Width / 2, bmp.Height / 2))
        .Where(a => a.Length() < 750)
        .ToArray();

    var triangles1 = Triangle.Create(labels, new Triangle[100]);
    var triangles2 = Triangle.Create(Projection(vec, atti, scale, rect).ToArray(), new Triangle[100]);

    foreach (var a in labels)
    { g.DrawEllipse(Pens.Orange, a.X - 10, a.Y - 10, 20, 20); }

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

    foreach (var (i, j, mtx) in Triangle.Search(triangles1, triangles2, dotMax, ratioMax))
    {
        var t1 = triangles1[i];
        var t2 = triangles2[j];

        var atti2 = atti * ToAttitude(mtx, scale);
        var points = Projection(vec, atti2, scale, rect).ToArray();

        var pairs = Pairing(labels, points, pairingTolerance, out var pairingRatio).ToArray();

        Console.WriteLine($"{i} {j} {pairingRatio:0.0000} {t1.Dot / t2.Dot:0.000} {t1.Ratio / t2.Ratio:0.000}");
        if (pairingRatio < 0.3) { continue; }

        atti2 *= ToAttitude(ToMatrix(pairs), scale);
        points = Projection(vec, atti2, scale, rect).ToArray();

        foreach (var a in points) { g.DrawEllipse(Pens.White, a.X - 5, a.Y - 5, 11, 11); }

        g.DrawLine(Pens.Orange, (PointF)t1.B, new(t1.B + t1.BA));
        g.DrawLine(Pens.Lime, (PointF)t2.B, new(t2.B + t2.BA));
        g.DrawLine(Pens.White, (PointF)t1.B, new(t2.B));
        g.DrawLine(Pens.White, (PointF)(t1.B + t1.BA), new(t2.B + t2.BA));

        {
            var (ra, dec, roll) = AttitudeToRaDecRoll(atti2);
            Console.WriteLine($"{ra / float.Pi * 12} {dec / float.Pi * 180} {roll / float.Pi * 180}");
        }

        break;
    }

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

static Quaternion ToAttitude(Matrix3x2 matrix, float scale)
 => Quaternion.CreateFromAxisAngle(Vector3.UnitX, float.Atan2(matrix.M21, matrix.M22)) *
    Quaternion.CreateFromAxisAngle(Vector3.UnitY, -float.Atan(matrix.M32 / scale)) *
    Quaternion.CreateFromAxisAngle(Vector3.UnitZ, float.Atan(matrix.M31 / scale));

static (float RA, float Dec, float Roll) AttitudeToRaDecRoll(Quaternion atti)
{
    var a = Vector3.Transform(Vector3.UnitX, atti);
    var ra = float.Atan2(a.Y, a.X);
    var dec = float.Atan2(a.Z, float.Sqrt(a.X * a.X + a.Y * a.Y));

    var b = Vector3.Transform(Vector3.UnitZ,
        Quaternion.CreateFromAxisAngle(Vector3.UnitY, dec) *
        Quaternion.CreateFromAxisAngle(Vector3.UnitZ, -ra) *
        atti);
    var roll = float.Atan2(-b.Y, b.Z);

    if (ra < 0) { ra += float.Pi * 2; }

    return (ra, dec, roll);
}
```

　基本的には前回までに使用した手法を組み合わせているだけで、画像から抽出した輝点と、初期姿勢で作成した特徴点から三角形を作成し、それらの三角形を元に大まかな姿勢の修正値を得て、続いて隣接する輝点からさらに細かい修正を行っている。
　画像周辺部のレンズ歪みの影響を除去するため、画像中央部の輝点のみを使用している。
　今回は細かい修正は1回のみ行ったが、場合によっては誤差を吸収しきれない場合があるため、必要に応じて複数回ループする等の工夫が必要になる。また、誤った推定値を破棄するために、pairingRatioに閾値を設定するなどの処理も必要な場合が多い。

***

![](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264452/1730820e-841b-a5a6-c106-d70136684beb.png)

　橙色のマーカーは画像から抽出した輝点を、緑色のマーカーは初期姿勢での特徴点を、白色のマーカーは姿勢更新後の特徴点を示している。橙色と白色のマーカーがほぼ重なっており、正しく姿勢が推定されていることがわかる。
　この手法では視野が重なっていれば検出できるため、光軸周りに180度近い回転があっても正常に処理できる。

　姿勢は`23.69391 70.93062 152.35083`と出力され、前回使用した`23.694077 70.93123 152.3225`とほぼ同じ値が得られている。