changing the timestamps that are written to fullwave files. (#71)

* changing the timestamps that are written to fullwave files.

* adding fullwave conversion script in first version

* added explanations

* version bump

* fixing shifted position due to waveform time shift,
added one bin safety margin

pyhelios_demo/txt2las_wdp.py

@@ -0,0 +1,78 @@
+#!/usr/bin/python
+# txt2las_wdp
+# copyright 2021 3DGeo Heidelberg
+# Script to transform a HELIOS++ xyz point cloud and fullwave file into las1.4 and external WDP file
+# Result can be viewed e.g. with CloudCompare (open file with "open"-> select "las 1.3 and las 1.4" as file type)
+
+import struct
+import laspy
+import numpy as np
+from collections import OrderedDict
+
+wfs = OrderedDict()
+with open(r"..\output\Survey Playback\als_hd_demo\points\leg000_fullwave.txt", 'r') as inf:
+    for line in inf.readlines():
+        data = [float(x) for x in line.split(" ")]
+        wf_idx = int(data[0])
+        tmin = data[7]
+        tmax = data[8]
+        fwf = np.abs(np.array(data[10:]))
+        fwf = fwf/np.max(fwf) * 255.
+        fwf = fwf.astype(int)
+        dx, dy, dz = data[4], data[5], data[6]
+        wfs[wf_idx] = [tmin, tmax, dx, dy, dz, fwf]
+
+max_len = max([len(x[-1][-1]) for x in wfs.items()])
+wfs_lut = {idx: val for idx, val in enumerate(np.unique(wfs.keys())[0])}
+wfs_lut_rev = {val: idx for idx, val in wfs_lut.items()}
+pts = np.loadtxt(r"..\output\Survey Playback\als_hd_demo\points\leg000_points.xyz", delimiter=" ")
+
+las = laspy.create(file_version="1.4", point_format=4)
+xyz = pts[:, :3]
+las.header.offsets = np.min(xyz, axis=0)
+las.header.scales = [0.001, 0.001, 0.001]
+las.x = xyz[:, 0]
+las.y = xyz[:, 1]
+las.z = xyz[:, 2]
+las.intensity = pts[:, 3]
+las.wavepacket_index=np.ones((xyz.shape[0]))  # 1 refers to the value 100 (99+1) in the VLR below
+las.wavepacket_offset = 60 + np.array([wfs_lut_rev[idx] for idx in pts[:, 7]]) * (max_len)
+las.wavepacket_size = [max_len] * xyz.shape[0]  # always write a wavepacket of the longest size
+las.return_point_wave_location = [0] * xyz.shape[0]  # adapt this if you need the information on which position this point
+                                                     # represents in the waveform
+las.x_t = [wfs[idx][2] for idx in pts[:, 7]]
+las.y_t = [wfs[idx][3] for idx in pts[:, 7]]
+las.z_t = [wfs[idx][4] for idx in pts[:, 7]]
+
+VLR_data = struct.pack("<BBLLdd",
+                        8,  # bits per sample
+                        0,  # compression type
+                        max_len,  # number of samples
+                        10,  # temporal sample spacing (ps)
+                        1,  # digitizer gain
+                        0,  # digitizer offset
+                       )
+new_vlr = laspy.VLR(user_id="LASF_Spec", record_id=100, record_data=VLR_data)
+las.vlrs.append(new_vlr)
+
+las.write("test.las")
+
+wf_data = bytes()  # WDP header (same as EVLR header)
+wf_data = b''.join([wf_data, struct.pack("<H", 0)])  #reserved
+wf_data = b''.join([wf_data, struct.pack("<s", b"LASF_Spec       ")])  #user id
+wf_data = b''.join([wf_data, struct.pack("<H", 65535)])  #record id
+wf_data = b''.join([wf_data, struct.pack("<Q", max_len*len(wfs))])  #record length after header
+wf_data = b''.join([wf_data, struct.pack("<s", b" "*32)])  #descripton
+
+with open("test.wdp", 'wb') as f:
+    f.write(wf_data)
+
+    for wf in wfs.items():
+        byt = np.zeros((max_len,), dtype=int)
+        byt[:len(wf[1][-1])] = [int(d) for d in wf[1][-1]]
+        f.write(struct.pack("<%dB" % max_len, *byt))
+
+
+with open("test.las", 'r+b') as f:
+    f.seek(6)
+    f.write(struct.pack("<H", int('0000000000000100', 2))) # set "external waveform" bit to 1 in las file (not possible with laspy)

src/helios_version.cpp

@@ -1,6 +1,6 @@
 #include <helios_version.h>
 
-const char * HELIOS_VERSION = "1.0.8";
+const char * HELIOS_VERSION = "1.0.9";
 
 const char * getHeliosVersion(){
     return HELIOS_VERSION;

src/scanner/detector/FullWaveformPulseRunnable.cpp

@@ -269,7 +269,8 @@ void FullWaveformPulseRunnable::digestIntersections(
         fullwave,
         distanceThreshold,
         minHitTime_ns,
-        nsPerBin
+        nsPerBin,
+        peakIntensityIndex
     );
 
     // Digest full waveform, generating measurements
@@ -337,12 +338,17 @@ bool FullWaveformPulseRunnable::initializeFullWaveform(
 ){
     // Calc time at minimum and maximum distance
     // (i.e. total beam time in fwf signal)
-    minHitTime_ns = minHitDist_m / cfg_speedOfLight_mPerNanosec;
-    maxHitTime_ns = maxHitDist_m / cfg_speedOfLight_mPerNanosec;
+
+    peakIntensityIndex = detector->scanner->peakIntensityIndex;
+    nsPerBin = detector->scanner->FWF_settings.binSize_ns;
+
+    minHitTime_ns = minHitDist_m / cfg_speedOfLight_mPerNanosec - peakIntensityIndex * nsPerBin; // time until first maximum minus rising flank
+    maxHitTime_ns = maxHitDist_m / cfg_speedOfLight_mPerNanosec + // time until last maximum
+            detector->scanner->getPulseLength_ns() - peakIntensityIndex * nsPerBin + // time for signal decay
+            nsPerBin; // 1 bin for buffer
 
     // Calc ranges and threshold
-    double hitTimeDelta_ns = maxHitTime_ns - minHitTime_ns +
-                             detector->scanner->getPulseLength_ns();
+    double hitTimeDelta_ns = maxHitTime_ns - minHitTime_ns;
     double maxFullwaveRange_ns =
         detector->scanner->FWF_settings.maxFullwaveRange_ns;
     distanceThreshold = maxHitDist_m;
@@ -359,9 +365,7 @@ bool FullWaveformPulseRunnable::initializeFullWaveform(
     }
 
     // Compute fullwave variables
-    nsPerBin = detector->scanner->FWF_settings.binSize_ns;
     numFullwaveBins = (int)(hitTimeDelta_ns / nsPerBin);
-    peakIntensityIndex = detector->scanner->peakIntensityIndex;
 
     return true;
 }
@@ -371,7 +375,8 @@ void FullWaveformPulseRunnable::populateFullWaveform(
     std::vector<double> &fullwave,
     double distanceThreshold,
     double minHitTime_ns,
-    double nsPerBin
+    double nsPerBin,
+    int peakIntensityIndex
 ){
     // Multiply each sub-beam intensity with time_wave and
     // add to the full waveform
@@ -382,7 +387,7 @@ void FullWaveformPulseRunnable::populateFullWaveform(
         if(entryDistance_m > distanceThreshold) continue;
         double entryIntensity = it->second;
         double wavePeakTime_ns = (entryDistance_m / cfg_speedOfLight_mPerNanosec); // [ns]
-        int binStart = (int)((wavePeakTime_ns-minHitTime_ns) / nsPerBin);
+        int binStart = (int)((wavePeakTime_ns-minHitTime_ns) / nsPerBin) - peakIntensityIndex;
         for (size_t i = 0; i < time_wave.size(); i++) {
             fullwave[binStart + i] += time_wave[i] * entryIntensity;
         }
@@ -440,7 +445,7 @@ void FullWaveformPulseRunnable::digestFullWaveform(
 
         // Compute distance
         double distance = cfg_speedOfLight_mPerNanosec *
-            ((i-peakIntensityIndex) * nsPerBin + minHitTime_ns);
+            (i * nsPerBin + minHitTime_ns);
 
         // Build list of objects that produced this return
         double minDifference = numeric_limits<double>::max();

src/scanner/detector/FullWaveformPulseRunnable.h

@@ -229,7 +229,8 @@ class FullWaveformPulseRunnable : public AbstractPulseRunnable {
         std::vector<double> &fullwave,
         double distanceThreshold,
         double minHitTime_ns,
-        double nsPerBin
+        double nsPerBin,
+        int peakIntensityIndex
     );
     /**
      * @brief Digest a previously populated full waveform vector,