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.
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\).
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.
Figura 3 - Abordagem "divisão e conquista" para o fecho convexo: saída do algoritmo.