full upload so not to lose anything important

tools/animate_score.py

@@ -0,0 +1,151 @@
+from functools import partial
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib import animation
+
+# from util import calculate_average_frame_rate, load_dataset
+from rich.progress import Progress, track
+
+scores_path = Path(
+    "/home/fedex/mt/projects/thesis-kowalczyk-jan/Deep-SAD-PyTorch/infer/DeepSAD/subter_selective_test/inference/3_smoke_human_walking_2023-01-23.npy"
+)
+
+# scores_path = Path(
+# "/home/fedex/mt/projects/thesis-kowalczyk-jan/Deep-SAD-PyTorch/infer/DeepSAD/subter_split_test/inference/3_smoke_human_walking_2023-01-23.npy"
+# )
+# scores_path = Path(
+# "/home/fedex/mt/projects/thesis-kowalczyk-jan/Deep-SAD-PyTorch/infer/DeepSAD/subter_split_test/inference/1_loop_closure_illuminated_two_LiDARs_2023-01-23.npy"
+# )
+# scores_path = Path(
+# "/home/fedex/mt/projects/thesis-kowalczyk-jan/Deep-SAD-PyTorch/infer/DeepSAD/esmera_split_test/inference/experiment_3.npy"
+# )
+
+# dataset_path = Path(
+# "/home/fedex/mt/data/subter/1_loop_closure_illuminated_two_LiDARs_2023-01-23.bag"
+# )
+# dataset_path = Path("/home/fedex/mt/data/subter/3_smoke_human_walking_2023-01-23.bag")
+
+# dataset = load_dataset(dataset_path)
+fps = 10
+
+all_scores = np.load(scores_path)
+y_limit = 1.1 * np.max(all_scores[np.isfinite(all_scores)])
+# y_limit = 10
+# all_scores = np.where(all_scores > 10, 10, all_scores)
+# all_scores = all_scores.reshape(-1, 16).T #SUBTER
+# all_scores = all_scores.reshape(-1, 8).T  # ESMERA
+all_scores = all_scores.reshape(-1, 1).T  # ESMERA
+
+print(all_scores.shape, y_limit)
+
+fig, axes = plt.subplots(
+    1,
+    1,
+    figsize=(7.68, 7.2),
+    # 1, 1, figsize=(7.68, 7.2), gridspec_kw={"wspace": 0.10, "hspace": 0.10}
+)
+
+axes = [axes]
+# fig, axes = plt.subplots(
+# 1, 8, figsize=(20.48, 7.2), gridspec_kw={"wspace": 0.05, "hspace": 0.05}
+# )
+
+# Flatten the axes for easier indexing
+# axes = axes.flatten()
+
+last_running_avg = [None] * len(all_scores)
+
+
+# Function to calculate the running average with a dynamic window size
+def running_average_dynamic(data, current_frame, max_window=20):
+    """Calculate the running average with dynamic window size up to max_window."""
+    window = min(
+        current_frame, max_window
+    )  # Use the minimum of current frame or the max window size
+    return np.convolve(data, np.ones(window) / window, mode="valid")
+
+
+# Function to animate each subplot
+def animate(i, progress_bar, progress_task):
+    for score_id, ax in enumerate(axes):
+        ax.cla()  # Clear the axes for the current frame
+
+        ax.plot(
+            all_scores[score_id][:i], label="Anomaly Score"
+        )  # Plot the scores up to frame `i`
+        if i > 0:
+            avg_data = running_average_dynamic(
+                all_scores[score_id][:i], current_frame=i, max_window=20
+            )
+            ax.plot(range(len(avg_data)), avg_data, color="orange", label="Running Avg")
+
+        ax.set_xlim([0, all_scores[score_id].shape[0]])  # Set x limits
+        ax.set_ylim([0, y_limit])  # Set y limits
+        # ax.set_title(f"Score {score_id + 1}")  # Add a title for each subplot
+
+        ax.set_ylabel("Score", fontsize=10)  # Add y-axis label
+        ax.set_xlabel("Frame", fontsize=10)  # Add y-axis label
+        ax.legend(loc="upper right", fontsize=10)  # Add a legend
+        # Only show y-axis tick labels for the leftmost subplots
+        # if score_id % 8 == 0:  # First column of subplots
+        #     ax.set_ylabel("Score", fontsize=8)  # Add y-axis label
+        # else:
+        #     ax.set_yticklabels([])  # Remove y-axis tick labels
+        # if score_id < 8:
+        #     ax.set_xticklabels([])  # Remove x-axis tick labels
+
+        # ax.tick_params(labelsize=6)
+
+        # Update the running average text every 10 frames
+        if i % fps == 0 and i > 0:
+            # Calculate the current running average value (up to last 20 frames)
+            current_window = min(i, 20)
+            last_running_avg[score_id] = np.mean(
+                all_scores[score_id][i - current_window : i]
+            )
+
+        # Display the last updated running average value (if available)
+        if last_running_avg[score_id] is not None:
+            ax.text(
+                0.05,
+                0.95,
+                f"Current Avg: {last_running_avg[score_id]:>2.1f}",
+                transform=ax.transAxes,
+                fontsize=10,
+                verticalalignment="top",
+                horizontalalignment="left",
+                color="black",
+                fontfamily="monospace",
+            )
+    progress_bar.update(progress_task, completed=i)
+
+
+# plt.subplots_adjust(
+# left=0.02, right=0.98, top=0.95, bottom=0.05, wspace=0.05, hspace=0.05
+# )
+
+
+with Progress() as progress:
+    total = all_scores[0].shape[0] + 1
+    progress_task = progress.add_task("[cyan]Animating...", total=total)
+    progress.update(progress_task, completed=0)
+
+    animate_partial = partial(
+        animate, progress_bar=progress, progress_task=progress_task
+    )
+
+    # anim = animation.FuncAnimation(fig, animate_partial, frames=50, interval=1, blit=False)
+    anim = animation.FuncAnimation(
+        fig, animate_partial, frames=total, interval=1, blit=False
+    )
+
+    # Save the animation as a single video
+    animated_score_filename = f"{scores_path.stem}_selective.mp4"
+    anim.save(animated_score_filename, writer=animation.FFMpegWriter(fps=fps))
+
+    progress.update(progress_task, completed=all_scores[0].shape[0] + 1)
+
+# Clear the figure after saving
+fig.clear()

tools/anomaly_scatter_plot.py

@@ -0,0 +1,93 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Fix the random seed for reproducibility
+np.random.seed(0)
+
+# 1. Generate NORMAL DATA (e.g. points roughly around the origin)
+#    We'll keep them relatively close together so that a circle can enclose them easily.
+normal_data = np.random.randn(50, 2) * 0.75
+
+# 2. Generate ANOMALOUS DATA
+#    - Cluster 1: 3 points close together
+anomaly_cluster_1 = np.array([[3.0, 3.0], [3.2, 3.1], [2.8, 2.9], [0.4, 4.0]])
+
+#    - Cluster 2: A single point
+
+# 3. Compute the center and radius for a boundary circle around normal data
+center = normal_data.mean(axis=0)
+distances = np.linalg.norm(normal_data - center, axis=1)
+radius = (
+    np.max(distances) + 0.2
+)  # Add a small margin to ensure all normal points are inside
+
+# Create coordinates for plotting the circular boundary
+theta = np.linspace(0, 2 * np.pi, 200)
+circle_x = center[0] + radius * np.cos(theta)
+circle_y = center[1] + radius * np.sin(theta)
+
+# 4. Plot the data
+plt.figure(figsize=(7, 7))
+
+# Scatter normal points with 'o'
+plt.scatter(
+    normal_data[:, 0], normal_data[:, 1], marker="o", color="blue", label="Normal Data"
+)
+
+# Scatter anomalous points with 'x', but separate them by cluster for clarity
+plt.scatter(
+    anomaly_cluster_1[:, 0],
+    anomaly_cluster_1[:, 1],
+    marker="x",
+    color="red",
+    label="Anomalies",
+)
+
+# Plot the boundary (circle) around the normal data
+plt.plot(circle_x, circle_y, linestyle="--", color="black", label="Boundary")
+
+# 5. Annotate/label the clusters
+# Label the normal cluster near its center
+# plt.text(center[0], center[1],
+#          'Normal Cluster',
+#          horizontalalignment='center',
+#          verticalalignment='center',
+#          fontsize=9,
+#          bbox=dict(facecolor='white', alpha=0.7))
+#
+# # Label anomaly cluster 1 near its centroid
+# ac1_center = anomaly_cluster_1.mean(axis=0)
+# plt.text(ac1_center[0], ac1_center[1],
+#          'Anomaly Cluster 1',
+#          horizontalalignment='center',
+#          verticalalignment='center',
+#          fontsize=9,
+#          bbox=dict(facecolor='white', alpha=0.7))
+#
+# # Label anomaly cluster 2
+# ac2_point = anomaly_cluster_2[0]
+# plt.text(ac2_point[0]+0.2, ac2_point[1],
+#          'Anomaly Cluster 2',
+#          horizontalalignment='left',
+#          verticalalignment='center',
+#          fontsize=9,
+#          bbox=dict(facecolor='white', alpha=0.7))
+
+# Add legend and make plot look nice
+plt.legend(loc="upper left")
+# plt.title('2D Scatter Plot Showing Normal and Anomalous Clusters')
+plt.xlabel("x")
+plt.ylabel("y")
+plt.tick_params(
+    axis="both",  # changes apply to the x-axis
+    which="both",  # both major and minor ticks are affected
+    bottom=False,  # ticks along the bottom edge are off
+    top=False,  # ticks along the top edge are off
+    left=False,  # ticks along the top edge are off
+    right=False,  # ticks along the top edge are off
+    labelbottom=False,
+    labelleft=False,
+)  #
+# plt.grid(True)
+plt.axis("equal")  # Makes circles look circular rather than elliptical
+plt.savefig("scatter_plot.png")

tools/data_analyze.py

@@ -22,21 +22,6 @@ from rosbags.typesys.stores.ros1_noetic import (
 
 from util import existing_path
 
-
-def get_o3d_pointcloud(points: np.ndarray) -> o3d.geometry.PointCloud:
-    xyz_array = np.stack(
-        [
-            points["x"].astype(np.float64),
-            points["y"].astype(np.float64),
-            points["z"].astype(np.float64),
-        ],
-        axis=-1,
-    )
-    filtered_xyz = xyz_array[~np.all(xyz_array == 0, axis=1)]
-    o3d_vector = o3d.utility.Vector3dVector(filtered_xyz)
-    return o3d.geometry.PointCloud(o3d_vector)
-
-
 # Mapping of PointField datatypes to NumPy dtypes
 POINTFIELD_DATATYPES = {
     1: np.int8,  # INT8
@@ -178,9 +163,9 @@ def main() -> int:
     with AnyReader([args.input_experiment_path]) as reader:
         connections = reader.connections
         topics = dict(sorted({conn.topic: conn for conn in connections}.items()))
-        assert (
-            args.pointcloud_topic in topics
-        ), f"Topic {args.pointcloud_topic} not found"
+        assert args.pointcloud_topic in topics, (
+            f"Topic {args.pointcloud_topic} not found"
+        )
 
         original_types = {}
         typestore = get_typestore(Stores.ROS1_NOETIC)
@@ -207,16 +192,78 @@ def main() -> int:
                     "Processing all messages", total=reader.message_count
                 )
 
+                i = 0
                 for connection, timestamp, rawdata in reader.messages():
                     if connection.topic == args.pointcloud_topic:
+                        i += 1
                         # For the pointcloud topic, we need to modify the data
                         msg = reader.deserialize(rawdata, connection.msgtype)
                         original_pointcloud = read_pointcloud(msg)
                         cleaned_pointcloud = clean_pointcloud(original_pointcloud)
 
-                        o3d_pcd = get_o3d_pointcloud(cleaned_pointcloud)
+                        # o3d_pointcloud = o3d.geometry.PointCloud()
+
+                        # xyz_points = np.zeros((cleaned_pointcloud.shape[0], 3))
+                        # xyz_points[:, 0] = cleaned_pointcloud["x"]
+                        # xyz_points[:, 1] = cleaned_pointcloud["y"]
+                        # xyz_points[:, 2] = cleaned_pointcloud["z"]
+
+                        # o3d_pointcloud.points = o3d.utility.Vector3dVector(xyz_points)
+
+                        range_dtypes = {
+                            500: ("range_smaller_500", "u1"),
+                            800: ("range_smaller_800", "u1"),
+                            1000: ("range_smaller_1000", "u1"),
+                            1200: ("range_smaller_1200", "u1"),
+                        }
+
+                        # radius_outlier_dtypes = {
+                        #     # 10: ("radius_outlier_10", "u1"),
+                        #     50: ("radius_outlier_50", "u1"),
+                        #     # 100: ("radius_outlier_100", "u1"),
+                        # }
+
+                        # statistical_outlier_dtypes = {
+                        #     # (20, 1.0): ("statisstical_outlier_20_1", "u1"),
+                        #     # (20, 2.0): ("statisstical_outlier_20_2", "u1"),
+                        #     (20, 4.0): ("statisstical_outlier_20_4", "u1"),
+                        # }
+
+                        edited_dtype = np.dtype(
+                            cleaned_pointcloud.dtype.descr + list(range_dtypes.values())
+                            # + list(radius_outlier_dtypes.values())
+                            # + list(statistical_outlier_dtypes.values())
+                        )
+
+                        edited_pointcloud = np.zeros(
+                            cleaned_pointcloud.shape, dtype=edited_dtype
+                        )
+                        for name in cleaned_pointcloud.dtype.names:
+                            edited_pointcloud[name] = cleaned_pointcloud[name]
+
+                        for key, val in range_dtypes.items():
+                            edited_pointcloud[val[0]][
+                                edited_pointcloud["range"] < key
+                            ] = 255
+
+                        # for key, val in radius_outlier_dtypes.items():
+                        #     _, mask_indices = o3d_pointcloud.remove_radius_outlier(
+                        #         nb_points=2, radius=key
+                        #     )
+                        #     mask = np.zeros(edited_pointcloud.shape[0], dtype=bool)
+                        #     mask[mask_indices] = True
+                        #     edited_pointcloud[val[0]][mask] = 255
+
+                        # for key, val in statistical_outlier_dtypes.items():
+                        #     _, mask_indices = o3d_pointcloud.remove_statistical_outlier(
+                        #         nb_neighbors=key[0], std_ratio=key[1]
+                        #     )
+                        #     mask = np.zeros(edited_pointcloud.shape[0], dtype=bool)
+                        #     mask[mask_indices] = True
+                        #     edited_pointcloud[val[0]][mask] = 255
+
+                        msg = create_pointcloud2_msg(msg, edited_pointcloud)
 
-                        msg = create_pointcloud2_msg(msg, cleaned_pointcloud)
                         writer.write(
                             connection=new_connections[connection.id],
                             timestamp=timestamp,
@@ -224,6 +271,16 @@ def main() -> int:
                                 message=msg, typename=msg.__msgtype__
                             ),
                         )
+
+                        # Create a boolean mask where the condition is met
+                        # mask = cleaned_pointcloud["range"] > 1
+
+                        # Use the mask to set fields to zero
+                        # cleaned_pointcloud["x"][mask] = 0
+                        # cleaned_pointcloud["y"][mask] = 0
+                        # cleaned_pointcloud["z"][mask] = 0
+                        # cleaned_pointcloud["range"][mask] = 0
+
                     else:
                         # For all other topics, we can write rawdata directly, no need to deserialize
                         writer.write(new_connections[connection.id], timestamp, rawdata)

tools/evaluate_prc.py

@@ -0,0 +1,127 @@
+import pickle
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.stats import sem, t
+from sklearn.metrics import PrecisionRecallDisplay, auc
+
+
+def confidence_interval(data, confidence=0.95):
+    """Compute mean and margin of error for a given list of scores."""
+    n = len(data)
+    mean = np.mean(data)
+    # Standard error of the mean:
+    std_err = sem(data)
+    # Confidence interval radius
+    h = std_err * t.ppf((1 + confidence) / 2.0, n - 1)
+    return mean, h
+
+
+# 1) LOAD PRECISION-RECALL DATA
+prc_data = []  # Stores (precision, recall) for each DeepSAD fold
+ap_scores = []  # Average Precision for each DeepSAD fold
+
+isoforest_prc_data = []  # Stores (precision, recall) for each IsoForest fold
+isoforest_ap_scores = []  # Average Precision for each IsoForest fold
+
+results_path = Path(
+    "/home/fedex/mt/projects/thesis-kowalczyk-jan/Deep-SAD-PyTorch/log/DeepSAD/subter_kfold_800_3000_new"
+)
+
+# We assume we have 5 folds (adjust if you have a different number)
+for i in range(5):
+    with (results_path / f"results_{i}.pkl").open("rb") as f:
+        data = pickle.load(f)
+        precision, recall, _ = data["test_prc"]  # (precision, recall, thresholds)
+        prc_data.append((precision, recall))
+        # Compute Average Precision (AP) via AUC of the (recall, precision) curve
+        ap_scores.append(auc(recall, precision))
+
+    with (results_path / f"results_isoforest_{i}.pkl").open("rb") as f:
+        data = pickle.load(f)
+        precision, recall, _ = data["test_prc"]
+        isoforest_prc_data.append((precision, recall))
+        isoforest_ap_scores.append(auc(recall, precision))
+
+# 2) CALCULATE PER-FOLD STATISTICS
+mean_ap, ap_ci = confidence_interval(ap_scores)
+isoforest_mean_ap, isoforest_ap_ci = confidence_interval(isoforest_ap_scores)
+
+# 3) INTERPOLATE EACH FOLD'S PRC ON A COMMON RECALL GRID FOR MEAN CURVE
+mean_recall = np.linspace(0, 1, 100)
+
+# -- DeepSAD
+deep_sad_precisions_interp = []
+for precision, recall in prc_data:
+    # Interpolate precision values at mean_recall
+    interp_precision = np.interp(mean_recall, precision, recall)
+    deep_sad_precisions_interp.append(interp_precision)
+
+mean_precision = np.mean(deep_sad_precisions_interp, axis=0)
+std_precision = np.std(deep_sad_precisions_interp, axis=0)
+
+# -- IsoForest
+isoforest_precisions_interp = []
+for precision, recall in isoforest_prc_data:
+    interp_precision = np.interp(mean_recall, precision, recall)
+    isoforest_precisions_interp.append(interp_precision)
+
+isoforest_mean_precision = np.mean(isoforest_precisions_interp, axis=0)
+isoforest_std_precision = np.std(isoforest_precisions_interp, axis=0)
+
+# 4) CREATE PLOT USING PrecisionRecallDisplay
+fig, ax = plt.subplots(figsize=(8, 6))
+
+# (A) Plot each fold (optional) for DeepSAD
+for i, (precision, recall) in enumerate(prc_data):
+    disp = PrecisionRecallDisplay(precision=precision, recall=recall)
+    # Label only the first fold to avoid legend clutter
+    disp.plot(
+        ax=ax, color="b", alpha=0.3, label=f"DeepSAD Fold {i+1}" if i == 0 else None
+    )
+
+# (B) Plot each fold (optional) for IsoForest
+for i, (precision, recall) in enumerate(isoforest_prc_data):
+    disp = PrecisionRecallDisplay(precision=precision, recall=recall)
+    disp.plot(
+        ax=ax, color="r", alpha=0.3, label=f"IsoForest Fold {i+1}" if i == 0 else None
+    )
+
+# (C) Plot mean curve for DeepSAD
+mean_disp_deepsad = PrecisionRecallDisplay(precision=mean_precision, recall=mean_recall)
+mean_disp_deepsad.plot(
+    ax=ax, color="b", label=f"DeepSAD Mean PR (AP={mean_ap:.2f} ± {ap_ci:.2f})"
+)
+ax.fill_between(
+    mean_recall,
+    mean_precision - std_precision,
+    mean_precision + std_precision,
+    color="b",
+    alpha=0.2,
+)
+
+# (D) Plot mean curve for IsoForest
+mean_disp_isoforest = PrecisionRecallDisplay(
+    precision=isoforest_mean_precision, recall=mean_recall
+)
+mean_disp_isoforest.plot(
+    ax=ax,
+    color="r",
+    label=f"IsoForest Mean PR (AP={isoforest_mean_ap:.2f} ± {isoforest_ap_ci:.2f})",
+)
+ax.fill_between(
+    mean_recall,
+    isoforest_mean_precision - isoforest_std_precision,
+    isoforest_mean_precision + isoforest_std_precision,
+    color="r",
+    alpha=0.2,
+)
+
+# 5) FINAL PLOT ADJUSTMENTS
+ax.set_xlabel("Recall")
+ax.set_ylabel("Precision")
+ax.set_title("Precision-Recall Curve with 5-Fold Cross-Validation")
+ax.legend(loc="upper right")
+
+plt.savefig("pr_curve_800_3000_2.png")

tools/evaluate_prc_2.py

@@ -0,0 +1,135 @@
+import pickle
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.stats import sem, t
+from sklearn.metrics import auc
+
+
+# Confidence interval function
+def confidence_interval(data, confidence=0.95):
+    n = len(data)
+    mean = np.mean(data)
+    std_err = sem(data)
+    h = std_err * t.ppf((1 + confidence) / 2.0, n - 1)
+    return mean, h
+
+
+# Load PRC (precision-recall) data and compute AP (average precision)
+prc_data = []
+ap_scores = []
+isoforest_prc_data = []
+isoforest_ap_scores = []
+
+results_path = Path(
+    "/home/fedex/mt/projects/thesis-kowalczyk-jan/Deep-SAD-PyTorch/log/DeepSAD/subter_kfold_3000_800_2"
+)
+
+for i in range(5):
+    with (results_path / f"results_{i}.pkl").open("rb") as f:
+        data = pickle.load(f)
+        # data["test_prc"] should be (precision, recall, thresholds)
+        precision, recall, _ = data["test_prc"]
+        prc_data.append((precision, recall))
+
+        # Compute AP using area under the precision-recall curve
+        ap_scores.append(auc(recall, precision))
+
+    with (results_path / f"results_isoforest_{i}.pkl").open("rb") as f:
+        data = pickle.load(f)
+        precision, recall, _ = data["test_prc"]
+        isoforest_prc_data.append((precision, recall))
+
+        isoforest_ap_scores.append(auc(recall, precision))
+
+# Calculate mean and confidence interval for DeepSAD AP scores
+mean_ap, ap_ci = confidence_interval(ap_scores)
+
+# Interpolate precision over a common recall range for DeepSAD
+mean_recall = np.linspace(0, 1, 100)
+precisions = []
+
+for precision, recall in prc_data:
+    # Make sure recall is sorted (usually is from sklearn)
+    # Interpolate precision at the points in mean_recall
+    interp_prec = np.interp(mean_recall, np.flip(recall), np.flip(precision))
+    precisions.append(interp_prec)
+
+mean_precision = np.mean(precisions, axis=0)
+std_precision = np.std(precisions, axis=0)
+
+# Calculate mean and confidence interval for IsoForest AP scores
+isoforest_mean_ap, isoforest_ap_ci = confidence_interval(isoforest_ap_scores)
+
+# Interpolate precision over a common recall range for IsoForest
+isoforest_precisions = []
+for precision, recall in isoforest_prc_data:
+    interp_prec = np.interp(mean_recall, np.flip(recall), np.flip(precision))
+    isoforest_precisions.append(interp_prec)
+
+isoforest_mean_precision = np.mean(isoforest_precisions, axis=0)
+isoforest_std_precision = np.std(isoforest_precisions, axis=0)
+
+# Plot Precision-Recall curves with confidence margins
+plt.figure(figsize=(8, 6))
+
+# DeepSAD curve
+plt.plot(
+    mean_recall,
+    mean_precision,
+    color="b",
+    label=f"DeepSAD Mean PR (AP = {mean_ap:.2f} ± {ap_ci:.2f})",
+)
+plt.fill_between(
+    mean_recall,
+    mean_precision - std_precision,
+    mean_precision + std_precision,
+    color="b",
+    alpha=0.2,
+    label="DeepSAD ± 1 std. dev.",
+)
+
+# IsoForest curve
+plt.plot(
+    mean_recall,
+    isoforest_mean_precision,
+    color="r",
+    label=f"IsoForest Mean PR (AP = {isoforest_mean_ap:.2f} ± {isoforest_ap_ci:.2f})",
+)
+plt.fill_between(
+    mean_recall,
+    isoforest_mean_precision - isoforest_std_precision,
+    isoforest_mean_precision + isoforest_std_precision,
+    color="r",
+    alpha=0.2,
+    label="IsoForest ± 1 std. dev.",
+)
+
+# Optional: plot each fold's curve for DeepSAD
+for i, (precision, recall) in enumerate(prc_data):
+    plt.plot(
+        recall,
+        precision,
+        lw=1,
+        alpha=0.3,
+        color="b",
+        label=f"DeepSAD Fold {i + 1} PR" if i == 0 else "",
+    )
+
+# Optional: plot each fold's curve for IsoForest
+for i, (precision, recall) in enumerate(isoforest_prc_data):
+    plt.plot(
+        recall,
+        precision,
+        lw=1,
+        alpha=0.3,
+        color="r",
+        label=f"IsoForest Fold {i + 1} PR" if i == 0 else "",
+    )
+
+plt.xlabel("Recall")
+plt.ylabel("Precision")
+plt.title("Precision-Recall Curve with 5-Fold Cross-Validation")
+plt.legend(loc="upper right")
+plt.savefig("pr_curve_800_3000_4.png")

tools/evaluate_prc_singlerun.py

@@ -0,0 +1,50 @@
+import pickle
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+from sklearn.metrics import auc
+
+results_path = Path(
+    "/home/fedex/mt/projects/thesis-kowalczyk-jan/Deep-SAD-PyTorch/log/DeepSAD/subter_selective"
+)
+
+# Paths to your result files
+deepsad_results_file = results_path / "results.pkl"
+isoforest_results_file = results_path / "results_isoforest.pkl"
+
+# Load DeepSAD precision-recall data
+with deepsad_results_file.open("rb") as f:
+    data = pickle.load(f)
+    deep_precision, deep_recall, _ = data["test_prc"]
+
+# Compute AP for DeepSAD
+deep_ap = auc(deep_recall, deep_precision)
+
+# Load IsoForest precision-recall data
+with isoforest_results_file.open("rb") as f:
+    data = pickle.load(f)
+    iso_precision, iso_recall, _ = data["test_prc"]
+
+# Compute AP for IsoForest
+iso_ap = auc(iso_recall, iso_precision)
+
+# Create plot
+plt.figure(figsize=(8, 6))
+
+# Plot DeepSAD PR curve
+plt.plot(deep_recall, deep_precision, color="b", label=f"DeepSAD (AP = {deep_ap:.2f})")
+
+# Plot IsoForest PR curve
+plt.plot(iso_recall, iso_precision, color="r", label=f"IsoForest (AP = {iso_ap:.2f})")
+
+# Labels
+plt.xlabel("Recall")
+plt.ylabel("Precision")
+plt.title("Precision-Recall Curve (Single Run)")
+
+# Add legend
+plt.legend(loc="upper right")
+
+# Save and/or show plot
+plt.savefig("pr_curve_single_run.png")
+plt.show()

tools/evaluate_roc.py

@@ -0,0 +1,82 @@
+import pickle
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.stats import sem, t
+from sklearn.metrics import auc
+
+
+# Confidence interval function
+def confidence_interval(data, confidence=0.95):
+    n = len(data)
+    mean = np.mean(data)
+    std_err = sem(data)
+    h = std_err * t.ppf((1 + confidence) / 2.0, n - 1)
+    return mean, h
+
+
+# Load ROC and AUC values from pickle files
+roc_data = []
+auc_scores = []
+isoforest_roc_data = []
+isoforest_auc_scores = []
+
+results_path = Path(
+    "/home/fedex/mt/projects/thesis-kowalczyk-jan/Deep-SAD-PyTorch/log/DeepSAD/subter_kfold_0_0"
+)
+
+for i in range(5):
+    with (results_path / f"results_{i}.pkl").open("rb") as f:
+        data = pickle.load(f)
+        roc_data.append(data["test_roc"])
+        auc_scores.append(data["test_auc"])
+    with (results_path / f"results.isoforest_{i}.pkl").open("rb") as f:
+        data = pickle.load(f)
+        isoforest_roc_data.append(data["test_roc"])
+        isoforest_auc_scores.append(data["test_auc"])
+
+# Calculate mean and confidence interval for AUC scores
+mean_auc, auc_ci = confidence_interval(auc_scores)
+
+# Combine ROC curves
+mean_fpr = np.linspace(0, 1, 100)
+tprs = []
+
+for fpr, tpr, _ in roc_data:
+    interp_tpr = np.interp(mean_fpr, fpr, tpr)
+    interp_tpr[0] = 0.0
+    tprs.append(interp_tpr)
+
+mean_tpr = np.mean(tprs, axis=0)
+mean_tpr[-1] = 1.0
+std_tpr = np.std(tprs, axis=0)
+
+# Plot ROC curves with confidence margins
+plt.figure()
+plt.plot(
+    mean_fpr,
+    mean_tpr,
+    color="b",
+    label=f"Mean ROC (AUC = {mean_auc:.2f} ± {auc_ci:.2f})",
+)
+plt.fill_between(
+    mean_fpr,
+    mean_tpr - std_tpr,
+    mean_tpr + std_tpr,
+    color="b",
+    alpha=0.2,
+    label="± 1 std. dev.",
+)
+
+# Plot each fold's ROC curve (optional)
+ for i, (fpr, tpr, _) in enumerate(roc_data):
+     plt.plot(fpr, tpr, lw=1, alpha=0.3, label=f"Fold {i + 1} ROC")
+
+# Labels and legend
+plt.plot([0, 1], [0, 1], "k--", label="Chance")
+plt.xlabel("False Positive Rate")
+plt.ylabel("True Positive Rate")
+plt.title("ROC Curve with 5-Fold Cross-Validation")
+plt.legend(loc="lower right")
+plt.savefig("roc_curve_0_0.png")

tools/evaluate_roc_all.py

@@ -0,0 +1,133 @@
+import pickle
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.stats import sem, t
+from sklearn.metrics import auc
+
+
+# Confidence interval function
+def confidence_interval(data, confidence=0.95):
+    n = len(data)
+    mean = np.mean(data)
+    std_err = sem(data)
+    h = std_err * t.ppf((1 + confidence) / 2.0, n - 1)
+    return mean, h
+
+
+# Load ROC and AUC values from pickle files
+roc_data = []
+auc_scores = []
+isoforest_roc_data = []
+isoforest_auc_scores = []
+
+results_path = Path(
+    "/home/fedex/mt/projects/thesis-kowalczyk-jan/Deep-SAD-PyTorch/log/DeepSAD/subter_kfold_800_3000_new"
+)
+
+for i in range(5):
+    with (results_path / f"results_{i}.pkl").open("rb") as f:
+        data = pickle.load(f)
+        roc_data.append(data["test_roc"])
+        auc_scores.append(data["test_auc"])
+    with (results_path / f"results_isoforest_{i}.pkl").open("rb") as f:
+        data = pickle.load(f)
+        isoforest_roc_data.append(data["test_roc"])
+        isoforest_auc_scores.append(data["test_auc"])
+
+# Calculate mean and confidence interval for DeepSAD AUC scores
+mean_auc, auc_ci = confidence_interval(auc_scores)
+
+# Combine ROC curves for DeepSAD
+mean_fpr = np.linspace(0, 1, 100)
+tprs = []
+
+for fpr, tpr, _ in roc_data:
+    interp_tpr = np.interp(mean_fpr, fpr, tpr)
+    interp_tpr[0] = 0.0
+    tprs.append(interp_tpr)
+
+mean_tpr = np.mean(tprs, axis=0)
+mean_tpr[-1] = 1.0
+std_tpr = np.std(tprs, axis=0)
+
+# -- ADDED: Calculate mean and confidence interval for IsoForest AUC scores
+isoforest_mean_auc, isoforest_auc_ci = confidence_interval(isoforest_auc_scores)
+
+# -- ADDED: Combine ROC curves for IsoForest
+isoforest_mean_fpr = np.linspace(0, 1, 100)
+isoforest_tprs = []
+
+for fpr, tpr, _ in isoforest_roc_data:
+    interp_tpr = np.interp(isoforest_mean_fpr, fpr, tpr)
+    interp_tpr[0] = 0.0
+    isoforest_tprs.append(interp_tpr)
+
+isoforest_mean_tpr = np.mean(isoforest_tprs, axis=0)
+isoforest_mean_tpr[-1] = 1.0
+isoforest_std_tpr = np.std(isoforest_tprs, axis=0)
+
+# Plot ROC curves with confidence margins for DeepSAD
+plt.figure(figsize=(8, 6))
+plt.plot(
+    mean_fpr,
+    mean_tpr,
+    color="b",
+    label=f"DeepSAD Mean ROC (AUC = {mean_auc:.2f} ± {auc_ci:.2f})",
+)
+plt.fill_between(
+    mean_fpr,
+    mean_tpr - std_tpr,
+    mean_tpr + std_tpr,
+    color="b",
+    alpha=0.2,
+    label="DeepSAD ± 1 std. dev.",
+)
+
+# -- ADDED: Plot ROC curves with confidence margins for IsoForest
+plt.plot(
+    isoforest_mean_fpr,
+    isoforest_mean_tpr,
+    color="r",
+    label=f"IsoForest Mean ROC (AUC = {isoforest_mean_auc:.2f} ± {isoforest_auc_ci:.2f})",
+)
+plt.fill_between(
+    isoforest_mean_fpr,
+    isoforest_mean_tpr - isoforest_std_tpr,
+    isoforest_mean_tpr + isoforest_std_tpr,
+    color="r",
+    alpha=0.2,
+    label="IsoForest ± 1 std. dev.",
+)
+
+# Plot each fold's ROC curve (optional) for DeepSAD
+for i, (fpr, tpr, _) in enumerate(roc_data):
+    plt.plot(
+        fpr,
+        tpr,
+        lw=1,
+        alpha=0.3,
+        color="b",
+        label=f"DeepSAD Fold {i+1} ROC" if i == 0 else "",
+    )
+
+# -- ADDED: Plot each fold's ROC curve (optional) for IsoForest
+for i, (fpr, tpr, _) in enumerate(isoforest_roc_data):
+    plt.plot(
+        fpr,
+        tpr,
+        lw=1,
+        alpha=0.3,
+        color="r",
+        label=f"IsoForest Fold {i+1} ROC" if i == 0 else "",
+    )
+
+# Labels and legend
+plt.plot([0, 1], [0, 1], "k--", label="Chance")
+plt.xlabel("False Positive Rate")
+plt.ylabel("True Positive Rate")
+plt.title("ROC Curve with 5-Fold Cross-Validation")
+plt.legend(loc="lower right")
+plt.savefig("roc_curve_800_3000_isoforest.png")
+plt.show()

tools/list_topics.py

@@ -0,0 +1,37 @@
+from pathlib import Path
+
+import numpy as np
+from PIL import Image
+from rosbags.highlevel import AnyReader
+
+
+def list_topics(bag_file):
+    frame_count = 0
+    output_folder = Path("./test") / bag_file.stem
+    output_folder.mkdir(exist_ok=True, parents=True)
+
+    with AnyReader([bag_file]) as reader:
+        connections = reader.connections
+        topics = dict(sorted({conn.topic: conn for conn in connections}.items()))
+
+        topic = topics["/stereo_inertial_publisher/left/image_rect"]
+
+        for connection, timestamp, rawdata in reader.messages(connections=[topic]):
+            img_msg = reader.deserialize(rawdata, connection.msgtype)
+
+            img_data = np.frombuffer(img_msg.data, dtype=np.uint8).reshape(
+                (img_msg.height, img_msg.width)
+            )
+            img = Image.fromarray(img_data, mode="L")
+            frame_count += 1
+            img.save(output_folder / f"frame_{frame_count:04d}.png")
+
+        print(f"Saved {frame_count} frames to {output_folder}")
+
+
+if __name__ == "__main__":
+    bag_path = Path(
+        "/home/fedex/mt/data/subter/3_smoke_robot_stationary_static_excess_smoke_2023-01-23.bag"
+    )
+    print(f"{bag_path.exists()=}")
+    list_topics(bag_path)

tools/plot_score.py

@@ -0,0 +1,100 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+if __name__ == "__main__":
+    # Path to your .npy file with anomaly scores
+    scores_path = Path(
+        "/home/fedex/mt/projects/thesis-kowalczyk-jan/Deep-SAD-PyTorch/infer/"
+        "DeepSAD/subter_test/inference/3_smoke_human_walking_2023-01-23.npy"
+    )
+
+    # Frames per second
+    fps = 10
+
+    # Load all scores
+    all_scores = np.load(scores_path)
+
+    # (Optional) If your scores are a 1D array, remove or adjust this reshape
+    # The final shape is expected to be (1, total_frames) or (N,) if you want only one line.
+    all_scores = all_scores.reshape(1, -1)  # shape becomes (1, total_frames)
+
+    # We only want frames [0..299] => 300 frames
+    max_frames = 300
+    scores = all_scores[0, :max_frames]  # Take the first row's first 300 frames
+
+    # Convert frames to time in seconds: time[i] = i / fps
+    time_values = np.arange(len(scores)) / fps  # shape is (max_frames,)
+
+    # Recalculate y-limit based on the selected slice
+    y_limit = 1.1 * np.max(scores[np.isfinite(scores)])
+    # If you want to ensure y_limit >= 10, uncomment below:
+    # y_limit = max(y_limit, 10)
+
+    print("Selected Scores Shape:", scores.shape)
+    print("Time Range:", time_values[0], "->", time_values[-1])
+    print("Y-limit:", y_limit)
+
+    # --------------------------
+    # Improve default font sizes
+    # --------------------------
+    plt.rcParams.update(
+        {
+            "figure.dpi": 150,
+            "font.size": 14,
+            "axes.labelsize": 16,
+            "axes.titlesize": 16,
+            "xtick.labelsize": 14,
+            "ytick.labelsize": 14,
+            "legend.fontsize": 14,
+        }
+    )
+
+    # --------------------------------------------------
+    # Create the figure (half the previous width = 3.84)
+    # --------------------------------------------------
+    fig, ax = plt.subplots(figsize=(3.84, 7))
+
+    # Plot the scores vs. time
+    ax.plot(time_values, scores, label="Anomaly Score", color="blue")
+
+    # Set axis labels
+    ax.set_xlabel("Time [s]", fontsize=16)
+    ax.set_ylabel("Score", fontsize=16)
+
+    # Restrict X-axis from 0 to last time value (or 30 s if exactly 300 frames)
+    ax.set_xlim([0, time_values[-1]])
+    ax.set_ylim([0, y_limit])
+
+    # -------------------------
+    # Customize the ticks
+    # -------------------------
+    # Example: only 4 ticks on the x-axis (0, 10, 20, 30) if 300 frames => 30 s
+    ax.set_xticks([0, 10, 20, 30])
+
+    # Guarantee at least one tick at y=10.
+    # Adjust other ticks as you like; here we use 0 and the y_limit.
+    # If y_limit <= 10, you might want to override it.
+    if y_limit > 10:
+        ax.set_yticks([0, 10, round(y_limit, 1)])  # round for cleaner display
+    else:
+        # If your data doesn't go that high, just keep 0 and y_limit
+        ax.set_yticks([0, round(y_limit, 1)])
+
+    # Optionally add a legend
+    ax.legend(loc="upper right")
+
+    plt.axvline(x=22, color="r", label="shown image")
+
+    # Tight layout for small figure
+    plt.tight_layout()
+
+    # ----------------------
+    # Save to disk, no show
+    # ----------------------
+    output_filename = f"{scores_path.stem}_static_time_plot.png"
+    plt.savefig(output_filename, dpi=150)
+
+    # Clear the figure if you continue using matplotlib
+    plt.clf()

tools/plot_scripts/data_count_lidar_frames.py

@@ -0,0 +1,179 @@
+# this script loads the numpy array files containing the lidar frames and counts the number of frames in each file
+# the number of frames is then printed to the console per file as well as a congregated sum of all frames in files
+# containing the word smoke and ones that do not contain that word, as well as an overall sum of all frames
+# We also plot a pie chart of the distribution of data points in normal and anomalous experiments
+
+import shutil
+from datetime import datetime
+from pathlib import Path
+
+import numpy as np
+from tabulate import tabulate
+
+# define data path containing the numpy array files and output path for the plots
+data_path = Path("/home/fedex/mt/data/subter")
+
+output_path = Path("/home/fedex/mt/plots/data_count_lidar_frames")
+datetime_folder_name = datetime.now().strftime("%Y-%m-%d_%H-%M-%S")
+
+latest_folder_path = output_path / "latest"
+archive_folder_path = output_path / "archive"
+output_datetime_path = output_path / datetime_folder_name
+
+# if output does not exist, create it
+output_path.mkdir(exist_ok=True, parents=True)
+output_datetime_path.mkdir(exist_ok=True, parents=True)
+latest_folder_path.mkdir(exist_ok=True, parents=True)
+archive_folder_path.mkdir(exist_ok=True, parents=True)
+
+
+# find all numpy array files and sort them correctly by name
+normal_experiment_paths, anomaly_experiment_paths = [], []
+for npy_file_path in data_path.iterdir():
+    if npy_file_path.suffix != ".npy":
+        continue
+    if "smoke" in npy_file_path.name:
+        anomaly_experiment_paths.append(npy_file_path)
+    else:
+        normal_experiment_paths.append(npy_file_path)
+
+
+# function that counts the number of frames in one experiment
+def count_frames(npy_file_path):
+    frames = np.load(npy_file_path).shape[0]
+    return frames
+
+
+# we want to print the numbers of frames in a table so we first gather all the data in two maps
+normal_experiment_frames = {
+    npy_file_path.stem: count_frames(npy_file_path)
+    for npy_file_path in normal_experiment_paths
+}
+anomaly_experiment_frames = {
+    npy_file_path.stem: count_frames(npy_file_path)
+    for npy_file_path in anomaly_experiment_paths
+}
+
+# prepare data for tabulate
+normal_experiment_table = [
+    (experiment, frames) for experiment, frames in normal_experiment_frames.items()
+]
+anomaly_experiment_table = [
+    (experiment, frames) for experiment, frames in anomaly_experiment_frames.items()
+]
+
+# sort the tables by experiment name
+normal_experiment_table.sort(key=lambda x: x[0])
+anomaly_experiment_table.sort(key=lambda x: x[0])
+
+# add the sum of all frames to the tables
+normal_experiment_table.append(("Sum", sum(normal_experiment_frames.values())))
+anomaly_experiment_table.append(("Sum", sum(anomaly_experiment_frames.values())))
+
+# print the number of frames in each file using tabulate
+print("Normal experiments:")
+print(
+    tabulate(normal_experiment_table, headers=["Experiment", "Frames"], tablefmt="grid")
+)
+
+# print the smallest, largest, mean and median time of the normal experiments assuming 10 frames per second
+normal_experiment_frames_values = list(normal_experiment_frames.values())
+print(
+    f"Smallest time: {min(normal_experiment_frames_values) / 10} seconds, Largest time: {max(normal_experiment_frames_values) / 10} seconds, Mean time: {np.mean(normal_experiment_frames_values) / 10} seconds, Median time: {np.median(normal_experiment_frames_values) / 10} seconds"
+)
+
+print("Anomaly experiments:")
+print(
+    tabulate(
+        anomaly_experiment_table, headers=["Experiment", "Frames"], tablefmt="grid"
+    )
+)
+
+# print the smallest, largest, mean and median time of the anomalous experiments assuming 10 frames per second
+anomaly_experiment_frames_values = list(anomaly_experiment_frames.values())
+print(
+    f"Smallest time: {min(anomaly_experiment_frames_values) / 10} seconds, Largest time: {max(anomaly_experiment_frames_values) / 10} seconds, Mean time: {np.mean(anomaly_experiment_frames_values) / 10} seconds, Median time: {np.median(anomaly_experiment_frames_values) / 10} seconds"
+)
+
+# print the sum of all frames in all experiments
+total_frames = sum(normal_experiment_frames.values()) + sum(
+    anomaly_experiment_frames.values()
+)
+print(f"Total frames in all (normal and anmoaly) experiments: {total_frames} frames")
+
+# print the sum of normal and anomalous experiments as percentage of the total frames
+print(
+    f"Percentage of normal experiments: {sum(normal_experiment_frames.values()) / total_frames * 100}%"
+)
+print(
+    f"Percentage of anomaly experiments: {sum(anomaly_experiment_frames.values()) / total_frames * 100}%"
+)
+sum(normal_experiment_frames.values()) + sum(anomaly_experiment_frames.values())
+
+
+# define function to plot pie chart of the distribution of data points in normal and anomalous experiments
+def plot_data_points_pie(normal_experiment_frames, anomaly_experiment_frames):
+    import matplotlib.pyplot as plt
+
+    # we want to plot the sum of all frames in normal and anomaly experiments as total values and also percentages
+    total_normal_frames = sum(normal_experiment_frames.values())
+    total_anomaly_frames = sum(anomaly_experiment_frames.values())
+    total_frames = total_normal_frames + total_anomaly_frames
+
+    # prepare data for pie chart
+    labels = [
+        "Normal Lidar Frames\nNon-Degraded Pointclouds",
+        "Anomalous Lidar Frames\nDegraded Pointclouds",
+    ]
+    sizes = [total_normal_frames, total_anomaly_frames]
+    explode = (0.1, 0)  # explode the normal slice
+
+    # define an autopct function that shows percentage and total number of frames per slice
+    def make_autopct(pct):
+        return f"{pct:.1f}%\n({int(pct * total_frames / 100)} frames)"
+
+    # plot pie chart without shadow and with custom autopct
+    fig1, ax1 = plt.subplots()
+    # set a figure size of 10x5
+    fig1.set_size_inches(10, 5)
+    ax1.pie(sizes, explode=explode, labels=labels, autopct=make_autopct, shadow=False)
+    # for labels use center alignment
+
+    ax1.axis("equal")  # Equal aspect ratio ensures that pie is drawn as a circle.
+    # display the total number of frames in the center of the pie chart (adjusted vertically)
+    plt.text(
+        0,
+        0.2,
+        f"Total:\n{total_frames} frames",
+        fontsize=12,
+        ha="center",
+        va="center",
+        color="black",
+    )
+    plt.title(
+        "Distribution of Normal and Anomalous\nPointclouds in all Experiments (Lidar Frames)"
+    )
+    plt.tight_layout()
+
+    # save the plot
+    plt.savefig(output_datetime_path / "data_points_pie.png")
+
+
+plot_data_points_pie(normal_experiment_frames, anomaly_experiment_frames)
+
+# delete current latest folder
+shutil.rmtree(latest_folder_path, ignore_errors=True)
+
+# create new latest folder
+latest_folder_path.mkdir(exist_ok=True, parents=True)
+
+# copy contents of output folder to the latest folder
+for file in output_datetime_path.iterdir():
+    shutil.copy2(file, latest_folder_path)
+
+# copy this python script to the output datetime folder to preserve the code used to generate the plots
+shutil.copy2(__file__, output_datetime_path)
+shutil.copy2(__file__, latest_folder_path)
+
+# move output date folder to archive
+shutil.move(output_datetime_path, archive_folder_path)

tools/plot_scripts/data_missing_points.py

@@ -0,0 +1,243 @@
+import pickle
+import shutil
+from datetime import datetime
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+from pointcloudset import Dataset
+
+# define data path containing the bag files
+all_data_path = Path("/home/fedex/mt/data/subter")
+
+output_path = Path("/home/fedex/mt/plots/data_missing_points")
+datetime_folder_name = datetime.now().strftime("%Y-%m-%d_%H-%M-%S")
+
+latest_folder_path = output_path / "latest"
+archive_folder_path = output_path / "archive"
+output_datetime_path = output_path / datetime_folder_name
+
+# if output does not exist, create it
+output_path.mkdir(exist_ok=True, parents=True)
+output_datetime_path.mkdir(exist_ok=True, parents=True)
+latest_folder_path.mkdir(exist_ok=True, parents=True)
+archive_folder_path.mkdir(exist_ok=True, parents=True)
+
+data_resolution = 32 * 2048
+
+normal_experiment_paths, anomaly_experiment_paths = [], []
+
+# find all bag files and sort them correctly by name (experiments with smoke in the name are anomalies)
+for bag_file_path in all_data_path.iterdir():
+    if bag_file_path.suffix != ".bag":
+        continue
+    if "smoke" in bag_file_path.name:
+        anomaly_experiment_paths.append(bag_file_path)
+    else:
+        normal_experiment_paths.append(bag_file_path)
+
+# sort anomaly and normal experiments by filesize, ascending
+anomaly_experiment_paths.sort(key=lambda path: path.stat().st_size)
+normal_experiment_paths.sort(key=lambda path: path.stat().st_size)
+
+
+# function that plots histogram of how many points are missing in pointclouds for both normal and anomaly experiments
+def plot_data_points(normal_experiment_paths, anomaly_experiment_paths, title):
+    # function that finds the number of missing points in list of experiments (used for both normal and anomalous)
+    def find_missing_points(experiment_paths):
+        missing_points = []
+        for dataset in (
+            Dataset.from_file(experiment_path, topic="/ouster/points")
+            for experiment_path in experiment_paths
+        ):
+            missing_points_per_pc = []
+            for pc in dataset:
+                missing_points_per_pc.append(data_resolution - pc.data.shape[0])
+            missing_points.append(missing_points_per_pc)
+            # FIXME temporary break to test code on only one experiment
+            # break
+        return missing_points
+
+    # check if the data has already been calculated and saved to a pickle file
+    if (output_path / "missing_points.pkl").exists():
+        with open(output_path / "missing_points.pkl", "rb") as file:
+            missing_points_normal, missing_points_anomaly = pickle.load(file)
+    else:
+        missing_points_normal = find_missing_points(normal_experiment_paths)
+        missing_points_anomaly = find_missing_points(anomaly_experiment_paths)
+
+    # for faster subsequent runs, save the data to a pickle file
+    with open(output_path / "missing_points.pkl", "wb") as file:
+        pickle.dump(
+            (missing_points_normal, missing_points_anomaly),
+            file,
+            protocol=pickle.HIGHEST_PROTOCOL,
+        )
+
+    # combine all missing points into one list for each type of experiment
+    missing_points_normal = np.concatenate(missing_points_normal)
+    missing_points_anomaly = np.concatenate(missing_points_anomaly)
+
+    # create histogram of missing points for normal and anomaly experiments
+    plt.figure(figsize=(10, 5))
+    plt.hist(missing_points_normal, bins=100, alpha=0.5, label="Normal Experiments")
+    plt.hist(missing_points_anomaly, bins=100, alpha=0.5, label="Anomaly Experiments")
+    plt.title(title)
+    plt.xlabel("Number of Missing Points")
+    plt.ylabel("Number of Pointclouds")
+    plt.legend()
+    plt.tight_layout()
+    plt.savefig(output_datetime_path / "missing_points.png")
+
+    plt.clf()
+
+    # alternatively plot curves representing the data
+
+    # create alternative version with missing points on y axis and number of pointclouds on x axis
+    plt.figure(figsize=(10, 5))
+    plt.hist(
+        missing_points_normal,
+        bins=100,
+        alpha=0.5,
+        label="Normal Experiments",
+        orientation="horizontal",
+    )
+    plt.hist(
+        missing_points_anomaly,
+        bins=100,
+        alpha=0.5,
+        label="Anomaly Experiments",
+        orientation="horizontal",
+    )
+    plt.title(title)
+    plt.xlabel("Number of Pointclouds")
+    plt.ylabel("Number of Missing Points")
+    plt.legend()
+    plt.tight_layout()
+    plt.savefig(output_datetime_path / "missing_points_alternative.png")
+
+    # find min and max of both categories so we can set the same limits for both plots
+    min = np.min([np.min(missing_points_normal), np.min(missing_points_anomaly)])
+    max = np.max([np.max(missing_points_normal), np.max(missing_points_anomaly)])
+
+    # create bins array with min and max values
+    bins = np.linspace(min, max, 100)
+
+    # since the two histograms (normal and anomalous) have different scales, normalize their amplitude and plot a density version as well
+    plt.clf()
+    plt.figure(figsize=(10, 5))
+    plt.hist(
+        missing_points_normal,
+        bins=bins,
+        alpha=0.5,
+        label="Normal Experiments",
+        color="blue",
+        density=True,
+    )
+    plt.hist(
+        missing_points_anomaly,
+        bins=bins,
+        alpha=0.5,
+        color="red",
+        label="Anomaly Experiments",
+        density=True,
+    )
+    plt.title(title)
+    plt.xlabel("Number of Missing Points")
+    plt.ylabel("Density")
+    plt.legend()
+    plt.tight_layout()
+    plt.savefig(output_datetime_path / "missing_points_density.png")
+
+    # create another density version which does not plot number of missing points but percentage of measurements that are missing (total number of points is 32*2048)
+    bins = np.linspace(0, 1, 100)
+    plt.clf()
+    plt.figure(figsize=(10, 5))
+    plt.hist(
+        missing_points_normal / data_resolution,
+        bins=bins,
+        alpha=0.5,
+        label="Normal Experiments (No Artifical Smoke)",
+        color="blue",
+        density=True,
+    )
+    plt.hist(
+        missing_points_anomaly / data_resolution,
+        bins=bins,
+        alpha=0.5,
+        color="red",
+        label="Anomaly Experiments (With Artifical Smoke)",
+        density=True,
+    )
+    plt.title(title)
+    plt.xlabel("Percentage of Missing Lidar Measurements")
+    plt.ylabel("Density")
+    # display the x axis as percentages
+    plt.gca().set_xticklabels(
+        ["{:.0f}%".format(x * 100) for x in plt.gca().get_xticks()]
+    )
+    plt.legend()
+    plt.tight_layout()
+    plt.savefig(output_datetime_path / "missing_points_density_percentage.png")
+
+    # mathplotlib does not support normalizing the histograms to the same scale, so we do it manually using numpy
+    num_bins = 100
+    bin_lims = np.linspace(0, 40000, num_bins + 1)
+    bin_centers = 0.5 * (bin_lims[:-1] + bin_lims[1:])
+    bin_widths = bin_lims[1:] - bin_lims[:-1]
+    # calculate the histogram for normal and anomaly experiments
+    normal_hist, _ = np.histogram(missing_points_normal, bins=bin_lims)
+    anomaly_hist, _ = np.histogram(missing_points_anomaly, bins=bin_lims)
+    # normalize the histograms to the same scale
+    normal_hist_normalized = np.array(normal_hist, dtype=float) / np.max(normal_hist)
+    anomaly_hist_normalized = np.array(anomaly_hist, dtype=float) / np.max(anomaly_hist)
+    # plot the normalized histograms
+    plt.clf()
+    plt.figure(figsize=(10, 5))
+    plt.bar(
+        bin_centers,
+        normal_hist_normalized,
+        width=bin_widths,
+        align="center",
+        alpha=0.5,
+        label="Normal Experiments",
+    )
+    plt.bar(
+        bin_centers,
+        anomaly_hist_normalized,
+        width=bin_widths,
+        align="center",
+        alpha=0.5,
+        label="Anomaly Experiments",
+    )
+    plt.title(title)
+    plt.xlabel("Number of Missing Points")
+    plt.ylabel("Normalized Density")
+    plt.legend()
+    plt.tight_layout()
+    plt.savefig(output_datetime_path / "missing_points_normalized.png")
+
+
+# plot histogram of missing points for normal and anomaly experiments
+plot_data_points(
+    normal_experiment_paths,
+    anomaly_experiment_paths,
+    "Missing Lidar Measurements per Scan",
+)
+
+# delete current latest folder
+shutil.rmtree(latest_folder_path, ignore_errors=True)
+
+# create new latest folder
+latest_folder_path.mkdir(exist_ok=True, parents=True)
+
+# copy contents of output folder to the latest folder
+for file in output_datetime_path.iterdir():
+    shutil.copy2(file, latest_folder_path)
+
+# copy this python script to the output datetime folder to preserve the code used to generate the plots
+shutil.copy2(__file__, output_datetime_path)
+shutil.copy2(__file__, latest_folder_path)
+
+# move output date folder to archive
+shutil.move(output_datetime_path, archive_folder_path)

tools/plot_scripts/data_particles_near_sensor.py

@@ -0,0 +1,220 @@
+import pickle
+import shutil
+from datetime import datetime
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+from pointcloudset import Dataset
+
+# define data path containing the bag files
+all_data_path = Path("/home/fedex/mt/data/subter")
+
+output_path = Path("/home/fedex/mt/plots/data_particles_near_sensor")
+datetime_folder_name = datetime.now().strftime("%Y-%m-%d_%H-%M-%S")
+
+latest_folder_path = output_path / "latest"
+archive_folder_path = output_path / "archive"
+output_datetime_path = output_path / datetime_folder_name
+
+# if output does not exist, create it
+output_path.mkdir(exist_ok=True, parents=True)
+output_datetime_path.mkdir(exist_ok=True, parents=True)
+latest_folder_path.mkdir(exist_ok=True, parents=True)
+archive_folder_path.mkdir(exist_ok=True, parents=True)
+
+normal_experiment_paths, anomaly_experiment_paths = [], []
+
+# find all bag files and sort them correctly by name (experiments with smoke in the name are anomalies)
+for bag_file_path in all_data_path.iterdir():
+    if bag_file_path.suffix != ".bag":
+        continue
+    if "smoke" in bag_file_path.name:
+        anomaly_experiment_paths.append(bag_file_path)
+    else:
+        normal_experiment_paths.append(bag_file_path)
+
+# print out the names of the normal and anomaly experiments that we found
+print("Normal experiments:")
+for path in normal_experiment_paths:
+    print(path.name)
+print("\nAnomaly experiments:")
+for path in anomaly_experiment_paths:
+    print(path.name)
+
+# sort anomaly and normal experiments by filesize, ascending
+anomaly_experiment_paths.sort(key=lambda path: path.stat().st_size)
+normal_experiment_paths.sort(key=lambda path: path.stat().st_size)
+
+
+# function that plots histogram of how many measurements have a range smaller than 0.5m for both normal and anomaly experiments
+def plot_data_points(normal_experiment_paths, anomaly_experiment_paths, title):
+    # function that finds the number of measurements with a range smaller than 0.5m in list of experiments (used for both normal and anomalous)
+    def find_particles_near_sensor(experiment_paths, range_threshold):
+        particles_near_sensor = []
+        for dataset in (
+            Dataset.from_file(experiment_path, topic="/ouster/points")
+            for experiment_path in experiment_paths
+        ):
+            particles_near_sensor_per_pc = []
+            for pc in dataset:
+                particles_near_sensor_per_pc.append(
+                    pc.data[pc.data["range"] < range_threshold].shape[0]
+                )
+            particles_near_sensor.append(particles_near_sensor_per_pc)
+        return particles_near_sensor
+
+    range_thresholds = [500, 750, 1000, 1250, 1500]
+
+    for rt in range_thresholds:
+        print(f"Processing range threshold {rt}...")
+
+        # check if the data has already been calculated and saved to a pickle file
+        if (output_path / f"particles_near_sensor_counts_{rt}.pkl").exists():
+            with open(
+                output_path / f"particles_near_sensor_counts_{rt}.pkl", "rb"
+            ) as file:
+                particles_near_sensor_normal, particles_near_sensor_anomaly = (
+                    pickle.load(file)
+                )
+        else:
+            particles_near_sensor_normal = find_particles_near_sensor(
+                normal_experiment_paths,
+                rt,
+            )
+            particles_near_sensor_anomaly = find_particles_near_sensor(
+                anomaly_experiment_paths,
+                rt,
+            )
+
+        # for faster subsequent runs, save the data to a pickle file
+        with open(output_path / f"particles_near_sensor_counts_{rt}.pkl", "wb") as file:
+            pickle.dump(
+                (particles_near_sensor_normal, particles_near_sensor_anomaly),
+                file,
+                protocol=pickle.HIGHEST_PROTOCOL,
+            )
+
+        # combine all counts of how many particles are close to the sensor into one list for each type of experiment
+        particles_near_sensor_normal = np.concatenate(particles_near_sensor_normal)
+        particles_near_sensor_anomaly = np.concatenate(particles_near_sensor_anomaly)
+
+        # find min and max of both categories so we can set the same limits for both plots
+        min = np.min(
+            [
+                np.min(particles_near_sensor_normal),
+                np.min(particles_near_sensor_anomaly),
+            ]
+        )
+        max = np.max(
+            [
+                np.max(particles_near_sensor_normal),
+                np.max(particles_near_sensor_anomaly),
+            ]
+        )
+
+        # create bins array with min and max values
+        bins = np.linspace(min, max, 100)
+
+        # since the two histograms (normal and anomalous) have different scales, normalize their amplitude and plot a density version as well
+        # commented out since boxplot is more informative
+        # plt.clf()
+        # plt.figure(figsize=(10, 5))
+        # plt.hist(
+        #     particles_near_sensor_normal,
+        #     bins=bins,
+        #     alpha=0.5,
+        #     label="Normal Experiments (No Artifical Smoke)",
+        #     color="blue",
+        #     density=True,
+        # )
+        # plt.hist(
+        #     particles_near_sensor_anomaly,
+        #     bins=bins,
+        #     alpha=0.5,
+        #     color="red",
+        #     label="Anomaly Experiments (With Artifical Smoke)",
+        #     density=True,
+        # )
+        # plt.title(title)
+        # plt.xlabel("Number of Particles Near Sensor")
+        # plt.ylabel("Density")
+        # plt.legend()
+        # plt.tight_layout()
+        # plt.savefig(output_datetime_path / f"particles_near_sensor_density_{rt}.png")
+
+        # alternatively create a box plot to show the distribution of the data
+        # instead of plotting the frequency of particles near sensor we'll plot the percentage of points (compared to the total number of points in the pointcloud)
+
+        data_resolution = 32 * 2048
+
+        plt.clf()
+        plt.figure(figsize=(10, 5))
+        plt.boxplot(
+            [
+                particles_near_sensor_normal / data_resolution,
+                particles_near_sensor_anomaly / data_resolution,
+            ],
+            tick_labels=[
+                "Normal Experiments (No Artifical Smoke)",
+                "Anomaly Experiments (With Artifical Smoke)",
+            ],
+        )
+        # format the y axis labels as percentages
+        plt.gca().set_yticklabels(
+            ["{:.0f}%".format(y * 100) for y in plt.gca().get_yticks()]
+        )
+        plt.title("Particles Closer than 0.5m to the Sensor")
+        plt.ylabel("Percentage of measurements closer than 0.5m")
+        plt.tight_layout()
+        plt.savefig(output_datetime_path / f"particles_near_sensor_boxplot_{rt}.png")
+
+        # we create the same boxplot but limit the y-axis to 5% to better see the distribution of the data
+        plt.clf()
+        plt.figure(figsize=(10, 5))
+        plt.boxplot(
+            [
+                particles_near_sensor_normal / data_resolution,
+                particles_near_sensor_anomaly / data_resolution,
+            ],
+            tick_labels=[
+                "Normal Experiments (No Artifical Smoke)",
+                "Anomaly Experiments (With Artifical Smoke)",
+            ],
+        )
+        # format the y axis labels as percentages
+        plt.gca().set_yticklabels(
+            ["{:.0f}%".format(y * 100) for y in plt.gca().get_yticks()]
+        )
+        plt.title("Particles Closer than 0.5m to the Sensor")
+        plt.ylabel("Percentage of measurements closer than 0.5m")
+        plt.ylim(0, 0.05)
+        plt.tight_layout()
+        plt.savefig(
+            output_datetime_path / f"particles_near_sensor_boxplot_zoomed_{rt}.png"
+        )
+
+
+# plot histogram of how many measurements have a range smaller than 0.5m for both normal and anomaly experiments
+plot_data_points(
+    normal_experiment_paths,
+    anomaly_experiment_paths,
+    "Density of Number of Particles Near Sensor",
+)
+
+# delete current latest folder
+shutil.rmtree(latest_folder_path, ignore_errors=True)
+
+# create new latest folder
+latest_folder_path.mkdir(exist_ok=True, parents=True)
+
+# copy contents of output folder to the latest folder
+for file in output_datetime_path.iterdir():
+    shutil.copy2(file, latest_folder_path)
+
+# copy this python script to the output datetime folder to preserve the code used to generate the plots
+shutil.copy2(__file__, output_datetime_path)
+shutil.copy2(__file__, latest_folder_path)
+
+# move output date folder to archive
+shutil.move(output_datetime_path, archive_folder_path)

tools/plot_scripts/data_spherical_projection.py

@@ -0,0 +1,188 @@
+import argparse
+import shutil
+from datetime import datetime
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib import colormaps
+from matplotlib.cm import ScalarMappable
+from matplotlib.colors import Colormap, ListedColormap
+from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable
+from PIL import Image  # new import
+
+
+def get_colormap_with_special_missing_color(
+    colormap_name: str, missing_data_color: str = "black", reverse: bool = False
+) -> Colormap:
+    colormap = (
+        colormaps[colormap_name] if not reverse else colormaps[f"{colormap_name}_r"]
+    )
+    colormap.set_bad(missing_data_color)
+    colormap.set_over("white")
+    return colormap
+
+
+# --- Setup output folders (similar to data_missing_points.py) ---
+# Change the output path as needed
+output_path = Path("/home/fedex/mt/plots/data_2d_projections")
+datetime_folder_name = datetime.now().strftime("%Y-%m-%d_%H-%M-%S")
+output_datetime_path = output_path / datetime_folder_name
+latest_folder_path = output_path / "latest"
+archive_folder_path = output_path / "archive"
+
+for folder in (
+    output_path,
+    output_datetime_path,
+    latest_folder_path,
+    archive_folder_path,
+):
+    folder.mkdir(exist_ok=True, parents=True)
+
+# --- Parse command-line arguments ---
+parser = argparse.ArgumentParser(
+    description="Plot two 2D projections (from .npy files) in vertical subplots"
+)
+parser.add_argument(
+    "--input1",
+    type=Path,
+    default=Path(
+        "/home/fedex/mt/data/subter/new_projection/1_loop_closure_illuminated_2023-01-23.npy"
+    ),
+    help="Path to first .npy file containing 2D projection data",
+)
+parser.add_argument(
+    "--input2",
+    type=Path,
+    default=Path(
+        "/home/fedex/mt/data/subter/new_projection/3_smoke_human_walking_2023-01-23.npy"
+    ),
+    help="Path to second .npy file containing 2D projection data",
+)
+parser.add_argument(
+    "--frame1",
+    type=int,
+    default=955,
+    help="Frame index to plot from first file (0-indexed)",
+)
+parser.add_argument(
+    "--frame2",
+    type=int,
+    default=242,
+    help="Frame index to plot from second file (0-indexed)",
+)
+parser.add_argument(
+    "--colormap",
+    default="viridis",
+    type=str,
+    help="Name of matplotlib colormap to use",
+)
+parser.add_argument(
+    "--missing-data-color",
+    default="black",
+    type=str,
+    help="Color to use for missing data in projection",
+)
+parser.add_argument(
+    "--reverse-colormap",
+    action="store_true",
+    help="Reverse the colormap if specified",
+)
+args = parser.parse_args()
+
+# --- Load the numpy projection data from the provided files ---
+# Each file is assumed to be a 3D array: (num_frames, height, width)
+proj_data1 = np.load(args.input1)
+proj_data2 = np.load(args.input2)
+
+# Choose the desired frames
+try:
+    frame1 = proj_data1[args.frame1]
+except IndexError:
+    raise ValueError(f"Frame index {args.frame1} out of range for file: {args.input1}")
+try:
+    frame2 = proj_data2[args.frame2]
+except IndexError:
+    raise ValueError(f"Frame index {args.frame2} out of range for file: {args.input2}")
+
+# Determine shared range across both frames (ignoring NaNs)
+global_vmin = 0  # min(np.nanmin(frame1), np.nanmin(frame2))
+global_vmax = 0.8  # max(np.nanmax(frame1), np.nanmax(frame2))
+
+# Create colormap using the utility (to mimic create_2d_projection)
+cmap = get_colormap_with_special_missing_color(
+    args.colormap, args.missing_data_color, args.reverse_colormap
+)
+
+# --- Create a figure with 2 vertical subplots ---
+fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(10, 5))
+for ax, frame, title in zip(
+    (ax1, ax2),
+    (frame1, frame2),
+    (
+        "Projection of Lidar Frame without Degradation",
+        "Projection of Lidar Frame with Degradation (Artifical Smoke)",
+    ),
+):
+    im = ax.imshow(frame, cmap=cmap, aspect="auto", vmin=global_vmin, vmax=global_vmax)
+    ax.set_title(title)
+    ax.axis("off")
+
+# Adjust layout to fit margins for a paper
+plt.tight_layout(rect=[0, 0.05, 1, 1])
+# Add a colorbar with the colormap below the subplots
+cbar = fig.colorbar(im, ax=[ax1, ax2], orientation="vertical", fraction=0.05)
+cbar.set_label("Normalized Range")
+
+# Add a separate colorbar for NaN values
+sm = ScalarMappable(cmap=ListedColormap([cmap.get_bad(), cmap.get_over()]))
+divider = make_axes_locatable(cbar.ax)
+nan_ax = divider.append_axes(
+    "bottom", size="15%", pad="3%", aspect=3, anchor=cbar.ax.get_anchor()
+)
+nan_ax.grid(visible=False, which="both", axis="both")
+nan_cbar = fig.colorbar(sm, cax=nan_ax, orientation="vertical")
+nan_cbar.set_ticks([0.3, 0.7])
+nan_cbar.set_ticklabels(["NaN", f"> {global_vmax}"])
+nan_cbar.ax.tick_params(length=0)
+
+# Save the combined plot
+output_file = output_datetime_path / "data_2d_projections.png"
+plt.savefig(output_file, dpi=300, bbox_inches="tight", pad_inches=0.1)
+plt.close()
+
+print(f"Plot saved to: {output_file}")
+
+# --- Create grayscale images (high precision) from the numpy frames using Pillow ---
+# Convert NaN values to 0 and ensure the array is in float32 for 32-bit precision.
+for degradation_status, frame_number, frame in (
+    ("normal", args.frame1, frame1),
+    ("smoke", args.frame2, frame2),
+):
+    frame_gray = np.nan_to_num(frame, nan=0).astype(np.float32)
+    gray_image = Image.fromarray(frame_gray, mode="F")
+    gray_output_file = (
+        output_datetime_path
+        / f"frame_{frame_number}_grayscale_{degradation_status}.tiff"
+    )
+    gray_image.save(gray_output_file)
+    print(f"Grayscale image saved to: {gray_output_file}")
+
+# --- Handle folder structure: update latest folder and archive the output folder ---
+
+# Delete current latest folder and recreate it
+shutil.rmtree(latest_folder_path, ignore_errors=True)
+latest_folder_path.mkdir(exist_ok=True, parents=True)
+
+# Copy contents of the current output datetime folder to latest folder
+for file in output_datetime_path.iterdir():
+    shutil.copy2(file, latest_folder_path)
+
+# Copy this python script to both output datetime folder and latest folder for preservation
+script_path = Path(__file__)
+shutil.copy2(script_path, output_datetime_path)
+shutil.copy2(script_path, latest_folder_path)
+
+# Move the output datetime folder to the archive folder
+shutil.move(output_datetime_path, archive_folder_path)
+print(f"Output archived to: {archive_folder_path}")

tools/plot_scripts/deepsad_tsne_latent_space.py

@@ -0,0 +1,160 @@
+import shutil
+from datetime import datetime
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+from sklearn.decomposition import PCA
+from sklearn.manifold import TSNE
+
+# define data path containing the npy output files from your method
+all_data_path = Path(
+    "/home/fedex/mt/projects/thesis-kowalczyk-jan/Deep-SAD-PyTorch/infer/DeepSAD/all_infer/inference"
+)
+
+output_path = Path("/home/fedex/mt/plots/deepsad_reduced_latent_space")
+datetime_folder_name = datetime.now().strftime("%Y-%m-%d_%H-%M-%S")
+
+latest_folder_path = output_path / "latest"
+archive_folder_path = output_path / "archive"
+output_datetime_path = output_path / datetime_folder_name
+
+# create required output directories if they do not exist
+output_path.mkdir(exist_ok=True, parents=True)
+output_datetime_path.mkdir(exist_ok=True, parents=True)
+latest_folder_path.mkdir(exist_ok=True, parents=True)
+archive_folder_path.mkdir(exist_ok=True, parents=True)
+
+normal_experiment_paths = []
+anomaly_experiment_paths = []
+
+# locate and sort the npy files (experiment outputs) based on file size
+for file_path in all_data_path.iterdir():
+    if file_path.suffix != ".npy":
+        continue
+    # check if the file name contains "output" to ensure it's an experiment output file
+    if "output" not in file_path.stem:
+        continue
+    if "smoke" in file_path.name:
+        anomaly_experiment_paths.append(file_path)
+    else:
+        normal_experiment_paths.append(file_path)
+
+print("Normal experiments:")
+for path in normal_experiment_paths:
+    print(path.name)
+
+print("\nAnomaly experiments:")
+for path in anomaly_experiment_paths:
+    print(path.name)
+
+normal_experiment_paths.sort(key=lambda path: path.stat().st_size)
+anomaly_experiment_paths.sort(key=lambda path: path.stat().st_size)
+
+
+def load_latent_space_data(experiment_paths):
+    """
+    Load latent space data from npy files and return a single numpy array.
+    Modify this function if your file structure is different.
+    """
+    data_list = []
+    for path in experiment_paths:
+        latent_data = np.load(path)
+        data_list.append(latent_data)
+    return np.vstack(data_list)
+
+
+def reduce_dimensionality(data, n_components=50):
+    """
+    Reduce the dimensionality of the data using PCA.
+    This function can be re-used by TSNE or other methods for an initial reduction.
+    """
+    pca = PCA(n_components=n_components, random_state=42)
+    return pca.fit_transform(data)
+
+
+def plot_tsne_latent_space(normal_data, anomaly_data, title="TSNE of Latent Space"):
+    """
+    Plot the TSNE representation of the latent space.
+    This function first applies a PCA-based dimensionality reduction for efficiency.
+    """
+    # Combine normal and anomaly data
+    combined_data = np.vstack((normal_data, anomaly_data))
+
+    # Initial dimensionality reduction with PCA
+    reduced_data = reduce_dimensionality(combined_data, n_components=50)
+
+    # Apply TSNE transformation on the PCA-reduced data
+    tsne = TSNE(n_components=2, random_state=42)
+    tsne_results = tsne.fit_transform(reduced_data)
+
+    # Split the TSNE results back into normal and anomaly arrays
+    tsne_normal = tsne_results[: len(normal_data)]
+    tsne_anomaly = tsne_results[len(normal_data) :]
+
+    # Plotting TSNE results
+    plt.clf()
+    plt.figure(figsize=(10, 5))
+    plt.scatter(
+        tsne_anomaly[:, 0], tsne_anomaly[:, 1], label="Anomaly", alpha=0.6, marker="x"
+    )
+    plt.scatter(
+        tsne_normal[:, 0], tsne_normal[:, 1], label="Normal", alpha=0.6, marker="o"
+    )
+    plt.title(title)
+    plt.legend()
+    plt.tight_layout()
+    plt.savefig(output_datetime_path / "tsne_latent_space_plot.png")
+
+
+def plot_pca_scatter(normal_data, anomaly_data, title="PCA Scatter Plot"):
+    """
+    Plot a 2-dimensional scatterplot of the latent space using PCA.
+    This is useful for visualization and can be easily extended.
+    """
+    # Combine normal and anomaly data
+    combined_data = np.vstack((normal_data, anomaly_data))
+    pca = PCA(n_components=2, random_state=42)
+    pca_results = pca.fit_transform(combined_data)
+
+    # Split the PCA results back into normal and anomaly arrays
+    pca_normal = pca_results[: len(normal_data)]
+    pca_anomaly = pca_results[len(normal_data) :]
+
+    # Plotting PCA scatter results
+    plt.clf()
+    plt.figure(figsize=(10, 5))
+    plt.scatter(
+        pca_anomaly[:, 0], pca_anomaly[:, 1], label="Anomaly", alpha=0.6, marker="x"
+    )
+    plt.scatter(
+        pca_normal[:, 0], pca_normal[:, 1], label="Normal", alpha=0.6, marker="o"
+    )
+    plt.title(title)
+    plt.legend()
+    plt.tight_layout()
+    plt.savefig(output_datetime_path / "pca_latent_space_plot.png")
+
+
+# load latent space data for both normal and anomalous experiments
+normal_data = load_latent_space_data(normal_experiment_paths)
+anomaly_data = load_latent_space_data(anomaly_experiment_paths)
+
+# call your plotting functions
+plot_tsne_latent_space(normal_data, anomaly_data)
+plot_pca_scatter(normal_data, anomaly_data)
+
+# update the 'latest' results folder: delete previous and copy current outputs
+shutil.rmtree(latest_folder_path, ignore_errors=True)
+latest_folder_path.mkdir(exist_ok=True, parents=True)
+
+for file in output_datetime_path.iterdir():
+    shutil.copy2(file, latest_folder_path)
+
+# copy this script to the output folder and to the latest folder to preserve the used code
+script_path = Path(__file__)
+shutil.copy2(script_path, output_datetime_path)
+shutil.copy2(script_path, latest_folder_path)
+
+# move the output date folder to the archive folder for record keeping
+shutil.move(output_datetime_path, archive_folder_path)

tools/pyproject.toml

@@ -20,6 +20,7 @@ dask = "^2024.4.2"
 dask-expr = "^1.1.3"
 pandas = "^2.2.2"
 pathvalidate = "^3.2.0"
+tabulate = "^0.9.0"
 
 [build-system]
 requires = ["poetry-core"]

tools/render2d.py

@@ -137,47 +137,47 @@ def process_frame(args) -> Tuple[int, np.ndarray, Optional[Path]]:
     lidar_data["horizontal_position"] = (
         lidar_data["original_id"] % horizontal_resolution
     )
-    lidar_data["horizontal_position_yaw_f"] = (
-        0.5
-        * horizontal_resolution
-        * (np.arctan2(lidar_data["y"], lidar_data["x"]) / pi + 1.0)
-    )
-    lidar_data["horizontal_position_yaw"] = np.floor(
-        lidar_data["horizontal_position_yaw_f"]
-    )
+    # lidar_data["horizontal_position_yaw_f"] = (
+    #     0.5
+    #     * horizontal_resolution
+    #     * (np.arctan2(lidar_data["y"], lidar_data["x"]) / pi + 1.0)
+    # )
+    # lidar_data["horizontal_position_yaw"] = np.floor(
+    #     lidar_data["horizontal_position_yaw_f"]
+    # )
     lidar_data["vertical_position"] = np.floor(
         lidar_data["original_id"] / horizontal_resolution
     )
     # fov = 32 * pi / 180
     # fov_down = 17 * pi / 180
-    fov = 31.76 * pi / 180
-    fov_down = 17.3 * pi / 180
-    lidar_data["vertical_angle"] = np.arcsin(
-        lidar_data["z"]
-        / np.sqrt(lidar_data["x"] ** 2 + lidar_data["y"] ** 2 + lidar_data["z"] ** 2)
-    )
-    lidar_data["vertical_angle_degree"] = lidar_data["vertical_angle"] * 180 / pi
+    # fov = 31.76 * pi / 180
+    # fov_down = 17.3 * pi / 180
+    # lidar_data["vertical_angle"] = np.arcsin(
+    #     lidar_data["z"]
+    #     / np.sqrt(lidar_data["x"] ** 2 + lidar_data["y"] ** 2 + lidar_data["z"] ** 2)
+    # )
+    # lidar_data["vertical_angle_degree"] = lidar_data["vertical_angle"] * 180 / pi
 
-    lidar_data["vertical_position_pitch_f"] = vertical_resolution * (
-        1 - ((lidar_data["vertical_angle"] + fov_down) / fov)
-    )
-    lidar_data["vertical_position_pitch"] = np.floor(
-        lidar_data["vertical_position_pitch_f"]
-    )
+    # lidar_data["vertical_position_pitch_f"] = vertical_resolution * (
+    #     1 - ((lidar_data["vertical_angle"] + fov_down) / fov)
+    # )
+    # lidar_data["vertical_position_pitch"] = np.floor(
+    #     lidar_data["vertical_position_pitch_f"]
+    # )
 
-    duplicates = lidar_data[
-        lidar_data.duplicated(
-            subset=["vertical_position_pitch", "horizontal_position_yaw"],
-            keep=False,
-        )
-    ].sort_values(by=["vertical_position_pitch", "horizontal_position_yaw"])
+    # duplicates = lidar_data[
+    #     lidar_data.duplicated(
+    #         subset=["vertical_position_pitch", "horizontal_position_yaw"],
+    #         keep=False,
+    #     )
+    # ].sort_values(by=["vertical_position_pitch", "horizontal_position_yaw"])
 
     lidar_data["normalized_range"] = 1 / np.sqrt(
         lidar_data["x"] ** 2 + lidar_data["y"] ** 2 + lidar_data["z"] ** 2
     )
     projection_data = lidar_data.pivot(
-        index="vertical_position_pitch",
-        columns="horizontal_position_yaw",
+        index="vertical_position",
+        columns="horizontal_position",
         values="normalized_range",
     )
     projection_data = projection_data.reindex(
@@ -208,8 +208,8 @@ def process_frame(args) -> Tuple[int, np.ndarray, Optional[Path]]:
         i,
         projection_data.to_numpy(),
         image_path,
-        lidar_data["vertical_position_pitch_f"].min(),
-        lidar_data["vertical_position_pitch_f"].max(),
+        # lidar_data["vertical_position_pitch_f"].min(),
+        # lidar_data["vertical_position_pitch_f"].max(),
     )
 
 
@@ -304,7 +304,7 @@ def create_projection_data(
     if render_images:
         return projection_data, rendered_images
     else:
-        return (projection_data,)
+        return projection_data, None
 
 
 def main() -> int:

tools/render3d.py

@@ -1,21 +1,24 @@
-from configargparse import ArgParser, YAMLConfigFileParser, ArgumentDefaultsRawHelpFormatter
+from pathlib import Path
 from sys import exit
-from open3d.visualization.rendering import OffscreenRenderer, MaterialRecord
+
+import matplotlib.pyplot as plt
+import numpy as np
+from configargparse import (
+    ArgParser,
+    ArgumentDefaultsRawHelpFormatter,
+    YAMLConfigFileParser,
+)
 from open3d.io import read_pinhole_camera_parameters, write_image
 from open3d.utility import Vector3dVector
-from pathlib import Path
+from open3d.visualization.rendering import MaterialRecord, OffscreenRenderer
 from pointcloudset import Dataset
 from rich.progress import track
 
-import matplotlib.pyplot as plt
-
-import numpy as np
-
 from util import (
-    load_dataset,
-    existing_file,
-    create_video_from_images,
     calculate_average_frame_rate,
+    create_video_from_images,
+    existing_file,
+    load_dataset,
 )
 
 
@@ -42,14 +45,18 @@ def render_3d_images(
         distances = np.linalg.norm(points, axis=1)
         max_distance = distances.max()
         min_distance = distances.min()
-        normalized_distances = (distances - min_distance) / (max_distance - min_distance)
+        normalized_distances = (distances - min_distance) / (
+            max_distance - min_distance
+        )
         colors = plt.get_cmap("jet")(normalized_distances)[:, :3]
         pcd.colors = Vector3dVector(colors)
         return pcd
 
     rendered_images = []
 
-    for i, pc in track(enumerate(dataset, 1), description="Rendering images...", total=len(dataset)):
+    for i, pc in track(
+        enumerate(dataset, 1), description="Rendering images...", total=len(dataset)
+    ):
         o3d_pc = pc.to_instance("open3d")
         o3d_pc = color_points_by_range(o3d_pc)
         renderer.scene.add_geometry("point_cloud", o3d_pc, MaterialRecord())
@@ -68,19 +75,35 @@ def main() -> int:
         formatter_class=ArgumentDefaultsRawHelpFormatter,
         description="Render a 3d representation of a point cloud",
     )
-    parser.add_argument("--render-config-file", is_config_file=True, help="yaml config file path")
-    parser.add_argument("--input-bag-path", required=True, type=existing_file, help="path to bag file")
     parser.add_argument(
-        "--tmp-files-path", default=Path("./tmp"), type=Path, help="path temporary files will be written to"
+        "--render-config-file", is_config_file=True, help="yaml config file path"
     )
     parser.add_argument(
-        "--output-images", type=bool, default=True, help="if rendered frames should be outputted as images"
+        "--input-bag-path", required=True, type=existing_file, help="path to bag file"
     )
     parser.add_argument(
-        "--output-images-path", default=Path("./output"), type=Path, help="path rendered frames should be written to"
+        "--tmp-files-path",
+        default=Path("./tmp"),
+        type=Path,
+        help="path temporary files will be written to",
     )
     parser.add_argument(
-        "--output-video", type=bool, default=True, help="if rendered frames should be outputted as a video"
+        "--output-images",
+        type=bool,
+        default=True,
+        help="if rendered frames should be outputted as images",
+    )
+    parser.add_argument(
+        "--output-images-path",
+        default=Path("./output"),
+        type=Path,
+        help="path rendered frames should be written to",
+    )
+    parser.add_argument(
+        "--output-video",
+        type=bool,
+        default=True,
+        help="if rendered frames should be outputted as a video",
     )
     parser.add_argument(
         "--output-video-path",
@@ -88,7 +111,12 @@ def main() -> int:
         type=Path,
         help="path rendered video should be written to",
     )
-    parser.add_argument("--output-images-prefix", default="2d_render", type=str, help="filename prefix for output")
+    parser.add_argument(
+        "--output-images-prefix",
+        default="2d_render",
+        type=str,
+        help="filename prefix for output",
+    )
     parser.add_argument(
         "--camera-config-input-json-path",
         default="./saved_camera_settings.json",
@@ -110,12 +138,21 @@ def main() -> int:
     dataset = load_dataset(args.input_bag_path)
 
     images = render_3d_images(
-        dataset, args.camera_config_input_json_path, args.tmp_files_path, args.output_images_prefix
+        dataset,
+        args.camera_config_input_json_path,
+        args.tmp_files_path,
+        args.output_images_prefix,
     )
 
     if args.output_video:
-        input_images_pattern = f"{args.tmp_files_path / args.output_images_prefix}_%04d.png"
-        create_video_from_images(input_images_pattern, args.output_video_path, calculate_average_frame_rate(dataset))
+        input_images_pattern = (
+            f"{args.tmp_files_path / args.output_images_prefix}_%04d.png"
+        )
+        create_video_from_images(
+            input_images_pattern,
+            args.output_video_path,
+            calculate_average_frame_rate(dataset),
+        )
 
     if not args.output_images:
         for image in images:

tools/util.py

@@ -14,6 +14,12 @@ def load_dataset(
     return Dataset.from_file(bag_file_path, topic=pointcloud_topic)
 
 
+def save_dataset(dataset: Dataset, output_file_path: Path):
+    if not output_file_path.is_dir():
+        raise ArgumentTypeError(f"{output_file_path} has to be a valid folder!")
+    dataset.to_file(output_file_path)
+
+
 def calculate_average_frame_rate(dataset: Dataset):
     timestamps = dataset.timestamps
     time_deltas = [