Change detection con M3C2 su nuvole di punti

L'algoritmo M3C2 viene utlizzato da Cloudcompare per effettuare il change detection ...volevo trovare una soluzione per rendere automatico il processing e la pubblicazione

Si parte da due LAS e tramite la libreria py4dgeo si ottiene un file las con i punti colorati secondo una scala colore basata sulla distanza tra i due LAS si origine

import numpy as np

import py4dgeo

import laspy

import sys

import math

py4dgeo.set_num_threads(4) # numero di threads attivi

epoch1, epoch2 = py4dgeo.read_from_las("2019.las", "2022.las")

corepoints = epoch1.cloud[::30]

m3c2 = py4dgeo.M3C2(epochs=(epoch1, epoch2),corepoints=corepoints,cyl_radii=(2.0,),normal_radii=(0.5, 1.0, 2.0),)

distances, uncertainties = m3c2.run()

dist_max = np.nanmax(distances)

dist_min = np.nanmin(distances)

reds = np.empty(corepoints.shape)

blues = np.empty(corepoints.shape)

greens = np.empty(corepoints.shape)

fattore =(dist_max-dist_min)/(dist_max+dist_min)

Red =0

Blue = 0

Green = 0

for i in range (distances.shape[0]):

if math.isnan(distances[i]):

t = 0

else:

#print(distances[i])

t = int(((distances[i]-dist_min)/(dist_max-dist_min)) * 16581375)

#print(t)

if math.isnan(t):

Blue = 0

Green = 0

Red = 0

else:

if (distances[i]< -1.5):

Red = 255

Green = 0

Blue = 0

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

Red = 255

Green = 255

Blue = 0

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

Red = 0

Green = 255

Blue = 0

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

Red = 0

Green = 255

Blue = 255

if (distances[i]> 1.5):

Red = 0

Green = 0

Blue = 255

#RGBint = t

#Red = RGBint & 255

#Blue = (RGBint >> 8) & 255

#Green = (RGBint >> 16) & 255'''

reds[i]= Red

greens[i]= Green

blues[i]=Blue

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

header.add_extra_dim(laspy.ExtraBytesParams(name="distance", type=np.float32))

xmin = np.floor(np.min(corepoints[:,0]))

ymin = np.floor(np.min(corepoints[:,1]))

zmin = np.floor(np.min(corepoints[:,2]))

header.offsets = [xmin,ymin,zmin]

#header.scales = np.array([0.1, 0.1, 0.1])

las = laspy.LasData(header)

las.x = corepoints[:,0]

las.y = corepoints[:,1]

las.z = corepoints[:,2]

las.red = reds[:,0]

las.blue = blues[:,0]

las.green = greens[:,0]

las.distance = np.copy(distances)

las.write("m3c2_result.las")

#./PotreeConverter m3c2_result.las -o /var/www/html/m3c2/ --generate-page index

Per la pubblicazione su web il file viene convertito e vestito tramite PoTree 

 

• Anders, K., Winiwarter, L., Lindenbergh, R., Williams, J. G., Vos, S. E., & Höfle, B. (2020). 4D objects-by-change: Spatiotemporal segmentation of geomorphic surface change from LiDAR time series. ISPRS Journal of Photogrammetry and Remote Sensing, 159, pp. 352-363. doi: 10.1016/j.isprsjprs.2019.11.025.

• Anders, K., Winiwarter, L., Mara, H., Lindenbergh, R., Vos, S. E., & Höfle, B. (2021). Fully automatic spatiotemporal segmentation of 3D LiDAR time series for the extraction of natural surface changes. ISPRS Journal of Photogrammetry and Remote Sensing, 173, pp. 297-308. doi: 10.1016/j.isprsjprs.2021.01.015.

• Truong, C., Oudre, L., Vayatis, N. (2018): ruptures: Change point detection in python. arXiv preprint: abs/1801.00826.