　今回から何回かのシリーズで天体画像を使用したプログラミングの遊びを行う。プログラミング環境はC#(.NET 7)を使用しており、一部の文法は古いC#では対応していない。また、環境依存のコードが含まれている可能性があり、すべての環境で動作することを保証しない。

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

　ほとんど同じ画像だが、微妙に位置が違う。どれくらい位置がずれているかを検出する。

　2枚の画像の青色を取り出し、1枚を赤色で、もう1枚を水色(シアン)で表現し、重ねると以下のようになる。

　位置がずれているので、それぞれの赤色と水色の輝点が少し離れた場所に表示される。

　手順としては

ラベリング(輝点の抽出)
輝点のペアリング
ペアからズレを求める

　というふうに行い、最終的にずれを補正する。
　今回はズレが微小だと仮定して処理を行う。ここで言う微小とは、輝点の密度に比べてずれが十分に小さい状況とし、2つの画像の一番近い輝点をペアとして採用する。

　ペアリングの処理は以下の様になる。

static List<(Vector2 a, Vector2 b)> Pairing(ReadOnlySpan<Vector2> a, ReadOnlySpan<Vector2> b, float tolerance, out float ratio)
{
    if (a.Length > b.Length) { a = a[..b.Length]; }
    if (a.Length < b.Length) { b = b[..a.Length]; }

    var sq = tolerance * tolerance;
    List<(Vector2, Vector2)> result = new(a.Length);
    var bUsed = new bool[b.Length];

    for (int i = 0; i < a.Length; i++)
    {
        for (int j = 0; j < b.Length; j++)
        {
            if (bUsed[j] ||
                sq < Vector2.DistanceSquared(a[i], b[j]))
            { continue; }

            result.Add((a[i], b[j]));
            bUsed[j] = true;
            break;
        }
    }

    ratio = (float)result.Count / a.Length;
    return result;
}

　System.Numerics.Vector2で与えられた2組の輝点位置を総当りでペアリングしていく。
　極端に数が違う輝点の組を与えられた場合、不正な組み合わせでペアリングしてしまうことがあるため、最初に先頭の同じ数だけを使うようしている。このため、輝点はあらかじめ輝度でソートしておき、明るい輝点を優先的に使うようにする。

　ペアが決まれば、続いて位置の差から移動と回転を求める。

static Matrix3x2 ToMatrix(ReadOnlySpan<(Vector2, Vector2)> pairs)
{
    double x = 0, y = 0;
    double u = 0, v = 0;
    foreach (var (a, b) in pairs)
    {
        var c = a + b;
        x += c.X;
        y += c.Y;

        var d = b - a;
        u += d.X;
        v += d.Y;
    }
    Vector2 center = new(
        (float)(x / pairs.Length / 2),
        (float)(y / pairs.Length / 2));
    Vector2 move = new(
        (float)(u / pairs.Length),
        (float)(v / pairs.Length));

    x = y = 0;
    foreach (var (a, b) in pairs)
    {
        var c = a - center;
        var d = b - center + move;
        x += c.X * d.X + c.Y * d.Y;
        y += c.Y * d.X - c.X * d.Y;
    }

    return
        Matrix3x2.CreateTranslation(-center) *
        Matrix3x2.CreateRotation((float)double.Atan2(y, x)) *
        Matrix3x2.CreateTranslation(center - move);
}

　もっと良いやり方があると思うが、とりあえずこの程度でもうまく動く。

　ペアを表現する行列が得られるから、この行列を使えば片方の輝点をもう片方の輝点の場所へ移動できる。

　橙・緑・白のマーカーがあるが、1枚目の輝点を橙で、2枚目の輝点を緑で示し、2枚目の輝点の位置を行列で変換した場所を白で示している。また、この画像では見づらいが、白線でペア間を接続している。橙と白のマーカーがほぼ同じ位置に重なっており、正しく変換行列が推定できていることがわかる。

　行列を使用して画像を変形し重ねると以下のようになる。

　赤色と水色の輝点が重なって白色の輝点として見えるようになっている。

　今回使用したプログラムの全文を以下に示す。

　ラベリングした結果は閾値未満の結果が含まれているため、最初の1要素をスキップし、また面積の狭いラベル(今回は4ピクセル未満)はノイズとして除去している。
　今回は手抜きのために、内部リソース(BitmapやGraphics等)の開放処理は省略している。

using System.Buffers;
using System.Drawing;
using System.Drawing.Imaging;
using System.Numerics;
using System.Runtime.InteropServices;

#pragma warning disable CA1416 // Validate platform compatibility

internal class Program
{
    static void Main()
    {
        var labelingThreshold = 200;
        var pairingTolerance = 15;

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

        var labels1 = Labeling(blue1, labelingThreshold, out _);
        var labels2 = Labeling(blue2, labelingThreshold, out _);

        labels1 = labels1.Skip(1).Where(a => 4 <= a.Area).OrderByDescending(a => a.Peak).ToArray();
        labels2 = labels2.Skip(1).Where(a => 4 <= a.Area).OrderByDescending(a => a.Peak).ToArray();

        var pairs = Pairing(
            labels1.Select(a => a.CenterOfGravity).ToArray(),
            labels2.Select(a => a.CenterOfGravity).ToArray(),
            pairingTolerance, out var pairRatio);
        Console.WriteLine(pairRatio);

        if (pairRatio < 0.7) { throw new Exception(); }

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

        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.Select(a => a.CenterOfGravity))
        { g.DrawEllipse(Pens.Orange, a.X - 10, a.Y - 10, 20, 20); }

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

        foreach (var (a, b) in pairs)
        { g.DrawLine(Pens.White, (PointF)a, new(b)); }

        foreach (var a in labels2.Select(a => Vector2.Transform(a.CenterOfGravity, mtx)))
        { g.DrawEllipse(Pens.White, a.X - 8, a.Y - 8, 16, 16); }

        dst.Save("log1.png");

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

    static Matrix3x2 ToMatrix(ReadOnlySpan<(Vector2, Vector2)> pairs)
    {
        double x = 0, y = 0;
        double u = 0, v = 0;
        foreach (var (a, b) in pairs)
        {
            var c = a + b;
            x += c.X;
            y += c.Y;

            var d = b - a;
            u += d.X;
            v += d.Y;
        }
        Vector2 center = new(
            (float)(x / pairs.Length / 2),
            (float)(y / pairs.Length / 2));
        Vector2 move = new(
            (float)(u / pairs.Length),
            (float)(v / pairs.Length));

        x = y = 0;
        foreach (var (a, b) in pairs)
        {
            var c = a - center;
            var d = b - center + move;
            x += c.X * d.X + c.Y * d.Y;
            y += c.Y * d.X - c.X * d.Y;
        }

        return
            Matrix3x2.CreateTranslation(-center) *
            Matrix3x2.CreateRotation((float)double.Atan2(y, x)) *
            Matrix3x2.CreateTranslation(center - move);
    }

    static List<(Vector2 a, Vector2 b)> Pairing(ReadOnlySpan<Vector2> a, ReadOnlySpan<Vector2> b, float tolerance, out float ratio)
    {
        if (a.Length > b.Length) { a = a[..b.Length]; }
        if (a.Length < b.Length) { b = b[..a.Length]; }

        var sq = tolerance * tolerance;
        List<(Vector2, Vector2)> result = new(a.Length);
        var bUsed = new bool[b.Length];

        for (int i = 0; i < a.Length; i++)
        {
            for (int j = 0; j < b.Length; j++)
            {
                if (bUsed[j] ||
                    sq < Vector2.DistanceSquared(a[i], b[j]))
                { continue; }

                result.Add((a[i], b[j]));
                bUsed[j] = true;
                break;
            }
        }

        ratio = (float)result.Count / a.Length;
        return result;
    }

    static (int Area, Vector2 CenterOfGravity, long Weight, int Peak)[] Labeling(int[][] src, int threshold, out int[][] map)
    {
        if (src.Any(a => a.Length != src[0].Length)) { throw new InvalidOperationException(); }

        var LUT = ArrayPool<int>.Shared.Rent(src[0].Length);
        var iLUT = 0;
        LUT[iLUT] = iLUT++;

        map = new int[src.Length][];
        for (int y = 0; y < src.Length; y++)
        {
            var line = src[y];
            map[y] = new int[line.Length];

            {
                var required = iLUT + line.Length;
                if (LUT.Length < required)
                {
                    var next = ArrayPool<int>.Shared.Rent(required);
                    LUT.CopyTo(next, 0);
                    ArrayPool<int>.Shared.Return(LUT);
                    LUT = next;
                }
            }

            for (int x = 0; x < line.Length; x++)
            {
                if (line[x] < threshold) { continue; }

                var a = x <= 0 ? 0 : map[y][x - 1];
                var b = y <= 0 ? 0 : map[y - 1][x];

                while (a != LUT[a]) { a = LUT[a]; }
                while (b != LUT[b]) { b = LUT[b]; }

                map[y][x] =
                    a == 0 && b == 0 ? LUT[iLUT] = iLUT++ :
                    a == 0 ? b :
                    b == 0 ? a :
                    a < b ?
                        LUT[b] = a :
                        LUT[a] = b;
            }
        }

        int labelN = 0;
        for (int i = 0; i < LUT.Length && i < iLUT; i++)
        {
            LUT[i] = LUT[i] == i ? labelN++ : LUT[LUT[i]];
        }

        var result = new (int area, Vector2 cg, long weight, int peak)[labelN];
        var xyAccum = new (long x, long y)[labelN];

        for (int y = 0; y < map.Length; y++)
        {
            var line = map[y];
            for (int x = 0; x < line.Length; x++)
            {
                var i = line[x] = LUT[line[x]];
                ref var a = ref result[i];
                ref var b = ref xyAccum[i];
                var weight = src[y][x];

                a.area++;
                a.weight += weight;
                if (a.peak < weight) { a.peak = weight; }

                b.x += x * weight;
                b.y += y * weight;
            }
        }

        ArrayPool<int>.Shared.Return(LUT);

        for (int i = 0; i < result.Length; i++)
        {
            var (x, y) = xyAccum[i];
            double w = result[i].weight;
            result[i].cg = new((float)(x / w), (float)(y / w));
        }

        return result;
    }

    static Bitmap LoadImage(string path)
    {
        using FileStream fs = new(path, FileMode.Open, FileAccess.Read);
        using var img = Image.FromStream(fs);
        return new(img);
    }

    static Bitmap BitmapAddSatulation(Bitmap a, Bitmap b)
    {
        if (a.Width != b.Width ||
            a.Height != b.Height) { throw new InvalidOperationException(); }

        Bitmap c = new(a.Width, a.Height);

        var u = a.LockBits(new Rectangle { Size = a.Size }, ImageLockMode.ReadOnly, PixelFormat.Format32bppArgb);
        var v = b.LockBits(new Rectangle { Size = b.Size }, ImageLockMode.ReadOnly, PixelFormat.Format32bppArgb);
        var w = c.LockBits(new Rectangle { Size = c.Size }, ImageLockMode.WriteOnly, PixelFormat.Format32bppArgb);

        var foo = new byte[u.Width * 4];
        var bar = new byte[foo.Length];
        var tmp = new byte[foo.Length];

        for (int y = 0; y < u.Height; y++)
        {
            Marshal.Copy(u.Scan0 + u.Stride * y, foo, 0, foo.Length);
            Marshal.Copy(v.Scan0 + v.Stride * y, bar, 0, bar.Length);

            for (int x = 0; x < tmp.Length; x++)
            {
                tmp[x] = (byte)int.Max(byte.MinValue, int.Min(byte.MaxValue, foo[x] + bar[x]));
            }

            Marshal.Copy(tmp, 0, w.Scan0 + w.Stride * y, tmp.Length);
        }

        a.UnlockBits(u);
        b.UnlockBits(v);
        c.UnlockBits(w);

        return c;
    }
}

static class MyExt
{
    public static Color[][] Colors(this Bitmap bmp)
    {
        var data = bmp.LockBits(
            new Rectangle { Size = bmp.Size },
            ImageLockMode.ReadOnly,
            PixelFormat.Format32bppArgb);
        var argb = new int[data.Width];
        var dst = new Color[data.Height][];

        for (int y = 0; y < dst.Length; y++)
        {
            Marshal.Copy(data.Scan0 + data.Stride * y, argb, 0, argb.Length);
            var line = dst[y] = new Color[argb.Length];
            for (int x = 0; x < argb.Length; x++)
            {
                line[x] = Color.FromArgb(argb[x]);
            }
        }

        bmp.UnlockBits(data);

        return dst;
    }

    public static Bitmap ToBitmap(this Color[][] src)
    {
        if (src.Any(line => line.Length != src[0].Length))
            throw new InvalidOperationException();

        Bitmap bmp = new(src[0].Length, src.Length);
        var data = bmp.LockBits(
            new Rectangle { Size = bmp.Size },
            ImageLockMode.WriteOnly,
            PixelFormat.Format32bppArgb);
        var argb = new int[data.Width];

        for (int y = 0; y < data.Height; y++)
        {
            var line = src[y];
            for (int x = 0; x < argb.Length; x++)
            {
                argb[x] = line[x].ToArgb();
            }
            Marshal.Copy(argb, 0, data.Scan0 + data.Stride * y, argb.Length);
        }

        bmp.UnlockBits(data);
        return bmp;
    }

    public static U[][] Select<T, U>(this T[][] src, Func<T, U> converter)
     => src.Select(a => a.Select(converter).ToArray()).ToArray();

    public static Bitmap Transform(this Bitmap src, Matrix3x2 matrix)
    {
        Bitmap dst = new(src.Width, src.Height);
        using var g = Graphics.FromImage(dst);
        g.MultiplyTransform(new(matrix));
        g.DrawImageUnscaled(src, 0, 0);
        return dst;
    }
}

#pragma warning restore CA1416 // Validate platform compatibility