拥有与坐标轴平行边的多边形内的点

Question

我在interviewstreet上解决了多边形问题，但感觉速度太慢了。这个问题的最佳解决方案是什么呢？

在X-Y平面上有N个点，它们的坐标都是整数（xi, yi）。你会得到一组多边形，这些多边形的所有边都与坐标轴平行（换句话说，多边形的所有角都是90度，所有线段都是水平或垂直的，没有对角线）。对于每个多边形，你的程序需要找出有多少个点位于它内部（如果一个点在多边形的边界上，也算在多边形内部）。

输入：

第一行有两个整数N和Q。接下来一行包含N个用空格分开的整数坐标（xi, yi）。然后是Q个查询。每个查询的第一行是一个整数Mi，接下来是Mi个用空格分开的整数坐标（x[i][j], y[i][j]），这些坐标按顺时针顺序指定查询多边形的边界。

多边形是由垂直线段和水平线段交替组成的。多边形有Mi条边，其中(x[i][j], y[i][j])与(x[i][(j+1)%Mi], y[i][(j+1)%Mi])相连。对于每个0 <= j < Mi，要么x[i][(j+1)%Mi] == x[i][j]，要么y[i][(j+1)%Mi] == y[i][j]，但不能同时满足。还保证多边形不会自交。

输出：

对于每个查询，输出该查询多边形内部的点的数量，每个结果占一行。

示例输入 #1：

16 2
0 0
0 1
0 2
0 3
1 0
1 1
1 2
1 3
2 0
2 1
2 2
2 3
3 0
3 1
3 2
3 3
8
0 0
0 1
1 1
1 2
0 2
0 3
3 3
3 0
4
0 0
0 1
1 1
1 0

示例输出 #1：

16
4

示例输入 #2：

6 1
1 1
3 3
3 5
5 2
6 3
7 4
10
1 3
1 6
4 6
4 3
6 3
6 1
4 1
4 2
3 2
3 3

示例输出 #2：

4

约束条件：

1 <= N <= 20,000
1 <= Q <= 20,000
4 <= Mi <= 20
每个坐标的值最多为200,000

我对提到的语言或伪代码的解决方案感兴趣。

编辑：这是我的代码，但它的复杂度是O(n^2)

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;

namespace Polygon
{
// avoding System.Drawing dependency
public struct Point
{
    public int X { get; private set; }
    public int Y { get; private set; }

    public Point(int x, int y)
        : this()
    {
        X = x;
        Y = y;
    }

    public override int GetHashCode()
    {
        return X ^ Y;
    }

    public override bool Equals(Object obj)
    {
        return obj is Point && this == (Point)obj;
    }

    public static bool operator ==(Point a, Point b)
    {
        return a.X == b.X && a.Y == b.Y;
    }

    public static bool operator !=(Point a, Point b)
    {
        return !(a == b);
    }
}

public class Solution
{
    static void Main(string[] args)
    {
        BasicTestCase();
        CustomTestCase();
        // to read from STDIN
        //string firstParamsLine = Console.ReadLine();
        //var separator = new char[] { ' ' };
        //var firstParams = firstParamsLine.Split(separator);
        //int N = int.Parse(firstParams[0]);
        //int Q = int.Parse(firstParams[1]);

        //List<Point> points = new List<Point>(N);

        //for (int i = 0; i < N; i++)
        //{
        //    var coordinates = Console.ReadLine().Split(separator);
        //    points.Add(new Point(int.Parse(coordinates[0]), int.Parse(coordinates[1])));
        //}

        //var polygons = new List<List<Point>>(Q); // to reduce realocation
        //for (int i = 0; i < Q; i++)
        //{

        //    var firstQ = Console.ReadLine().Split(separator);
        //    int coordinatesLength = int.Parse(firstQ[0]);
        //    var polygon = new List<Point>(coordinatesLength);
        //    for (int j = 0; j < coordinatesLength; j++)
        //    {
        //        var coordinates = Console.ReadLine().Split(separator);
        //        polygon.Add(new Point(int.Parse(coordinates[0]), int.Parse(coordinates[1])));
        //    }
        //    polygons.Add(polygon);
        //}

        //foreach (var polygon in polygons)
        //{
        //    Console.WriteLine(CountPointsInPolygon(points, polygon));
        //}


    }

    private static void BasicTestCase()
    {
        List<Point> points = new List<Point>(){ new Point(0, 0),
                                                new Point(0, 1),
                                                new Point(0, 2),
                                                new Point(0, 3),
                                                new Point(1, 0),
                                                new Point(1, 1),
                                                new Point(1, 2),
                                                new Point(1, 3),
                                                new Point(2, 0),
                                                new Point(2, 1),
                                                new Point(2, 2),
                                                new Point(2, 3),
                                                new Point(3, 0),
                                                new Point(3, 1),
                                                new Point(3, 2),
                                                new Point(3, 3) };

        List<Point> polygon1 = new List<Point>(){   new Point(0, 0),
                                                    new Point(0, 1),
                                                    new Point(2, 1),
                                                    new Point(2, 2),
                                                    new Point(0, 2),
                                                    new Point(0, 3),
                                                    new Point(3, 3),
                                                    new Point(3, 0)};

        List<Point> polygon2 = new List<Point>(){   new Point(0, 0),
                                                    new Point(0, 1),
                                                    new Point(1, 1),
                                                    new Point(1, 0),};

        Console.WriteLine(CountPointsInPolygon(points, polygon1));
        Console.WriteLine(CountPointsInPolygon(points, polygon2));

        List<Point> points2 = new List<Point>(){new Point(1, 1),
                                                new Point(3, 3),
                                                new Point(3, 5),
                                                new Point(5, 2),
                                                new Point(6, 3),
                                                new Point(7, 4),};

        List<Point> polygon3 = new List<Point>(){   new Point(1, 3),
                                                    new Point(1, 6),
                                                    new Point(4, 6),
                                                    new Point(4, 3),
                                                    new Point(6, 3),

                                                    new Point(6, 1),
                                                    new Point(4, 1),
                                                    new Point(4, 2),
                                                    new Point(3, 2),
                                                    new Point(3, 3),};

        Console.WriteLine(CountPointsInPolygon(points2, polygon3));
    }

    private static void CustomTestCase()
    {
        // generated 20 000 points and polygons
        using (StreamReader file = new StreamReader(@"in3.txt"))
        {
            string firstParamsLine = file.ReadLine();
            var separator = new char[] { ' ' };
            var firstParams = firstParamsLine.Split(separator);
            int N = int.Parse(firstParams[0]);
            int Q = int.Parse(firstParams[1]);

            List<Point> pointsFromFile = new List<Point>(N);

            for (int i = 0; i < N; i++)
            {
                var coordinates = file.ReadLine().Split(separator);
                pointsFromFile.Add(new Point(int.Parse(coordinates[0]), int.Parse(coordinates[1])));
            }

            var polygons = new List<List<Point>>(Q); // to reduce realocation
            for (int i = 0; i < Q; i++)
            {

                var firstQ = file.ReadLine().Split(separator);
                int coordinatesLength = int.Parse(firstQ[0]);
                var polygon = new List<Point>(coordinatesLength);
                for (int j = 0; j < coordinatesLength; j++)
                {
                    var coordinates = file.ReadLine().Split(separator);
                    polygon.Add(new Point(int.Parse(coordinates[0]), int.Parse(coordinates[1])));
                }
                polygons.Add(polygon);
            }

            foreach (var polygon in polygons)
            {
                Console.WriteLine(CountPointsInPolygon(pointsFromFile, polygon));
            }

        }
    }

    public static int CountPointsInPolygon(List<Point> points, List<Point> polygon)
    {
        // TODO input check
        polygon.Add(polygon[0]); // for simlicity
        // check if any point is outside of the bounding box of the polygon
        var minXpolygon = polygon.Min(p => p.X);
        var maxXpolygon = polygon.Max(p => p.X);
        var minYpolygon = polygon.Min(p => p.Y);
        var maxYpolygon = polygon.Max(p => p.Y);
        // ray casting algorithm (form max X moving to point)
        int insidePolygon = 0;
        foreach (var point in points)
        {
            if (point.X >= minXpolygon && point.X <= maxXpolygon && point.Y >= minYpolygon && point.Y <= maxYpolygon)
            {       // now points are inside the bounding box 
                    isPointsInside(polygon, point, ref insidePolygon);
            } // else outside
        }
        return insidePolygon;

    }
    private static void isPointsInside(List<Point> polygon, Point point, ref int insidePolygon)
    {
        int intersections = 0;    

        for (int i = 0; i < polygon.Count - 1; i++)
        {
            if (polygon[i] == point)
            {
                insidePolygon++;
                return;
            }
            if (point.isOnEdge(polygon[i], polygon[i + 1]))
            {
                insidePolygon++;
                return;
            }
            if (Helper.areIntersecting(polygon[i], polygon[i + 1], point))
            {
                intersections++;
            }
        }

        if (intersections % 2 != 0)
        {
            insidePolygon++;
        }   
    }
}

static class Helper
{
    public static bool isOnEdge(this Point point, Point first, Point next)
    {
            // onVertical 
            if (point.X == first.X && point.X == next.X && point.Y.InRange(first.Y, next.Y))
            {
                return true;
            }
            //onHorizontal 
            if (point.Y == first.Y && point.Y == next.Y && point.X.InRange(first.X, next.X))
            {
                return true;
            }
            return false;
    }

    public static bool InRange(this int value, int first, int second)
    {
        if (first <= second)
        {
            return value >= first && value <= second;
        }
        else
        {
            return value >= second && value <= first;
        }
    }

    public static bool areIntersecting(Point polygonPoint1, Point polygonPoint2, Point vector2End)
    {
        // "move" ray up for 0.5 to avoid problem with parallel edges
        if (vector2End.X < polygonPoint1.X )
        {
            var y = (vector2End.Y + 0.5);
            var first = polygonPoint1.Y;
            var second = polygonPoint2.Y;
            if (first <= second)
            {
                return y >= first && y <= second;
            }
            else
            {
                return y >= second && y <= first;
            }
        }
        return false;
    }
}

}

Answer 1

梯形分解对你有用吗？

Answer 2

一个更快的解决办法是把这些点放进一个四叉树里。

区域四叉树可能更容易编写代码，但点四叉树可能会更快。如果使用区域四叉树，当某个区域里的点数少于一个阈值（比如16个点）时，可以停止继续细分这个四叉树。

每个区域会存储它包含的点的数量，还有一个坐标列表（对于叶子节点）或者指向更小区域的指针。（当区域的大小达到1时，可以省略坐标列表，因为这些点一定是重合的）

要计算多边形内部的点数，你需要查看代表四叉树根节点的最大区域。

1. 把多边形裁剪到区域的边缘
2. 如果多边形和区域没有重叠，就返回0
3. 如果这个区域的大小是1x1，就返回这个区域里的点数
4. 如果多边形完全包围了这个区域，就返回这个区域里的点数。
5. 如果这是一个叶子节点，就用铅垂线算法测试每个点
6. 否则，递归地计算每个子区域里的点数

（如果一个区域只有一个非空的子区域，你可以跳过第1、2、3、4、5步，这样会快一点）

（第2和第4步的测试不需要完全准确）

Answer 3

这是作者提供的解决方案 - 有点复杂，不是吗？

#include <iostream>
#include <ctime>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <algorithm>
#include <vector>
#include <ctime>
using namespace std;

typedef long long int64;
const int N = 100000, X = 2000000001;
const int Q = 100000, PQ = 20;
struct Point {
    int x, y, idx;
    Point(int _x = 0, int _y = 0, int _idx = 0) {
        x = _x;
        y = _y;
        idx = _idx;
    }
} arr_x[N], arr_y[N];
struct VLineSegment {
    int x, y1, y2, idx, sign;
    VLineSegment(int _x = 0, int _y1 = 0, int _y2 = 0, int _sign = 1, int _idx = 0) {
        x = _x;
        y1 = _y1;
        y2 = _y2;
        sign = _sign;
        idx = _idx;
    }
    bool operator<(const VLineSegment& v) const {
        return x < v.x;
    }
} segs[Q * PQ];

struct TreeNode {
    int idx1, idx2, cnt;
    TreeNode *left, *right;
    TreeNode() { left = right = 0; cnt = 0; }
    ~TreeNode() { if(left) delete left; if(right) delete right; }

    void update_stat() {
        cnt = left->cnt + right->cnt;
    }

    void build(Point* arr, int from, int to, bool empty) {
        idx1 = from;
        idx2 = to;
        if(from == to) {
            if(!empty) {
                cnt = 1;
            } else { 
                cnt = 0;
            }
        } else {
            left = new TreeNode();
            right = new TreeNode();
            int mid = (from + to) / 2;
            left->build(arr, from, mid, empty);
            right->build(arr, mid + 1, to, empty);
            update_stat();
        }
    }
    void update(Point& p, bool add) {
        if(p.idx >= idx1 && p.idx <= idx2) {
            if(idx1 != idx2) {
                left->update(p, add);
                right->update(p, add);
                update_stat();
            } else {
                if(add) {
                    cnt = 1;
                } else {
                    cnt = 0;
                }
            }
        }
    }
    int query(int ya, int yb) {
        int y1 = arr_y[idx1].y, y2 = arr_y[idx2].y;
        if(ya <= y1 && y2 <= yb) {
            return cnt;
        } else if(max(ya, y1) <= min(yb, y2)) {
            return left->query(ya, yb) + right->query(ya, yb);
        }
        return 0;
    }
};

bool cmp_x(const Point& a, const Point& b) {
    return a.x < b.x;
}
bool cmp_y(const Point& a, const Point& b) {
    return a.y < b.y;
}

void calc_ys(int x1, int y1, int x2, int y2, int x3, int sign, int& ya, int& yb) {
    if(x2 < x3) {
        yb = 2 * y2 - sign;
    } else {
        yb = 2 * y2 + sign;
    }
    if(x2 < x1) {
        ya = 2 * y1 + sign;
    } else {
        ya = 2 * y1 - sign;
    }
}

bool process_polygon(int* x, int* y, int cnt, int &idx, int i) {
    for(int j = 0; j < cnt; j ++) {
        //cerr << x[(j + 1) % cnt] - x[j] << "," << y[(j + 1) % cnt] - y[j] << endl;
        if(x[j] == x[(j + 1) % cnt]) {
            int _x, y1, y2, sign;
            if(y[j] < y[(j + 1) % cnt]) {
                _x = x[j] * 2 - 1;
                sign = -1;
                calc_ys(x[(j + cnt - 1) % cnt], y[j], x[j], y[(j + 1) % cnt], x[(j + 2) % cnt], sign, y1, y2);
            } else {
                _x = x[j] * 2 + 1;
                sign = 1;
                calc_ys(x[(j + 2) % cnt], y[(j + 2) % cnt], x[j], y[j], x[(j + cnt - 1) % cnt], sign, y1, y2);
            }
            segs[idx++] = VLineSegment(_x, y1, y2, sign, i);
        }
    }
}
int results[Q];

int n, q, c;

int main() {
    int cl = clock();
    cin >> n >> q;
    for(int i = 0; i < n; i ++) {
        cin >> arr_y[i].x >> arr_y[i].y;
        arr_y[i].x *= 2;
        arr_y[i].y *= 2;
    }
    int idx = 0, cnt, x[PQ], y[PQ];
    for(int i = 0; i < q; i ++) {
        cin >> cnt;
        for(int j = 0; j < cnt; j ++) cin >> x[j] >> y[j];
        process_polygon(x, y, cnt, idx, i);
    }
    sort(segs, segs + idx);
    memset(results, 0, sizeof results);
    sort(arr_y, arr_y + n, cmp_y);
    for(int i = 0; i < n; i ++) {
        arr_y[i].idx = i;
        arr_x[i] = arr_y[i];
    }
    sort(arr_x, arr_x + n, cmp_x);
    TreeNode tleft;
    tleft.build(arr_y, 0, n - 1, true);
    for(int i = 0, j = 0; i < idx; i ++) {
        for(; j < n && arr_x[j].x <= segs[i].x; j ++) {
            tleft.update(arr_x[j], true);
        }
        int qcnt = tleft.query(segs[i].y1, segs[i].y2);
        //cerr << segs[i].x * 0.5 << ", " << segs[i].y1 * 0.5 << ", " << segs[i].y2 * 0.5 << " = " << qcnt << " * " << segs[i].sign << endl;
        results[segs[i].idx] += qcnt * segs[i].sign;
    }
    for(int i = 0; i < q; i ++) {
        cout << results[i] << endl;
    }
    cerr << (clock() - cl) * 0.001 << endl;
    return 0;
}