version used in the paper

.gitignore

@@ -0,0 +1,119 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+src/*
+dependencies/*/
+dependencies/*.sh
+rendering/*.*
+pointcloud/*.*
+# C extensions
+*.so
+*.ply
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+.hypothesis/
+.pytest_cache/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# pyenv
+.python-version
+
+# celery beat schedule file
+celerybeat-schedule
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+
+# PyCharm
+.idea/
+
+# 2to3 backup files
+*.bak
+
+# eclipse
+.project
+.pydevproject
+.settings/*
+
+# h5 files used in bv/imageStack.py
+*.h5

3D-Tools/.gitignore

@@ -0,0 +1 @@
+*.csv

3D-Tools/ComparePoses.py

@@ -0,0 +1,269 @@
+import numpy as np
+import cv2
+import glob
+import matplotlib.pyplot as plt
+import open3d as o3d
+import concurrent.futures
+import sys
+
+
+def get_files(path="../pointcloud", rotation=None, num_rotations=1):
+    # Get file paths for reflection images, point clouds, and transformation matrices
+    reflection_files = sorted(glob.glob(f"{path}**/Cam*.png"))
+    depth_files = sorted(glob.glob(f"{path}**/Depth*.png"))
+    height_map_files = sorted(glob.glob(f"{path}**/HeightMap*.png"))  # HeightMap files
+    point_cloud_files = sorted(glob.glob(f"{path}**/*.ply"))
+    transformation_file = sorted(glob.glob(f"{path}/values.csv"))[0]
+    backup_file = sorted(glob.glob(f"{path}/features.csv"))
+
+    # Load transformation matrices
+    transformations = np.genfromtxt(transformation_file, delimiter=";", dtype=float)
+    transformations = [np.reshape(T, (4, 4)) for T in transformations]
+
+    # Load backup features if available
+    backup = np.genfromtxt(backup_file[0], delimiter=";", dtype=float) if backup_file else None
+
+    # Determine the maximum number of files to process
+    max_count = min(len(height_map_files), len(reflection_files))
+
+    print(max_count)
+    # Handle rotation and number of rotations
+    if rotation is None:
+        rotation = 0
+    else:
+        rotation = int(rotation % num_rotations)
+
+    # Return file paths and backup data (if available)
+    if backup is not None and len(backup) == max_count:
+        return (reflection_files[rotation:max_count:num_rotations],
+                depth_files[rotation:max_count:num_rotations],
+                height_map_files[rotation:max_count:num_rotations],
+                point_cloud_files[rotation:max_count:num_rotations],
+                transformations[rotation:max_count:num_rotations],
+                backup)
+
+    return (reflection_files[rotation:max_count:num_rotations],
+            depth_files[rotation:max_count:num_rotations],
+            height_map_files[rotation:max_count:num_rotations],
+            point_cloud_files[rotation:max_count:num_rotations],
+            transformations[rotation:max_count:num_rotations],
+            None)
+
+
+def rotation_matrix_to_euler_angles(R):
+    # Calculate Euler angles from rotation matrix
+    sy = np.sqrt(R[0, 0]**2 + R[1, 0]**2)
+    singular = sy < 1e-6
+
+    if not singular:
+        x = np.arctan2(R[2, 1], R[2, 2])
+        y = np.arctan2(-R[2, 0], sy)
+        z = np.arctan2(R[1, 0], R[0, 0])
+    else:
+        x = np.arctan2(-R[1, 2], R[1, 1])
+        y = np.arctan2(-R[2, 0], sy)
+        z = 0
+
+    return np.array([x, y, z])
+
+
+def process_reflection(reflection_path, depth_path, height_map_path, point_cloud_path, transformation, index, scale, max_sensor_area, total_computations, read_pointcloud=False):
+    # Read and normalize reflection image
+    reflection = cv2.imread(reflection_path, -1) / scale
+    depth = cv2.imread(depth_path, -1)
+    reflection = depth>0 + reflection
+    height_map = hm = cv2.imread(height_map_path, -1)  # Load height map
+
+    depth_middle = depth[depth.shape[0]//2,depth.shape[1]//2]
+    hm_middle = hm[hm.shape[0]//2,hm.shape[1]//2]
+
+    if hm_middle>0 and depth_middle>0:
+            depth=depth/depth_middle
+            height_map=hm/hm_middle
+    else:
+            depth=depth/255
+            height_map=hm/255
+
+    # Compute areas
+    multireflection_area = np.sum(reflection > 1).astype(float)
+    normal_reflection_area = np.sum(reflection == 1).astype(float)
+
+    reconstructed_simulated_area = np.sum(height_map>0)/np.sum(depth>0)
+    # Calculate clean ratio
+    clean_ratio = normal_reflection_area / (multireflection_area + normal_reflection_area)
+
+
+    reconstructed_area=np.sum(depth>0)
+    if reconstructed_area>100:
+        diff=np.abs(height_map-depth)
+        diff[height_map==0]=0
+        height_map_ratio = np.mean(diff)  # Compute some ratio from height map
+    else:
+        height_map_ratio=1
+
+    features = (
+        index,
+        clean_ratio,
+        normal_reflection_area / max_sensor_area,
+        (normal_reflection_area + multireflection_area) / max_sensor_area,
+        reconstructed_simulated_area,
+        height_map_ratio  # Include height map ratio
+    )
+
+    print(f"\x1b[2K\r {index + 1}/{total_computations}", end="")
+    return features
+
+
+def extract_features(reflection_paths, depth_paths, height_map_paths, point_cloud_paths, transformations, backup, path, threaded=True):
+    # Read the first reflection image to determine the scale
+    first_image = cv2.imread(reflection_paths[0], -1)
+    scale = np.unique(first_image)[1]
+    max_sensor_area = int(first_image.shape[0] * first_image.shape[1])
+    features = []
+
+    if backup is None:
+        print("Creating features from data")
+        if not threaded:
+            # Process reflections sequentially
+            for index, (reflection_path, depth_path, height_map_path, point_cloud_path, transformation) in enumerate(zip(reflection_paths, depth_paths, height_map_paths, point_cloud_paths, transformations)):
+                features.append(process_reflection(reflection_path, depth_path, height_map_path, point_cloud_path, transformation, index, scale, max_sensor_area, len(reflection_paths)))
+        else:
+            # Process reflections in parallel
+            with concurrent.futures.ThreadPoolExecutor() as executor:
+                futures = [
+                    executor.submit(process_reflection, reflection_path, depth_path, height_map_path, point_cloud_path, transformation, index, scale, max_sensor_area, len(reflection_paths))
+                    for index, (reflection_path, depth_path, height_map_path, point_cloud_path, transformation) in enumerate(zip(reflection_paths, depth_paths, height_map_paths, point_cloud_paths, transformations))
+                ]
+                for future in concurrent.futures.as_completed(futures):
+                    features.append(future.result())
+
+        features = np.array(features)
+        np.savetxt(f"{path}/features.csv", features, delimiter=";", fmt="%32.32f")
+    else:
+        print("Loaded features from file")
+        features = backup
+
+    return features
+
+
+def show_poses(points, features, mesh_path=None):
+    print(f"Showing Poses - min: {np.min(features)} - max: {np.max(features)} Value of feature")
+    features = np.array(features)
+    norm_features = (features - np.min(features)) / (np.max(features) - np.min(features))
+
+    # Create a dictionary to store the highest feature value for each point
+    unique_points = {}
+    for point, feature in zip(points, norm_features):
+        point_tuple = tuple(point)
+        if point_tuple not in unique_points or unique_points[point_tuple] < feature:
+            unique_points[point_tuple] = feature
+
+    # Extract the points and features with the highest feature values
+    points = np.array([point for point in unique_points.keys()])
+    features = np.array([unique_points[point] for point in unique_points.keys()])
+
+    point_cloud = o3d.geometry.PointCloud()
+    point_cloud.points = o3d.utility.Vector3dVector(points)
+    colors = plt.get_cmap("cool")(features)[:, :3]
+    point_cloud.colors = o3d.utility.Vector3dVector(colors)
+
+    # Visualize the point cloud with bigger points
+    vis = o3d.visualization.Visualizer()
+    vis.create_window()
+    vis.add_geometry(point_cloud)
+    if mesh_path is not None:
+        mesh = o3d.io.read_triangle_mesh(mesh_path)
+        vis.add_geometry(mesh)
+
+    # Change the render option to make the points bigger
+    render_option = vis.get_render_option()
+    render_option.point_size = 15.0
+
+    # Enable shading
+    render_option.mesh_shade_option = o3d.visualization.MeshShadeOption.Color
+    vis.run()
+    vis.destroy_window()
+
+
+if __name__ == "__main__":
+    # Take first argument in the cli command as path string if available
+    path = sys.argv[1] if len(sys.argv) > 1 else "../pointcloud"
+
+    reflection_paths, depth_paths, height_map_paths, point_cloud_paths, transformations, backup = get_files(path)
+
+    features = extract_features(reflection_paths, depth_paths, height_map_paths, point_cloud_paths, transformations, backup, path)
+
+    #sort features by features[-2]:= Reconstructed/Visible
+    points = [np.reshape(T, (4, 4))[:3, 3] for T in transformations]
+    combined = list(zip(features, points))
+    combined = sorted(combined, key=lambda x: x[0][-1], reverse=True)
+
+    # Extract the sorted transformations
+    points = [np.array(item[1]) for item in combined]
+    features = [item[0] for item in combined]
+
+    first_image = cv2.imread(reflection_paths[0], -1)
+    scale = np.unique(first_image)[1]
+
+    clean_reflection_ratio = [feature[1] for feature in features]
+    clean_total_ratio = [feature[2] for feature in features]
+    height_map_ratio = [feature[-1] for feature in features]
+
+    plt.plot(height_map_ratio,clean_total_ratio,".")
+    plt.xlabel("Average Probing Error")
+    plt.ylabel("Completness")
+    plt.show()
+
+    # Show Poses
+    show_poses(points[:], clean_total_ratio[:], mesh_path="../pyrefloid/res/data/cad.stl")
+
+    # Show Histograms
+    plt.subplot(121)
+    plt.xlabel("Points without Multireflection vs all Points")
+    plt.hist(clean_reflection_ratio, int(np.sqrt(len(clean_reflection_ratio))))
+    plt.subplot(122)
+    plt.xlabel("Points without Multireflection vs Sensor Area")
+    plt.hist(clean_total_ratio, int(np.sqrt(len(clean_total_ratio))))
+    plt.show()
+
+    # Show best reflections
+    for index, clean_reflected_ratio, clean_total_ratio, reflected_total_area, simulated_reflected_area, height_map_ratio in features[::]:
+
+        reflection = cv2.imread(reflection_paths[int(index)], -1) / scale
+
+        depth = cv2.imread(depth_paths[int(index)], -1) / 255
+        hm = cv2.imread(height_map_paths[int(index)], -1) / 255
+
+        depth_middle = depth[depth.shape[0]//2,depth.shape[1]//2]
+        hm_middle = hm[hm.shape[0]//2,hm.shape[1]//2]
+
+        if hm_middle>0 and depth_middle>0:
+            depth=depth/depth_middle
+            hm=hm/hm_middle
+        else:
+            depth=depth/255
+            hm=hm/255
+
+        print(hm[hm.shape[0]//2,hm.shape[1]//2],depth[depth.shape[0]//2,depth.shape[1]//2])
+
+        print(f"{int(index)}\t-\t Clean/Reflected: {clean_reflected_ratio:3.3f}\tClean/Area: {clean_total_ratio:3.5f}\tReflected/Area: {reflected_total_area:3.3f}\tReconstructed/Visible: {simulated_reflected_area}\t Average Difference Sim: {height_map_ratio:3.3f}")
+
+        plt.figure("10 best poses (descending)")
+        ax=plt.subplot(221)
+        plt.title("Depthmap")
+        plt.imshow(depth,cmap="turbo")
+        plt.colorbar()
+        plt.subplot(222,sharex=ax,sharey=ax)
+        plt.title("Simulated Reconstruction")
+        plt.imshow(hm,cmap="turbo")
+        plt.subplot(223,sharex=ax,sharey=ax)
+        plt.title("Difference Reconstruction/GT")
+        diff=depth-hm
+        diff[hm==0]=0
+        plt.imshow(diff,cmap="turbo")
+        plt.colorbar()
+        plt.subplot(224,sharex=ax,sharey=ax)
+        plt.title("Mutlireflections")
+        plt.imshow(reflection,cmap="turbo")
+        plt.colorbar()
+        plt.show()