領域の内側か外側かを調べる

はじめに

複素関数の定理で閉領域内なら2πi、領域外なら0を返すものがあります。この定理を使うことで指定した点が領域内にあるかどうかを判定できます。参考書は「複素関数攻略の一本道,森北出版」です。

2πiの定理

コーシーの積分定理はwikipediaによると、以下の通りです。

D を単連結領域とし、f(z) は D 上で正則である複素関数とするとき、C を D 内にある長さを持つ閉曲線とすると、

この定理から「2πiの定理」が成り立つことが分かります。

z0がC外なら関数はC内で正則になるので、積分値はコーシーの積分定理により０となります。
一方、z0がC内なら
とおくと、

より、積分値は2πiになります。
つまり「2πiの定理」を評価すれば点が閉領域内なのか領域外なのかを調べることができます。

実装

実装をするためには積分の計算をする必要があります。積分の実装は難しいので、参考書では"G.S. Gipson, Use of the residue theorem in locating points within an arbitrary multiply-connected region, 1986."で提案されている式変換を使っています。
これは閉曲線CをN個のセグメントに分けて計算する手法です。

より、

が成り立ちます。

θはz-z0が実軸となす角です。r1をz0を始点として、セグメントの開始位置までのベクトルとすると、
より

となります。また、閉曲線の進む向きは外積で求めることができます。

実行例

適当に決めた点の内側にカーソルがあるかを判定します。マウスのMoveイベントで毎回評価しています。

ソースコード

private enum Result
{
    IN,
    OUT,
};

List<PointF> _points = new List<PointF>
        {
            new PointF(100f, 49.0f),
            new PointF(59f, 77.0f),
            new PointF(57.0f, 178.0f),
            new PointF(122.0f, 207.0f),
            new PointF(191.0f, 167.0f),
            new PointF(197.0f, 105.0f),
            new PointF(135.0f, 77.0f),
        };

public Form1()
{
    InitializeComponent();

    var bitmap = new Bitmap(pictureBox1.Size.Width, pictureBox1.Size.Height);
    pictureBox1.Image = bitmap;

    // 点と線を表示
    var g = Graphics.FromImage(bitmap);
    g.DrawPolygon(Pens.Blue, _points.Select(v => new Point((int)v.X, (int)v.Y)).ToArray());
    foreach (var p in _points)
    {
        //g.FillRectangle(Brushes.Black, new Rectangle((int)p.X, (int)p.Y, 10, 10));
        g.FillEllipse(Brushes.Black, new Rectangle((int)p.X-5, (int)p.Y-5, 10, 10));
    }
}

private Result InOrOut(List<PointF> points, PointF z0)
{
    // 内積
    Func<PointF, PointF, double> dp 
        = (r1, r2) => r1.X * r2.X + r1.Y * r2.Y;

    // ノルム
    Func<PointF, double> norm 
        = r => Math.Sqrt(r.X * r.X + r.Y * r.Y);

    // 一個ずらした要素をまとめる
    var pack = points.Zip(points.Skip(1), (z1, z2) => new {Z1=z1, Z2=z2}).ToList();
    pack.Add(new { Z1=points.Last(), Z2=points.First() });

    double integral = 0.0;
    foreach (var p in pack)
    {
        var r1 = new PointF(p.Z1.X - z0.X, p.Z1.Y - z0.Y);
        var r2 = new PointF(p.Z2.X - z0.X, p.Z2.Y - z0.Y);

        // 反時計回り+ or 時計回り-
        int sign = Math.Sign(r1.X*r2.Y - r1.Y*r2.X);

        integral += sign * Math.Acos(dp(r1, r2) / (norm(r1) * norm(r2)));
    };

    // 積分結果がだいたい2πiなら内側
    const double epsilon = 0.1;
    return Math.Abs(Math.Abs(integral / (2.0 * Math.PI)) - 1.0) < epsilon ? Result.IN : Result.OUT;
}

private void pictureBox1_MouseMove(object sender, MouseEventArgs e)
{
    var r = InOrOut(_points, new PointF(e.X, e.Y));
    this.label1.Text = string.Format("({0}, {1}) : {2}", e.X, e.Y, r);
}

複素関数の定理がプログラムに使えるのはおもしろいですね。関係ないですが「複素関数攻略の一本道」は途中式を省略せずに書いてあるので、数学の苦手は私でも分かった気になります。そういう本がふえてくれないかなー。