Comprobar si un punto pertenece al polígono convexo en $O(\log N)$

Consideremos el siguiente problema: nos dan un polígono convexo con vértices enteros y muchas consultas. Cada consulta es un punto, para el cual debemos determinar si yace dentro o en el borde del polígono o no. Supongamos que el polígono está ordenado en sentido antihorario. Responderemos cada consulta en $O(\log n)$ de forma online.

Algoritmo

Elijamos el punto con la menor coordenada x. Si hay varios, elegimos el de menor coordenada y. Denotémoslo como $p_0$. Ahora todos los demás puntos $p_1,\dots,p_n$ del polígono están ordenados por su ángulo polar desde el punto elegido (porque el polígono está ordenado en sentido antihorario).

Si el punto pertenece al polígono, pertenece a algún triángulo $p_0, p_i, p_{i + 1}$ (tal vez más de uno si yace en el borde de los triángulos). Consideremos el triángulo $p_0, p_i, p_{i + 1}$ tal que $p$ pertenece a este triángulo e $i$ es máximo entre todos esos triángulos.

Hay un caso especial. $p$ yace en el segmento $(p_0, p_n)$. Este caso lo comprobaremos por separado. En caso contrario todos los puntos $p_j$ con $j \le i$ están en sentido antihorario desde $p$ respecto de $p_0$, y todos los demás puntos no están en sentido antihorario desde $p$. Esto significa que podemos aplicar búsqueda binaria para el punto $p_i$, tal que $p_i$ no está en sentido antihorario desde $p$ respecto de $p_0$, e $i$ es máximo entre todos esos puntos. Y después comprobamos si el punto está realmente en el triángulo determinado.

El signo de $(a - c) \times (b - c)$ nos dirá si el punto $a$ está en sentido horario o antihorario desde el punto $b$ respecto del punto $c$. Si $(a - c) \times (b - c) > 0$, entonces el punto $a$ está a la derecha del vector que va de $c$ a $b$, lo que significa en sentido horario desde $b$ respecto de $c$. Y si $(a - c) \times (b - c) < 0$, entonces el punto está a la izquierda, o en sentido antihorario. Y está exactamente en la recta entre los puntos $b$ y $c$.

Volviendo al algoritmo: Consideremos un punto de consulta $p$. Primero, debemos comprobar si el punto yace entre $p_1$ y $p_n$. En caso contrario ya sabemos que no puede ser parte del polígono. Esto se puede hacer comprobando si el producto cruz $(p_1 - p_0)\times(p - p_0)$ es cero o tiene el mismo signo que $(p_1 - p_0)\times(p_n - p_0)$, y $(p_n - p_0)\times(p - p_0)$ es cero o tiene el mismo signo que $(p_n - p_0)\times(p_1 - p_0)$. Luego manejamos el caso especial en el que $p$ es parte de la recta $(p_0, p_1)$. Y luego podemos hacer búsqueda binaria del último punto de $p_1,\dots p_n$ que no está en sentido antihorario desde $p$ respecto de $p_0$. Para un solo punto $p_i$ esta condición se puede comprobar verificando que $(p_i - p_0)\times(p - p_0) \le 0$. Después de encontrar tal punto $p_i$, debemos probar si $p$ yace dentro del triángulo $p_0, p_i, p_{i + 1}$. Para probar si pertenece al triángulo, podemos simplemente comprobar que $|(p_i - p_0)\times(p_{i + 1} - p_0)| = |(p_0 - p)\times(p_i - p)| + |(p_i - p)\times(p_{i + 1} - p)| + |(p_{i + 1} - p)\times(p_0 - p)|$. Esto comprueba si el área del triángulo $p_0, p_i, p_{i+1}$ tiene exactamente el mismo tamaño que la suma de los tamaños del triángulo $p_0, p_i, p$, el triángulo $p_0, p, p_{i+1}$ y el triángulo $p_i, p_{i+1}, p$. Si $p$ está fuera, entonces la suma de esos tres triángulos será mayor que el tamaño del triángulo. Si está dentro, entonces será igual.

Implementación

La función prepare se asegurará de que el punto lexicográficamente más pequeño (menor valor x, y en empates menor valor y) sea $p_0$, y computa los vectores $p_i - p_0$. Después la función pointInConvexPolygon computa el resultado de una consulta. Además recordamos el punto $p_0$ y trasladamos todos los puntos consultados con él para computar la distancia correcta, ya que los vectores no tienen un punto inicial. Al trasladar los puntos de consulta podemos asumir que todos los vectores empiezan en el origen $(0, 0)$, y simplificar los cómputos de distancias y longitudes.

struct pt {
    long long x, y;
    pt() {}
    pt(long long _x, long long _y) : x(_x), y(_y) {}
    pt operator+(const pt &p) const { return pt(x + p.x, y + p.y); }
    pt operator-(const pt &p) const { return pt(x - p.x, y - p.y); }
    long long cross(const pt &p) const { return x * p.y - y * p.x; }
    long long dot(const pt &p) const { return x * p.x + y * p.y; }
    long long cross(const pt &a, const pt &b) const { return (a - *this).cross(b - *this); }
    long long dot(const pt &a, const pt &b) const { return (a - *this).dot(b - *this); }
    long long sqrLen() const { return this->dot(*this); }
};

bool lexComp(const pt &l, const pt &r) {
    return l.x < r.x || (l.x == r.x && l.y < r.y);
}

int sgn(long long val) { return val > 0 ? 1 : (val == 0 ? 0 : -1); }

vector<pt> seq;
pt translation;
int n;

bool pointInTriangle(pt a, pt b, pt c, pt point) {
    long long s1 = abs(a.cross(b, c));
    long long s2 = abs(point.cross(a, b)) + abs(point.cross(b, c)) + abs(point.cross(c, a));
    return s1 == s2;
}

void prepare(vector<pt> &points) {
    n = points.size();
    int pos = 0;
    for (int i = 1; i < n; i++) {
        if (lexComp(points[i], points[pos]))
            pos = i;
    }
    rotate(points.begin(), points.begin() + pos, points.end());

    n--;
    seq.resize(n);
    for (int i = 0; i < n; i++)
        seq[i] = points[i + 1] - points[0];
    translation = points[0];
}

bool pointInConvexPolygon(pt point) {
    point = point - translation;
    if (seq[0].cross(point) != 0 &&
            sgn(seq[0].cross(point)) != sgn(seq[0].cross(seq[n - 1])))
        return false;
    if (seq[n - 1].cross(point) != 0 &&
            sgn(seq[n - 1].cross(point)) != sgn(seq[n - 1].cross(seq[0])))
        return false;

    if (seq[0].cross(point) == 0)
        return seq[0].sqrLen() >= point.sqrLen();

    int l = 0, r = n - 1;
    while (r - l > 1) {
        int mid = (l + r) / 2;
        int pos = mid;
        if (seq[pos].cross(point) >= 0)
            l = mid;
        else
            r = mid;
    }
    int pos = l;
    return pointInTriangle(seq[pos], seq[pos + 1], pt(0, 0), point);
}

Problemas

• SGU253 Theodore Roosevelt 
• Codeforces 55E Very simple problem 
• Codeforces 166B Polygons