天体画像で遊ぶ (2) 大きな位置誤差を持つ2枚の画像の位置合わせ

　[前回](https://qiita.com/amutou/items/cc6d010c4df8a7637f8b)の続き

***

　以下に今回使用する2枚の画像を示す。

![DSC01336.JPG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264452/d544e587-b9dc-153d-f9ad-c176f15eb563.jpeg)

![DSC01337.JPG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264452/aa7728b7-5f97-e1db-6134-b5d579e8ecb7.jpeg)

　重ねた画像を以下に示す。

![log1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264452/784be370-5c45-0240-29d0-3cf96fcce473.png)

　天体(輝点)の分布に比べて大きなズレを持っている。

***

　大きな位置誤差の推定には[OPMアルゴリズム](http://stella.astron.s.u-tokyo.ac.jp/nmatsuna/Japanese/software/OPM.html)と呼ばれる手法を参考にした。ただし元の論文やソースコードは参照していないので、OPMアルゴリズムそれ自体とは異なる実装である。

　流れとしては

1) 輝点の抽出
1) 輝点を3個選び、3角形を作る
1) 三角形の短辺を挟む辺の内積と、長辺と短辺の長さの比の2つの値を三角形の比較に使用する
1) 2枚の画像から適当な数の三角形を作成し、三角形を総当りで比較する
1) 同じ形状の三角形が見つかった場合は、それらを変換する行列を作成し、すべての輝点に対してペアリングを行う
1) ペアリングして十分な一致が見られた場合はその行列を採用する
1) 十分な一致が得られない場合は別の三角形の組み合わせを試す

***

　前回の差分のコードを以下に示す。

```C#
static void Main()
{
    var labelingThreshold = 200;
    var dotMax = 1.1f;
    var ratioMax = 1.01f;
    var pairingTolerance = 15;

    var blue1 = LoadImage("DSC01336.JPG").Colors().Select(a => (int)a.B);
    var blue2 = LoadImage("DSC01337.JPG").Colors().Select(a => (int)a.B);

    var labels1 = Labeling(blue1, labelingThreshold, out _)
        .Skip(1).Where(a => 4 <= a.Area).OrderByDescending(a => a.Weight).Select(a => a.CenterOfGravity).ToArray();
    var labels2 = Labeling(blue2, labelingThreshold, out _)
        .Skip(1).Where(a => 4 <= a.Area).OrderByDescending(a => a.Weight).Select(a => a.CenterOfGravity).ToArray();

    var triangles1 = Triangle.Create(labels1, new Triangle[100]);
    var triangles2 = Triangle.Create(labels2, new Triangle[100]);

    var bmp1 = blue1.Select(a => Color.FromArgb(a, 0, 0)).ToBitmap();
    var bmp2 = blue2.Select(a => Color.FromArgb(0, a, a)).ToBitmap();

    var dst = BitmapAddSatulation(bmp1, bmp2);
    var g = Graphics.FromImage(dst);

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

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

    foreach (var (i, j, mtx) in Triangle.Search(triangles1, triangles2, dotMax, ratioMax))
    {
        var pairs = Pairing(
            labels1,
            labels2.Select(a => Vector2.Transform(a, mtx)).ToArray(),
            pairingTolerance, out var pairingRatio);

        var t1 = triangles1[i];
        var t2 = triangles2[j];
        Console.WriteLine($"{i} {j} {pairingRatio:0.0000} {t1.Dot / t2.Dot:0.000} {t1.Ratio / t2.Ratio:0.000}");

        if (pairingRatio < 0.5) { continue; }

        var matrix = mtx * ToMatrix(CollectionsMarshal.AsSpan(pairs));

        Pen pen = new(Color.FromArgb(96, Color.White));
        foreach (var (a, b) in labels2.Select(a => (a, Vector2.Transform(a, matrix))))
        {
            g.DrawEllipse(Pens.White, b.X - 8, b.Y - 8, 16, 16);
            var c = Vector2.Normalize(b - a);
            g.DrawLine(pen, (PointF)(a + c * 10), new(b - c * 8));
        }

        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));

        BitmapAddSatulation(bmp1, bmp2.Transform(matrix))
            .Save("log2.png");

        break;
    }

    dst.Save("log1.png");
}

record Triangle
{
    public float Ratio { get; private init; }
    public float Dot { get; private init; }
    public Vector2 B { get; private init; }
    public Vector2 BA { get; private init; }

    public static Triangle Create(Vector2 a, Vector2 b, Vector2 c)
    {
        var ab = Vector2.DistanceSquared(a, b);
        var bc = Vector2.DistanceSquared(b, c);
        var ca = Vector2.DistanceSquared(c, a);

        (a, b, c, ab, _, ca) =
            ab >= bc && bc >= ca ? (a, b, c, ab, bc, ca) :
            ca >= bc && bc >= ab ? (a, c, b, ca, bc, ab) :
            ab >= ca && ca >= bc ? (b, a, c, ab, ca, bc) :
            bc >= ca && ca >= ab ? (b, c, a, bc, ca, ab) :
            ca >= ab && ab >= bc ? (c, a, b, ca, ab, bc) :
            bc >= ab && ab >= ca ? (c, b, a, bc, ab, ca) :
            throw new Exception();

        return new()
        {
            Ratio = float.Sqrt(ca / ab),
            Dot = Vector2.Dot(a - b, c - b),
            B = b,
            BA = a - b
        };
    }

    public static ReadOnlySpan<Triangle> Create(ReadOnlySpan<Vector2> src, Span<Triangle> dst)
    {
        int n = 0;
        for (int i = 2; i < src.Length; i++)
        {
            for (int j = 1; j < i; j++)
            {
                for (int k = 0; k < j; k++)
                {
                    dst[n++] = Create(src[i], src[j], src[k]);
                    if (dst.Length <= n) { goto exit; }
                }
            }
        }
    exit:;
        return dst[..n];
    }

    public static IEnumerable<(int i, int j, Matrix3x2 mtx)> Search(
        ReadOnlySpan<Triangle> triangles1, ReadOnlySpan<Triangle> triangles2, float dotMax, float ratioMax)
        => Search(triangles1.ToArray(), triangles2.ToArray(), dotMax, ratioMax);

    public static IEnumerable<(int i, int j, Matrix3x2 mtx)> Search(
        Triangle[] triangles1, Triangle[] triangles2, float dotMax, float ratioMax)
    {
        if (dotMax <= 1) { throw new ArgumentException(nameof(dotMax)); }
        if (ratioMax <= 1) { throw new ArgumentException(nameof(ratioMax)); }

        var dotMin = 1 / dotMax;
        var ratioMin = 1 / ratioMax;
        var n = int.Min(triangles1.Length, triangles2.Length);
        for (int k = 0; k < n; k++)
        {
            for (int i = 0; i <= k; i++)
            {
                var j = k - i;
                var t1 = triangles1[i];
                var t2 = triangles2[j];

                var dot = t1.Dot / t2.Dot;
                if (dot < dotMin || dotMax < dot) { continue; }

                var ratio = t1.Ratio / t2.Ratio;
                if (ratio < ratioMin || ratioMax < ratio) { continue; }

                var (a, b) = (t1.BA, t2.BA);
                var rad = float.Atan2(
                    a.Y * b.X - a.X * b.Y,
                    a.X * b.X + a.Y * b.Y);

                yield return (i, j,
                    Matrix3x2.CreateTranslation(-t2.B) *
                    Matrix3x2.CreateRotation(rad) *
                    Matrix3x2.CreateTranslation(t1.B));
            }
        }
    }
}
```

　三角形を使用した位置補正では、理想的な状況ではかなり正確だが、実際の画像を使うとレンズ歪みなどによって正確なパラメータが得られないことがある。その場合でもある程度の精度ではパラメータが得られているので、前回の微小な位置誤差の補正を使用してパラメータの細かい修正を行っている。

　このようにして得られた画像を以下に示す。

![](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264452/559baf40-a318-acc1-008c-8e13b36b3c7a.png)

![](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264452/a701c331-1b74-8c78-9f4f-587c0b6da9bb.png)

　ある程度正しく推定できているが、レンズ歪みの影響によって(特に周辺部で)正しく重ならない。