Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como encontrar o fecho convexo de um conjunto de pontos usando a estratégia de "divisão e conquista"?

+1 voto
105 visitas

1 Resposta

0 votos
respondida Jun 17 por Pedro Antero (26 pontos)  
editado Jul 6 por Pedro Antero

O fecho convexo, \(S^*\), de um conjunto de pontos, \(S\), é o menor polígono convexo que contém S. Nesta resposta, exploraremos duas estratégias para encontrarmos o fecho convexo de um grupo de pontos: (1) força bruta e (2) divisão e conquista.

O método “força bruta” consiste em se traçar retas \(r(i,j)\) conectando todos os pares de pontos \((i,j)\) pertencentes a \(S\). Neste caso, \(r(i,j)\) pertencerá a \(S^*\) se todos os demais pontos de \(S\backslash\{i,j\}\) se postarem apenas em um dos lados da reta. É fácil perceber que esta estratégia tem complexidade \(O(n^3)\): gasta-se \(O(n^2)\) para traçar as retas \(r(i,j)\) e \(O(n)\) para se encontrar as posições dos demais pontos relativas a \(r(i,j)\).

No código abaixo é apresentada uma implementação do método “força bruta”. Para facilitar a compreensão, os comentários específicos foram mantidos nas linhas do código.

# 1. Convex Hull Brute Force implementation

import numpy as np  # imports numpy library
import matplotlib.pyplot as plt  # imports pyplot, a MATLAB-like plotting framework

np.random.seed(150)  # sets the initial seed of the RNG

fig = plt.figure(num=None, figsize=(8, 8), dpi=80, facecolor='w', edgecolor='k')  # creates a new figure and sets its
# size (in inches), resolution and colors

epsilon = 0.1
plt.axis([0 - epsilon, 1 + epsilon, 0 - epsilon, 1 + epsilon])  # sets both Y and X axes to [-0.1, 0.1]


def convex_hull_bf(vector):  # defines the function that will generate the convex hull (brute force implementation)
    n = np.shape(vector)[0]  # number of points
    ch = []  # empty array to store the points that belong to the convex hull

    for i in range(n - 1):  # explores all the possible pairs (i, j) of the given points

        for j in range(i + 1, n):
            a = (vector[i, 1] - vector[j, 1])/(vector[i, 0] - vector[j, 0])  # calculates "a" and "b", the parameters
            b = vector[i, 1] - a * vector[i, 0]  # of the connecting line between "i" and "j" (as in a * x + b = 0)
            line_is_convex_hull = True  # dummy that states if the line belongs to the convex hull
            first = True
            k = 0

            while line_is_convex_hull and k < n:  # tests convexity for all the "n" points...

                if k != i and k != j:  # ...that are different from the points defining the current line (i, j)

                    if first:
                        sign_test = np.sign(vector[k, 1] - a * vector[k, 0] - b)
                        first = False  # determines the position of the 1st point (right/above or left/below the line)

                    else:
                        sign = np.sign(vector[k, 1] - a * vector[k, 0] - b)
                        if sign * sign_test < 0:  # checks if all the other points are in the same position as the 1st
                            line_is_convex_hull = False  # ...if NOT, then the line "i-j" does not belong to the ch
                k = k + 1

            if line_is_convex_hull:  # if the line belongs to the convex hull...
                ch.append([vector[i, :], vector[j, :]])  # ...append its points to the complex hull array
    return ch

number_of_points = 30  # sets the number of points
points = np.random.uniform(size=[number_of_points, 2])  # defines the (x,y) coordinates randomly
convex_hull = convex_hull_bf(points)  # calls the brute force convex hull function, defined above
plt.plot(points[:, 0], points[:, 1], 'b.', markersize=20)  # plots the points with blue markers

for l in range(len(convex_hull)):
    plt.plot([convex_hull[l][0][0], convex_hull[l][1][0]], [convex_hull[l][0][1], convex_hull[l][2][1]], 'r',
             markersize=20)  # plots the convex hull lines, colored in red
    plt.savefig('convexHull.eps')  # saves the plot in the .eps format

Na Figura 1, apresenta-se o “fecho convexo” obtido a partir da abordagem “força bruta”, acima, para uma lista \(S\) com 30 pontos.
A imagem será apresentada aqui.
Figura 1 - Abordagem "força bruta" para o fecho convexo: saída do algoritmo.

A abordagem de “divisão e conquista” baseia-se na divisão do conjunto \(S\) em subconjuntos menores de pontos. Assim, reduz-se o problema ao cálculo do fecho para cada um destes subconjuntos. O passo final (e menos trivial) é a combinação de cada um destes fechos para formação de \(S^*\). Abaixo, apresentaremos uma explicação resumida de cada um dos passos, baseada na excelente aula do Prof. Srinivas Devadas, do MIT (disponível em https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-046j-design-and-analysis-of-algorithms-spring-2015/lecture-videos/lecture-2-divide-conquer-convex-hull-median-finding/)

  • Passo 1: Divida \(S\) em dois subconjuntos \(A\) e \(B\), utilizando uma reta
    vertical (representada pela reta “\(L\)” na Figura 2);
  • Passo 2: Compute e combine os fechos convexos de \(A\) e \(B\);
  • Passo 2.a) Encontre os pontos (\(a_i\), \(b_j\)) que definem a tangente
    superior, com \(a_i\) \(\in\) \(A\) e \(b_j\)) \(\in\) \(B\). A tangente
    superior é aquela que maximiza a altura \(y\) de corte da reta vertical
    definida no Passo 1. Na Figura 2, os pontos \(a_i\) = \(a_4\) e \(b_j\) = \(b_2\) definem a
    tangente superior;
  • Passo 2.b) Encontre os pontos (\(a_k\), \(b_m\)) que definem a tangente
    inferior, com \(a_k\) \(\in\) \(A\) e \(b_m\) \(\in\) \(B\). A tangente
    inferior é aquela que minimiza a altura \(y\) de corte da reta vertical
    definida no Passo 1. Na Figura 2, os pontos \(a_k\) = \(a_3\) e \(b_m\) = \(b_3\) definem a
    tangente inferior;
  • Passo 2.c) Conecte \(a_i\) a \(b_j\);
  • Passo 2.d) A partir de \(b_j\), percorra \(B\) até encontrar \(b_m\);
  • Passo 2.e) Conecte \(b_m\) a \(a_k\) e percorra \(A\) até retornar a \(a_i\).

A imagem será apresentada aqui.
Figura 2 - Estratégia "divisão e conquista" (retirada de https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-046j-design-and-analysis-of-algorithms-spring-2015/lecture-notes/MIT6_046JS15_lec02.pdf)

O método mais direto para se encontrar as tangentes superiores e inferiores descritas nos passos 2.a) e 2.b) consiste da comparação das alturas de corte (\(y\)) para todas as combinações de pontos (\(a\), \(b\)), com \(a\) \(\in\) \(A\) e \(b\) \(\in\) \(B\). Essa implementação teria complexidade de \(O(n^2)\).

Alternativamente, usaremos uma implementação mais eficiente, proposta pelo Prof. Srinivas Devadas. Nesta abordagem, \(a_i\) é encontrado percorrendo-se \(A\) no sentido anti-horário, a partir do ponto de \(A\) que estiver mais próximo da reta vertical divisora de \(S\). \(b_i\), por sua vez, é encontrado percorrendo-se \(B\) no sentido horário, a partir do ponto de \(B\) mais próximo da reta vertical. Os pontos \(a_k\) e \(b_m\) são encontrados de modo análago: caminhado sobre \(A\) no sentido horário, para \(a_k\), e sobre \(B\) no sentido anti-horário, para \(b_m\). Utilizando-se esta estratégia, é possível reduzir a complexidade do algoritmo para \(T(n) = 2T(n/2) + \Theta(n) = \Theta(n*log (n))\).

A implementação é apresentada no código abaixo. Os comentários específicos foram mantidos nas linhas do código.

# 2. Convex Hull Divide to Conquer implementation
#
import numpy as np  # imports numpy library
import math # imports math library
import matplotlib.pyplot as plt  # imports pyplot, a MATLAB-like plotting framework

np.random.seed(150)  # sets the initial seed of the RNG

fig = plt.figure(num=None, figsize=(8, 8), dpi=80, facecolor='w', edgecolor='k')  # creates a new figure and sets its
# size (in inches), resolution and colors

epsilon = 0.1
plt.axis([0 - epsilon, 1 + epsilon, 0 - epsilon, 1 + epsilon])  # sets both Y and X axes to [-0.1, 0.1]


def clockwiseangle_and_distance(point, origin, refvec): # defines the auxiliary function that sorts the points clockwise
    angle_vec = [0] * len(point)
    angle = 0
    for i in range(len(point)):
        vector = [point[i][0]-origin[0], point[i][1]-origin[1]] # vector between point and the origin: v = p - o
        lenvector = math.hypot(vector[0], vector[1]) # length of vector: ||v||
        if lenvector == 0: # if length is zero there is no angle
            angle = 0
        else:
            normalized = [vector[0]/lenvector, vector[1]/lenvector] # normalize vector: v/||v||
            dotprod  = normalized[0]*refvec[0] + normalized[1]*refvec[1] # x1*x2 + y1*y2
            if refvec[1] >= 0:
                diffprod = refvec[1]*normalized[0] - refvec[0]*normalized[1] # x1*y2 - y1*x2
            else:
                diffprod = refvec[0] * normalized[1] - refvec[1] * normalized[0]
            angle = math.atan2(diffprod, dotprod)
            if angle < 0:
                angle =  2 * math.pi + angle # negative angles represent counter-clockwise angles so we need to
                # subtract them from 2*pi (360 degrees)
        angle_vec[i] = angle
    r_angle_vec = angle_vec
    return r_angle_vec


def y(x1, x2, x_r): # defines the auxilary function that calculates the height where the tangents cross the "x_r" line
    a = (x2[1] - x1[1]) / (x2[0] - x1[0]) # tangent line "a" and "b" coefficients
    b = x1[1] - (a * x1[0])
    return (a * x_r) + b


def merge_ch(ch_s1, ch_s2): # defines the function that will merge the points to the left and to the right of "x_r",
    # (ch_s1 and ch_s2)
    ch = [] # declares "ch", to store the points that belong to the convex hull

    ch_s1 = sorted(ch_s1, key=lambda l: l[0]) # sorts points according to their "X" coordinates
    ch_s2 = sorted(ch_s2, key=lambda l: l[0])

    if len(ch_s1) == 0:
        x_r = ch_s2[0][0] # defines the vertical splitting line (special case when there are only points to the left)

    elif len(ch_s2) == 0:
        x_r = ch_s1[len(ch_s2)-1][0] # defines the vertical splitting line (special case when there are only points to
        # the right)

    else:
        x_r = (ch_s1[len(ch_s1)-1][0] + ch_s2[0][0]) / 2  # X coordinate of the vertical line that splits "s" into "s1"
        # and "s2"

    # ---------------------------------------- Upper tangent ---------------------------------------------------------#
    an_ch_s1 = clockwiseangle_and_distance(ch_s1, ch_s1[len(ch_s1) - 1], [0, -1])
    an_ch_s2 = clockwiseangle_and_distance(ch_s2, ch_s2[0], [0, 1])

    ch_s1 = np.vstack([x for (y, x) in sorted(zip(an_ch_s1, ch_s1))])  # sorts points to the left of "x_r" (ch_s1)
    # counter-clockwise
    ch_s2 = np.vstack([x for (y, x) in sorted(zip(an_ch_s2, ch_s2))]) # sorts points to the right of "x_r" (ch_s2)
    # clockwise

    i = 0 # "i" is the left point of the upper tangent
    j = 0 # "j" is the right point of the upper tangent

    while (j < len(ch_s2) - 1 or i < len(ch_s1) - 1) and \
            (y(ch_s1[min(i, len(ch_s1) - 1)], ch_s2[min(j + 1, len(ch_s2)-1)],
               x_r) > y(ch_s1[min(i, len(ch_s1) - 1)], ch_s2[min(j, len(ch_s2) - 1)], x_r)
             or y(ch_s1[min(i + 1, len(ch_s1) - 1)], ch_s2[min(j, len(ch_s2) - 1)],
                  x_r) > y(ch_s1[min(i, len(ch_s1) - 1)], ch_s2[min(j, len(ch_s2) - 1)], x_r)):  # finds upper tangent

        if y(ch_s1[min(i, len(ch_s1) - 1)], ch_s2[min(j + 1, len(ch_s2)-1)],
             x_r) > y(ch_s1[min(i, len(ch_s1) - 1)], ch_s2[min(j, len(ch_s2) - 1)], x_r):
            j = j + 1

        else:
            j = 0
            i = i + 1

    up_tg = [ch_s1[i], ch_s2[j]]

    # ---------------------------------------- Lower tangent ---------------------------------------------------------#
    an_ch_s1 = clockwiseangle_and_distance(ch_s1, ch_s1[0], [0, 1])
    an_ch_s2 = clockwiseangle_and_distance(ch_s2, ch_s2[0], [0, -1])

    ch_s1 = np.vstack([x for (y, x) in sorted(zip(an_ch_s1, ch_s1))]) # sorts points to the left of "x_r" (ch_s1)
    # clockwise
    ch_s2 = np.vstack([x for (y, x) in sorted(zip(an_ch_s2, ch_s2))]) # sorts points to the right of "x_r" (ch_s2)
    # counter-clockwise

    m = 0 # "m" is the left point of the lower tangent
    k = 0 # "k" is the right point of the upper tangent

    while (m < len(ch_s2) - 1 or k < len(ch_s1) - 1) and \
            (y(ch_s1[min(k, len(ch_s1) - 1)], ch_s2[min(m + 1, len(ch_s2) - 1)],
               x_r) < y(ch_s1[min(k, len(ch_s1) - 1)], ch_s2[min(m, len(ch_s2) - 1)], x_r) or
                     y(ch_s1[min(k + 1, len(ch_s1) - 1)], ch_s2[min(m, len(ch_s2) - 1)],
                       x_r) < y(ch_s1[min(k, len(ch_s1) - 1)], ch_s2[min(m, len(ch_s2) - 1)], x_r)):  # finds lower
        # tangent

        if y(ch_s1[min(k, len(ch_s1) - 1)], ch_s2[min(m + 1, len(ch_s2)-1)],
             x_r) < y(ch_s1[min(k, len(ch_s1) - 1)], ch_s2[min(m, len(ch_s2) - 1)], x_r):
            m = m + 1

        else:
            m = 0
            k = k + 1

    lw_tg = [ch_s1[k], ch_s2[m]]

    # ---------------------------------------- Merging process ------------------------------------------------------ #
    an_ch_s1 = clockwiseangle_and_distance(ch_s1, ch_s1[0], [0, -1])
    an_ch_s2 = clockwiseangle_and_distance(ch_s2, ch_s2[0], [0, 1])

    ch_s1 = np.vstack([x for (y, x) in sorted(zip(an_ch_s1, ch_s1))]) # sorts points to the left of "x_r" (ch_s1)
    # counter-clockwise
    ch_s2 = np.vstack([x for (y, x) in sorted(zip(an_ch_s2, ch_s2))]) # sorts points to the right of "x_r" (ch_s2)
    # clockwise

    i = np.where(ch_s1 == up_tg[0]) # retrieves "i", "j", "k", "m", after the above sorting process
    j = np.where(ch_s2 == up_tg[1])
    k = np.where(ch_s1 == lw_tg[0])
    m = np.where(ch_s2 == lw_tg[1])

    i = int(i[0][0])
    j = int(j[0][0])
    k = int(k[0][0])
    m = int(m[0][0])

    ch.append(ch_s1[i]) # appends the upper tangent points to the convex hull
    ch.append(ch_s2[j])
    t = 1
    u = 1

    if ch_s1[k][0] != ch_s1[i][0] or ch_s1[k][6] != ch_s1[i][7] or \
                    ch_s2[m][0] != ch_s2[j][0] or ch_s2[m][8] != ch_s2[j][9]:
            while j + t <= len(ch_s2) - 1:
                if ch_s2[j + t][0] != ch_s2[m][0] or ch_s2[j + t][10] != ch_s2[m][11]:
                    ch.append(ch_s2[j + t]) # appends the right side points that belong to the convex hull to "ch"
                    t = t + 1
                else:
                    t = len(ch_s2) - j


            if ch_s2[m][0] != ch_s2[j][0] or ch_s2[m][12] != ch_s2[j][13]:
                ch.append(ch_s2[m]) # appends the lower tangent right point to the convex hull

            if ch_s2[k][0] != ch_s2[i][0] or ch_s2[k][14] != ch_s2[i][15]:
                ch.append(ch_s1[k]) # appends the lower tangent left point to the convex hull

            an_ch_s1 = clockwiseangle_and_distance(ch_s1, ch_s1[0], [0, 1]) # sorts points to the left of "x_r" (ch_s1)
            #  clockwise, according to the left point of the lower tanget
            ch_s1_k = np.vstack([x for (y, x) in sorted(zip(an_ch_s1, ch_s1))])

            k = np.where(ch_s1_k == lw_tg[0])
            k = int(k[0][0])
            i = np.where(ch_s1 == up_tg[0])
            i = int(i[0][0])

            if k + 1 <= len(ch_s1) - 1:

                while u + k <= len(ch_s1) - 1:
                    if ch_s1_k[u + k][0] != ch_s1[i][0] or ch_s1_k[u + k][16] != ch_s1[i][17]:
                        ch.append(ch_s1_k[u + k]) # appends the left side points that belong to the convex hull to "ch"
                        u = u + 1
                    else:
                        u = len(ch_s1) - k


    s = np.vstack(ch)
    return s


def convex_hull_dq(s):  # defines the function that will generate the convex hull (divide to conquer implementation)
    n = np.size(s,0)
    if n > 1:
        s = np.vstack(sorted(s, key=lambda l: l[0])) # sorts the vector of points based on the X coordinate
        s1 = s[:len(s) // 2]  # divides the set of points "s" into two subsets "s1" and "s2": "s1" ("s2") is on the left
        s2 = s[len(s) // 2:]  # (right) side of the vertical line that divides the set "s"
        ch_s1 = convex_hull_dq(s1)  # calculates the convex hull of the subset "s1"
        ch_s2 = convex_hull_dq(s2)  # calculates the convex hull of the subset "s2"
        s = merge_ch(ch_s1, ch_s2)

    return s

number_of_points = 30 # sets the number of points
points = np.random.uniform(size=[number_of_points, 2])  # defines the (x,y) coordinates randomly

convex_hull = convex_hull_dq(points)  # calls the divide to conquer convex hull function, defined above
convex_hull = np.vstack((convex_hull, convex_hull[0]))
plt.plot(points[:, 0], points[:, 1], 'b.', markersize=20)  # plots the points with blue markers

for l in range(len(convex_hull)-1):
    plt.plot([convex_hull[l][0], convex_hull[l + 1][0]], [convex_hull[l][18], convex_hull[l + 1][19]], 'r',
             markersize=20)  # plots the convex hull lines, colored in red
    plt.savefig('convexHull.eps')  # saves the plot in the .eps format

A saída da abordagem de “divisão e conquista” é apresentada na Figura 3 para um conjunto \(S\) de 30 pontos. O resultado é idêntico àquele da Figura 1.

A imagem será apresentada aqui.
Figura 3 - Abordagem "divisão e conquista" para o fecho convexo: saída do algoritmo.

comentou Jul 5 por Caue (226 pontos)  
A solução está ótima. Gostaria apenas de fazer um comentário:
Na chamada vector = [point[i][0]-origin[0], point[i][5]-origin[1]], não seria vector = [point[i][0]-origin[0], point[i][1]-origin[1]], por se tratar de um vetor (x,y)?
Isso eu já identifiquei em mais duas respostas, então acho que o prorum em algumas situações está incrementando os números entre [].

Aproveito e deixo uma dica para quando precisar usar coordenadas (x,y): namedtuple
https://docs.python.org/3.6/library/collections.html#collections.namedtuple
comentou Jul 6 por Pedro Antero (26 pontos)  
Olá, Cauê!
De fato, por algum motivo o Prorum incrementou o número entre []. Já corrigi o código.
Obrigado!
...