Semana 14

Como mencionado na semana anterior, será usada a biblioteca Open3D para criar a aplicação de cálculo de volumes. Assim, esta semana, foi focada na construção do algoritmo para este propósito.

Para efetuarmos o cálculo de volumes em Open3D, devemos seguir os seguintes passos:

​

• Visualizar point cloud

 

Ler a point cloud

pc=o3d.io.read_point_cloud("data/… ")

Visualizar eixos

axes=o3d.geometry.TriangleMesh.create_coordinate_frame()

o3d.visualization.draw_geometries([pc, axes])

​

• Filtrar point cloud

 

Para filtrar a point cloud, ou seja, remover os pontos desnecessários, como os pontos representados no chão e outras interferências, podemos criar plano e eliminar os pontos que estejam a uma certa distância deste.

​

plane_model, inliers = pc.segment_plane(distance_threshold=…,

                                           ransac_n=…,

                                           num_iterations=…)

[a, b, c, d] = plane_model

 

plane_pc = pc.select_by_index(inliers)

 

stockpile_pc = pc.select_by_index(inliers, invert=True)

 

o3d.visualization.draw_geometries([stockpile_pc])

 

 

Para remover pontos próximos da point cloud:

 

cl, ind = stockpile_pc.remove_statistical_outlier(nb_neighbors=…,

                                                   std_ratio=…)

stockpile_pc = stockpile_pc.select_by_index(ind)

o3d.visualization.draw_geometries([stockpile_pc])

​

• ·       Criar malha triangular

 

Será necessário reconstruir a superfície da point cloud por triangulação. Esta triangulação será projetada no plano XY por remoção dos valores das coordenadas z. A triangulação será baseada no algoritmo de delaunay.

 

downpc = stockpile_pc.voxel_down_sample(voxel_size=…)

xyz = np.asarray(downpc.points)

xy = []

for point in xyz:

   xy.append([point[0], point[1]])

tri = Delaunay(np.array(xy))

 

surface = o3d.geometry.TriangleMesh()

surface.vertices = o3d.utility.Vector3dVector(xyz)

surface.triangles = o3d.utility.Vector3iVector(tri.simplices)

o3d.visualization.draw_geometries([surface], mesh_show_wireframe=True)

 

Através deste código obtemos os índices de cada triângulo e visualização da malha projetada.

 

• ·       Cálculo do Volume

 

Depois da reconstrução da malha, podemos começar com o cálculo de volumes. Diferentes abordagens podem ser seguidas. Vou calcular o volume de cada triângulo da superfície para o plano XY (altura 0). Depois é só somar os volumes de todos os triângulos.

Triângulos em malhas Open3D definem os índices dos vértices, não os vértices, então precisamos converter os triângulos para que os seus valores sejam, na verdade, pontos 3D e não índices para a lista de vértices:

 

def get_triangles_vertices(triangles, vertices):

   triangles_vertices = []

   for triangle in triangles:

       new_triangles_vertices = [vertices[triangle[0]], vertices[triangle[1]], vertices[triangle[2]]]

       triangles_vertices.append(new_triangles_vertices)

   return np.array(triangles_vertices)

​

​

Precisamos de uma função que faça computação do volume para cada triângulo.

​

def volume_under_triangle(triangle):

   p1, p2, p3 = triangle

   x1, y1, z1 = p1

   x2, y2, z2 = p2

   x3, y3, z3 = p3

   return abs((z1+z2+z3)*(x1*y2-x2*y1+x2*y3-x3*y2+x3*y1-x1*y3)/6)

​

volume = reduce(lambda a, b:  a + volume_under_triangle(b),

get_triangles_vertices(surface.triangles, surface.vertices), 0)

​

​

Até ao próximo relatório semanal! :)