Posizione date le distanze da due punti conosciuti

 Un problema che si ha con il posizionamento decawave e' quello di ottenere la posizione (coordinate) data la distanza da due punti conosciuti

Alla fine il problema si riduce ad ottenere gli angoli di un triangolo di cui sono conosciute le lunghezze di tutti i lati

diciamo che il lato A sia la distanza tra i due anchor mentre C e B sono le distanze del tag rispetto ai due anchor

Tutti i lati possono essere riscritti come

A = Ccos(alfa)+Bcos(beta)

B = Acos(beta)+Ccos(gamma)

C=Acos(alfa)+Bcos(gamma)

si tratta si un sistema di tre equazioni con 3 incognite e quindi risolvibile. Per semplice sostituzione si arriva alla fine che

cos(beta) = (A*A)+(B*B)-(C*C)/2AB

cos(alfa) = (A*A)+(C*C)-(B*B)/2AC

cos(gamma) = (B*B)+(C*C)-(A*A)/2BC

a questo punto abbiamo tutti gli angoli. Si puo' spostare un anchor sulle origine degli assi e ruotare gli assi in modo che il secondo anchor sia in coordinata (x2,0)

A questo punto le coordinate incognite del punto tra B e C sono

X = Bcos(alfa)

Y = SQRT(1 - cos(alfa)*cos(alfa))

Makerfabs DW3000

Aggiornamento

Frugando meglio nei sorgenti della libreria ho visto che sono previste 33 configurazioni dell'antenna in base al tranfer rate, alla lunghezza dei pacchetti ed altro. Di default e' attivato il profilo 33 che non risulta essere adatto per le lunghe distanze

Ho quindi modificato il file dw3000_config_options.h della libreria alla riga 240 con 

#define CONFIG_OPTION_16

che corrisponde alla impostazione

/* Configuration option 16.
 * Channel 9, PRF 64M, Preamble Length 1024, PAC 8, Preamble code 10, Data Rate 850k, STS Length 64
 */
dwt_config_t config_options = {
    9,                  /* Channel number. */
    DWT_PLEN_1024,      /* Preamble length. Used in TX only. */
    DWT_PAC8,           /* Preamble acquisition chunk size. Used in RX only. */
    10,                  /* TX preamble code. Used in TX only. */
    10,                  /* RX preamble code. Used in RX only. */
    3,                  /* 0 to use standard 8 symbol SFD, 1 to use non-standard 8 symbol, 2 for non-standard 16 symbol SFD and 3 for 4z 8 symbol SDF type */
    DWT_BR_850K,        /* Data rate. */
    DWT_PHRMODE_STD,    /* PHY header mode. */
    DWT_PHRRATE_STD,    /* PHY header rate. */
    (1024 + 1 + 8 - 8),  /* SFD timeout (preamble length + 1 + SFD length - PAC size). Used in RX only. */
    DWT_STS_MODE_1,      /* Mode 1 STS enabled */
    DWT_STS_LEN_64,      /* STS length*/
    DWT_PDOA_M0         /* PDOA mode off */
};
#endif

in pratica ho impostato il valore piu' basso di transfer rate (850Kbps con il pacchetto di preambolo della massima lunghezza ed abilitando STS. La modalita' STS (scrambled timestamp sequence) migliora la qualita' della stima della distanza impedendo interferenze

Ho utilizzato gli sketch di initiator e responder negli esempi della libreria ex_06a_ss_twr_initiator_sts e ex_06a_ss_twr_responder_sts senza modifiche

Da una serie di 1850 misure statiche la standard deviation e' risultata intorno ai 4 cm 

Ho provato a vedere spostando di 3 cm per volta se riuscivo ad osservare lo spostamento (teoricamente non doveva essere possibile perche' la sensibilita' e' data genericamente intorno ai 10 cm dalla ditta)

Pos 1 : 930 misure Media 2.66 m StDev 0.062 m

Pos 2  (Pos1 + 3 cm) : 670 misure Media 2.71 m StDev 0.025 m

Pos 2  (Pos1 + 7 cm) : 1480 misure Media 2.72 m StDev 0.027 m

Un primo aspetto da indicare e' che la distribuzione dei dati non risulta particolarmente simmetrica rispetto al valore medio 

Con questa configurazione sono arrivato a ricevere fino a circa 12 m di distanza  

==================================================================

Ho avuto modo di provare la scheda Makersfab Dw3000 basata su un Esp32 a cui e' collegata una antenna Decawave UWB Dw3000 (3110 per l'esattezza) che permette il calcolo della distanza tramite il calcolo del tempo di doppia percorrenza di un segnale tra due antenna (una chiamata initiator responsabile dell'invio dello scambio e ed un responser)...esiste una versione che usa la differenza di fase 

Questo metodo di calcolo del tempo di volo permette di non necessitare la sincronizzazione degli orologi delle due schede

Note estremamente dolenti...rispetto alla versione Dw1000 sempre di Makerfabs al momento il supporto software e di documentazione e' al limite dell'inesistente. Nella libreria Arduino ci sono alcuni esempi poco commentati (per avere informazioni meglio riferirsi all'SDK ufficiale )

Ho fatto qualche prova con lo sketch ex_06a_ss_twr_initiator_sts_no_data ma la distanza stimata e' molto instabile con variazione di circa 10 cm (e fino a qui rientra nell'errore dichiarato del sensore) ma molto sensibile anche dalla presenza di oggetti nell'intorno. Non sono stato in grado di avere una retta di calibrazione

Un altro aspetto che non riesco a gestire e' la potenza del segnale. A questo link viene indicato che si possano registrare pacchetti fino a circa 20 m di distanza mentre con i miei moduli non riesco ad allontanarmi oltre i 5 m

Alla fine non sono riuscito a capire se l'insuccesso e' relativo ad un difetto delle mie schede o qualche programmazione non corretta (ma di fatto ho usato gli sketch di esempio)

Telit GE864-PY

Mi e' stato regalato una volta dismesso questo dispositivo nato per leggere una seriale ed inviare i dati via rete cellulare 

Pensavo di poter recuperare qualche componente ma tutte le funzioni sono integrate nell'unico chip

Seagate ST-251

 Il mio piccolo museo dell'informatica ha avuto una nuova acquisizione, un hard disk Seagate ST-251 del 1990

Si tratta di un HD da 42 Mb formattato che ha connettori MTM e non IDE. Non ho il controller adatto per cui non ho la piu' pallida idea se sia funzionante (ma guardando su Internet sembra che questi dispositivi fossero particolarmente delicati e pochi sono sopravvissuti funzionanti)

PIP Python error: externally-managed-environment

 Da qualche tempo sulla mia Debian testing non riesco ad installare le librerie Python tramite PIP ma sono tramite APT a causa dell'errore

error: externally-managed-environment

× This environment is externally managed
╰─> To install Python packages system-wide, try apt install python3-xyz, where xyz is the package you are trying to install.

l'unico modo per risolvere e' quello di usare un venv od usare PyCharm che in automatico crea un venv per ogni progetto

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")

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.