Change detection con Open3D su nuvole di punti

Anche la libreria Open3D ha una funzione che permette di calcolare la distanza tra due nuvole di punti anche se non e' esplicitato l'algoritmo (inoltre per i punti in cui non c'e' sovrapposizione immediata il codice del precedente post riporta NaN mentre qui un valore di distanza viene sempre proposto. Inoltre la matrice delle distanze contiene solo valori positivi

I risultati mi sembrano migliori quelli dell'algoritmo del precedente post

import numpy as np

import laspy as lp

import open3d as o3d

p2019= lp.read("2019.las")

points_2019 = np.vstack((p2019.x, p2019.y, p2019.z)).transpose()

#colors_2019 = np.vstack((p2019.red, p2019.green,p2019.blue)).transpose()

pcd2019 = o3d.geometry.PointCloud()

pcd2019.points = o3d.utility.Vector3dVector(points_2019)

#pcd2019.colors = o3d.utility.Vector3dVector(colors_2019 / 65535)

p2022 = lp.read("2022.las")

points_2022 = np.vstack((p2022.x, p2022.y, p2022.z)).transpose()

#colors_2022 = np.vstack((p2022.red, p2022.green,p2022.blue)).transpose()

pcd2022 = o3d.geometry.PointCloud()

pcd2022.points = o3d.utility.Vector3dVector(points_2022)

distance=pcd2019.compute_point_cloud_distance(pcd2022)

rows, columns = points_2019.shape

reds = np.empty(rows)

blues = np.empty(rows)

greens = np.empty(rows)

for i in range(0,rows):

#print(t)

if (distance[i]< -1.5):

Red = 255

Green = 0

Blue = 0

if ((distance[i]>=-1.5) and (distance[i]<=-0.75)):

Red = 255

Green = 255

Blue = 0

if ((distance[i]>=-0.75) and (distance[i]<=0.75)):

Red = 0

Green = 255

Blue = 0

if ((distance[i]>=0.75) and (distance[i]<=1.5)):

Red = 0

Green = 255

Blue = 255

if (distance[i]> 1.5):

Red = 0

Green = 0

Blue = 255

if (distance[i]> 10):

Red = 0

Green = 0

Blue = 0

if (distance[i]< 10):

Red = 0

Green = 0

Blue = 0

reds[i]= Red

greens[i]= Green

blues[i]=Blue

header = lp.LasHeader(point_format=3, version="1.2")

las = lp.LasData(header)

las.x = points_2019[:,0]

las.y = points_2019[:,1]

las.z = points_2019[:,2]

las.red = reds[:]

las.blue = blues[:]

las.green = greens[:]

las.write("open3d_result.las")

Estrazione di piani da pointcloud

 Avevo gia' provato a fare una estrazione dei piani su un affioramento roccioso con PCL qualche tempo fa

Questa volta ho provato ad utilizzare la libreria Open3D con l'algoritmo Ransac per l'estrazione di piani multipli da una scansione fatta con il tablet Tango 

Per l'algoritmo e' sufficiente definire il numero massimo di piani attesi, la distanza di soglia ed il numero di iterazioni

I vari piani oltre ad essere visualizzati vengono anche estratti come ply

import sys
import numpy as np
import open3d as o3d
import matplotlib.pyplot as plt
input_path="./"
dataname="maiano.ply"
pcd = o3d.io.read_point_cloud(input_path+dataname)
o3d.visualization.draw_geometries([pcd])

#ransac piani multipli
segment_models={}
segments={}
max_plane_idx=3
rest=pcd
for i in range(max_plane_idx):
    colors = plt.get_cmap("tab20")(i)
    segment_models[i], inliers = rest.segment_plane(distance_threshold=0.1,ransac_n=3,num_iterations=5000)
    segments[i]=rest.select_by_index(inliers)
    segments[i].paint_uniform_color(list(colors[:3]))
    rest = rest.select_by_index(inliers, invert=True)
    print("pass",i,"/",max_plane_idx,"done.")
o3d.io.write_point_cloud("0.ply", segments[0])
o3d.visualization.draw_geometries([segments[0]] + [rest])
o3d.io.write_point_cloud("1.ply", segments[1])
o3d.visualization.draw_geometries([segments[1]] + [rest])
o3d.io.write_point_cloud("2.ply", segments[2])
o3d.visualization.draw_geometries([segments[2]] + [rest])