I have 3D points from noisy detection of one planar surface. I want to fit an oriented bounding box for these points.

I assume (it should hold in my case) that the edges are much lighter than the inner parts (i.e. higher intensity), so I approximately filter it by intensity to get only points close to the edges. But the detection is pretty noisy, there are also some surfaces and lidar noise around the edges so that makes the oriented bounding box a bit shifted and a bit rotated.

I want to apply optimization procedure to shift and rotate this original estimate so that it fits as many boundary points as possible. However, I am experience problem constructing such an optimization. I use scipy.optimize.minimize and the bounding box seems to just randomly just all over the place with the final result not matching the structure of interest if I use Powell method, and just not moving bounding box at all and staying at the original estimate if I use non-Powell method.

As I cannot by rules just dump my file here, I tried to generate toy example to show what I mean (I try to simulate both random noise within the structure + some neighboring structure; I also simulate my original non-perfect estimation of the bounding box by rotating and shifting ground truth bounding box):

import open3d as o3d

def generate_noisy_boundary_pcd(center, width, height, depth, margin, noise_level, num_points):
    """
    Generate a 3D Open3D point cloud with noisy points only on the boundary of a 3D rectangle.

    Args:
        center (tuple): Center of the rectangle (x, y, z).
        width (float): Width of the rectangle in the X direction.
        height (float): Height of the rectangle in the Y direction.
        depth (float): Thickness of the rectangle in the Z direction.
        margin (float): Thickness of the boundary region to retain.
        noise_level (float): Standard deviation of Gaussian noise.
        num_points (int): Number of points to generate.

    Returns:
        open3d.geometry.PointCloud: Open3D point cloud containing only noisy boundary points.
    """
    # Generate rectangle points in 3D space
    x_vals = np.linspace(-width / 2, width / 2, int(np.sqrt(num_points)))
    y_vals = np.linspace(-height / 2, height / 2, int(np.sqrt(num_points)))
    x_grid, y_grid = np.meshgrid(x_vals, y_vals)

    # Front and back face at ±depth/2
    z_front = np.full_like(x_grid, depth / 2)
    z_back = np.full_like(x_grid, -depth / 2)

    # Create front and back faces
    rect_points = np.column_stack((x_grid.ravel(), y_grid.ravel(), z_front.ravel()))

    # Filter only boundary points (margin region)
    boundary_mask = (
        (np.abs(x_grid) >= (width / 2 - margin)) |  # Vertical boundary
        (np.abs(y_grid) >= (height / 2 - margin))    # Horizontal boundary
    )
    
    boundary_points = rect_points[boundary_mask.ravel()]
    
    # Apply noise
    noise = np.random.normal(0, noise_level, boundary_points.shape)
    noisy_boundary_points = boundary_points + noise

    # Shift to the desired center
    noisy_boundary_points += np.array(center)

    # Convert to Open3D point cloud
    pcd = o3d.geometry.PointCloud()
    pcd.points = o3d.utility.Vector3dVector(noisy_boundary_points)

    return pcd

# Generate noisy boundary point cloud
pcd_noisy_boundary = generate_noisy_boundary_pcd(center=(0, 0, 0), width=0.8, height=1.0, depth=0.1, margin=0.1, noise_level=0.02, num_points=500)

# Visualize with Open3D
o3d.visualization.draw_geometries([pcd_noisy_boundary])


def apply_random_shift_and_rotation(obb, max_shift=0.05, max_rotation=np.radians(5)):
    """
    Apply a small random shift and rotation to an oriented bounding box.

    Args:
        obb (o3d.geometry.OrientedBoundingBox): The original oriented bounding box.
        max_shift (float): Maximum translation shift in each direction.
        max_rotation (float): Maximum rotation in radians for each axis.

    Returns:
        o3d.geometry.OrientedBoundingBox: The transformed OBB.
    """
    # Generate a small random shift
    shift = np.random.uniform(-max_shift, max_shift, 3)

    # Generate small random rotations around each axis
    delta_angle_x = np.random.uniform(-max_rotation, max_rotation)
    delta_angle_y = np.random.uniform(-max_rotation, max_rotation)
    delta_angle_z = np.random.uniform(-max_rotation, max_rotation)

    # Compute the new rotation matrix
    delta_Rx = o3d.geometry.get_rotation_matrix_from_xyz((delta_angle_x, 0, 0))
    delta_Ry = o3d.geometry.get_rotation_matrix_from_xyz((0, delta_angle_y, 0))
    delta_Rz = o3d.geometry.get_rotation_matrix_from_xyz((0, 0, delta_angle_z))

    # Apply rotation perturbation on top of existing rotation
    R_new = obb.R @ (delta_Rz @ delta_Ry @ delta_Rx)

    # Apply translation shift
    center_new = obb.center + shift

    # Create the transformed OBB
    transformed_obb = o3d.geometry.OrientedBoundingBox(center_new, R_new, obb.extent)

    return transformed_obb

# Compute the initial oriented bounding box from the noisy boundary point cloud
obb_original = o3d.geometry.OrientedBoundingBox.create_from_points(pcd_noisy_boundary.points)
obb_original.color = (1, 0, 0)  # Red for original

# Apply random shift and rotation
obb_transformed = apply_random_shift_and_rotation(obb_original)
obb_transformed.color = (0, 1, 0)  # Green for transformed

# Visualize both original and transformed bounding boxes
o3d.visualization.draw_geometries([pcd_noisy_with_cluster, obb_original, obb_transformed])

All above is just toy data generation. Below is code I am trying for optimization:

from scipy.optimize import minimize

def fit_obb_to_points(params, pcd, obb, margin=0.1):
    """
    Optimize the OBB center and refine its rotation so that it encloses the most points.

    Args:
        params (list): [cx, cy, cz, delta_angle_x, delta_angle_y, delta_angle_z].
        pcd (o3d.geometry.PointCloud): Input 3D point cloud.
        obb (o3d.geometry.OrientedBoundingBox): Initial OBB.
        expand_factor (float): Factor to expand the bounding box slightly.

    Returns:
        int: Negative count of points inside the OBB (to be minimized).
    """
    # Extract parameters
    cx, cy, cz, delta_angle_x, delta_angle_y, delta_angle_z = params
    center = np.array([cx, cy, cz])

    # Compute refined rotation by applying small delta angles to the existing rotation
    delta_Rx = o3d.geometry.get_rotation_matrix_from_xyz((delta_angle_x, 0, 0))
    delta_Ry = o3d.geometry.get_rotation_matrix_from_xyz((0, delta_angle_y, 0))
    delta_Rz = o3d.geometry.get_rotation_matrix_from_xyz((0, 0, delta_angle_z))

    # Apply small rotation adjustments on top of the original OBB rotation
    R_new = obb.R @ (delta_Rz @ delta_Ry @ delta_Rx)

    expanded_obb = o3d.geometry.OrientedBoundingBox(center, R_new, [obb.extent[0] + margin, obb.extent[1] + margin, obb.extent[2]])

    o3d.visualization.draw_geometries([expanded_obb, pcd])

    # Count inliers inside both bounding boxes
    inliers_large = expanded_obb.get_point_indices_within_bounding_box(pcd.points)

    # Compute number of boundary inliers
    boundary_inliers = len(inliers_large)
    print(boundary_inliers)
    return -boundary_inliers 
    

def optimize_bounding_box(pcd, initial_obb):
    """
    Adjusts the center and fine-tunes the rotation of the OBB so that it fits as many points as possible.

    Args:
        pcd (o3d.geometry.PointCloud): The point cloud.
        initial_obb (o3d.geometry.OrientedBoundingBox): Initial bounding box (rotation should be preserved).

    Returns:
        o3d.geometry.OrientedBoundingBox: Optimized OBB.
    """
    # Initial parameters: [center_x, center_y, center_z, delta_angle_x, delta_angle_y, delta_angle_z]
    initial_params = [
        *initial_obb.center,
        0, 0, 0  # Start with no change to rotation
    ]

    # Optimize the bounding box parameters
    result = minimize(
        fit_obb_to_points, initial_params,
        args=(pcd, initial_obb),
        # method="Powell"
    )

    # Extract optimized parameters
    optimized_center = np.array(result.x[:3])
    delta_angles = result.x[3:]

    # Compute the final refined rotation
    delta_Rx = o3d.geometry.get_rotation_matrix_from_xyz((delta_angles[0], 0, 0))
    delta_Ry = o3d.geometry.get_rotation_matrix_from_xyz((0, delta_angles[1], 0))
    delta_Rz = o3d.geometry.get_rotation_matrix_from_xyz((0, 0, delta_angles[2]))
    optimized_rotation = initial_obb.R @ (delta_Rz @ delta_Ry @ delta_Rx)

    # Create the final optimized OBB
    optimized_obb = o3d.geometry.OrientedBoundingBox(optimized_center, optimized_rotation, initial_obb.extent)

    return optimized_obb

optimized_obb = optimize_bounding_box(pcd_noisy_with_cluster, obb_transformed)
obb_transformed.color = (255, 0, 0)
optimized_obb.color = (0, 255, 0)
o3d.visualization.draw_geometries([pcd_noisy_with_cluster, obb_transformed, optimized_obb])

If you leave o3d.visualization.draw_geometries([expanded_obb, pcd]) as it is in the optimization code, you will see that either bounding box is jumping a lot, sometimes even completely away from any points (Powell method), or is just stuck for a non-Powell method.

It seems clear that this method is not really suitable. But I cannot think of some alternative. In 2D image I would apply some contour finding methods but I cannot seem to find something for 3D points. There are of course plane detection methods in open3d but these are exactly where I get my original non-perfect estimate from. I however want to find some method to refine the estimation based on what was found by those open3d methods.

1. Is it possible to fix the approach I am using?
2. Or if it is completely unsuitable, which approach should I pursue?

I have 3D points from noisy detection of one planar surface. I want to fit an oriented bounding box for these points.

I assume (it should hold in my case) that the edges are much lighter than the inner parts (i.e. higher intensity), so I approximately filter it by intensity to get only points close to the edges. But the detection is pretty noisy, there are also some surfaces and lidar noise around the edges so that makes the oriented bounding box a bit shifted and a bit rotated.

I want to apply optimization procedure to shift and rotate this original estimate so that it fits as many boundary points as possible. However, I am experience problem constructing such an optimization. I use scipy.optimize.minimize and the bounding box seems to just randomly just all over the place with the final result not matching the structure of interest if I use Powell method, and just not moving bounding box at all and staying at the original estimate if I use non-Powell method.

As I cannot by rules just dump my file here, I tried to generate toy example to show what I mean (I try to simulate both random noise within the structure + some neighboring structure; I also simulate my original non-perfect estimation of the bounding box by rotating and shifting ground truth bounding box):

import open3d as o3d

def generate_noisy_boundary_pcd(center, width, height, depth, margin, noise_level, num_points):
    """
    Generate a 3D Open3D point cloud with noisy points only on the boundary of a 3D rectangle.

    Args:
        center (tuple): Center of the rectangle (x, y, z).
        width (float): Width of the rectangle in the X direction.
        height (float): Height of the rectangle in the Y direction.
        depth (float): Thickness of the rectangle in the Z direction.
        margin (float): Thickness of the boundary region to retain.
        noise_level (float): Standard deviation of Gaussian noise.
        num_points (int): Number of points to generate.

    Returns:
        open3d.geometry.PointCloud: Open3D point cloud containing only noisy boundary points.
    """
    # Generate rectangle points in 3D space
    x_vals = np.linspace(-width / 2, width / 2, int(np.sqrt(num_points)))
    y_vals = np.linspace(-height / 2, height / 2, int(np.sqrt(num_points)))
    x_grid, y_grid = np.meshgrid(x_vals, y_vals)

    # Front and back face at ±depth/2
    z_front = np.full_like(x_grid, depth / 2)
    z_back = np.full_like(x_grid, -depth / 2)

    # Create front and back faces
    rect_points = np.column_stack((x_grid.ravel(), y_grid.ravel(), z_front.ravel()))

    # Filter only boundary points (margin region)
    boundary_mask = (
        (np.abs(x_grid) >= (width / 2 - margin)) |  # Vertical boundary
        (np.abs(y_grid) >= (height / 2 - margin))    # Horizontal boundary
    )
    
    boundary_points = rect_points[boundary_mask.ravel()]
    
    # Apply noise
    noise = np.random.normal(0, noise_level, boundary_points.shape)
    noisy_boundary_points = boundary_points + noise

    # Shift to the desired center
    noisy_boundary_points += np.array(center)

    # Convert to Open3D point cloud
    pcd = o3d.geometry.PointCloud()
    pcd.points = o3d.utility.Vector3dVector(noisy_boundary_points)

    return pcd

# Generate noisy boundary point cloud
pcd_noisy_boundary = generate_noisy_boundary_pcd(center=(0, 0, 0), width=0.8, height=1.0, depth=0.1, margin=0.1, noise_level=0.02, num_points=500)

# Visualize with Open3D
o3d.visualization.draw_geometries([pcd_noisy_boundary])


def apply_random_shift_and_rotation(obb, max_shift=0.05, max_rotation=np.radians(5)):
    """
    Apply a small random shift and rotation to an oriented bounding box.

    Args:
        obb (o3d.geometry.OrientedBoundingBox): The original oriented bounding box.
        max_shift (float): Maximum translation shift in each direction.
        max_rotation (float): Maximum rotation in radians for each axis.

    Returns:
        o3d.geometry.OrientedBoundingBox: The transformed OBB.
    """
    # Generate a small random shift
    shift = np.random.uniform(-max_shift, max_shift, 3)

    # Generate small random rotations around each axis
    delta_angle_x = np.random.uniform(-max_rotation, max_rotation)
    delta_angle_y = np.random.uniform(-max_rotation, max_rotation)
    delta_angle_z = np.random.uniform(-max_rotation, max_rotation)

    # Compute the new rotation matrix
    delta_Rx = o3d.geometry.get_rotation_matrix_from_xyz((delta_angle_x, 0, 0))
    delta_Ry = o3d.geometry.get_rotation_matrix_from_xyz((0, delta_angle_y, 0))
    delta_Rz = o3d.geometry.get_rotation_matrix_from_xyz((0, 0, delta_angle_z))

    # Apply rotation perturbation on top of existing rotation
    R_new = obb.R @ (delta_Rz @ delta_Ry @ delta_Rx)

    # Apply translation shift
    center_new = obb.center + shift

    # Create the transformed OBB
    transformed_obb = o3d.geometry.OrientedBoundingBox(center_new, R_new, obb.extent)

    return transformed_obb

# Compute the initial oriented bounding box from the noisy boundary point cloud
obb_original = o3d.geometry.OrientedBoundingBox.create_from_points(pcd_noisy_boundary.points)
obb_original.color = (1, 0, 0)  # Red for original

# Apply random shift and rotation
obb_transformed = apply_random_shift_and_rotation(obb_original)
obb_transformed.color = (0, 1, 0)  # Green for transformed

# Visualize both original and transformed bounding boxes
o3d.visualization.draw_geometries([pcd_noisy_with_cluster, obb_original, obb_transformed])

All above is just toy data generation. Below is code I am trying for optimization:

from scipy.optimize import minimize

def fit_obb_to_points(params, pcd, obb, margin=0.1):
    """
    Optimize the OBB center and refine its rotation so that it encloses the most points.

    Args:
        params (list): [cx, cy, cz, delta_angle_x, delta_angle_y, delta_angle_z].
        pcd (o3d.geometry.PointCloud): Input 3D point cloud.
        obb (o3d.geometry.OrientedBoundingBox): Initial OBB.
        expand_factor (float): Factor to expand the bounding box slightly.

    Returns:
        int: Negative count of points inside the OBB (to be minimized).
    """
    # Extract parameters
    cx, cy, cz, delta_angle_x, delta_angle_y, delta_angle_z = params
    center = np.array([cx, cy, cz])

    # Compute refined rotation by applying small delta angles to the existing rotation
    delta_Rx = o3d.geometry.get_rotation_matrix_from_xyz((delta_angle_x, 0, 0))
    delta_Ry = o3d.geometry.get_rotation_matrix_from_xyz((0, delta_angle_y, 0))
    delta_Rz = o3d.geometry.get_rotation_matrix_from_xyz((0, 0, delta_angle_z))

    # Apply small rotation adjustments on top of the original OBB rotation
    R_new = obb.R @ (delta_Rz @ delta_Ry @ delta_Rx)

    expanded_obb = o3d.geometry.OrientedBoundingBox(center, R_new, [obb.extent[0] + margin, obb.extent[1] + margin, obb.extent[2]])

    o3d.visualization.draw_geometries([expanded_obb, pcd])

    # Count inliers inside both bounding boxes
    inliers_large = expanded_obb.get_point_indices_within_bounding_box(pcd.points)

    # Compute number of boundary inliers
    boundary_inliers = len(inliers_large)
    print(boundary_inliers)
    return -boundary_inliers 
    

def optimize_bounding_box(pcd, initial_obb):
    """
    Adjusts the center and fine-tunes the rotation of the OBB so that it fits as many points as possible.

    Args:
        pcd (o3d.geometry.PointCloud): The point cloud.
        initial_obb (o3d.geometry.OrientedBoundingBox): Initial bounding box (rotation should be preserved).

    Returns:
        o3d.geometry.OrientedBoundingBox: Optimized OBB.
    """
    # Initial parameters: [center_x, center_y, center_z, delta_angle_x, delta_angle_y, delta_angle_z]
    initial_params = [
        *initial_obb.center,
        0, 0, 0  # Start with no change to rotation
    ]

    # Optimize the bounding box parameters
    result = minimize(
        fit_obb_to_points, initial_params,
        args=(pcd, initial_obb),
        # method="Powell"
    )

    # Extract optimized parameters
    optimized_center = np.array(result.x[:3])
    delta_angles = result.x[3:]

    # Compute the final refined rotation
    delta_Rx = o3d.geometry.get_rotation_matrix_from_xyz((delta_angles[0], 0, 0))
    delta_Ry = o3d.geometry.get_rotation_matrix_from_xyz((0, delta_angles[1], 0))
    delta_Rz = o3d.geometry.get_rotation_matrix_from_xyz((0, 0, delta_angles[2]))
    optimized_rotation = initial_obb.R @ (delta_Rz @ delta_Ry @ delta_Rx)

    # Create the final optimized OBB
    optimized_obb = o3d.geometry.OrientedBoundingBox(optimized_center, optimized_rotation, initial_obb.extent)

    return optimized_obb

optimized_obb = optimize_bounding_box(pcd_noisy_with_cluster, obb_transformed)
obb_transformed.color = (255, 0, 0)
optimized_obb.color = (0, 255, 0)
o3d.visualization.draw_geometries([pcd_noisy_with_cluster, obb_transformed, optimized_obb])

If you leave o3d.visualization.draw_geometries([expanded_obb, pcd]) as it is in the optimization code, you will see that either bounding box is jumping a lot, sometimes even completely away from any points (Powell method), or is just stuck for a non-Powell method.

It seems clear that this method is not really suitable. But I cannot think of some alternative. In 2D image I would apply some contour finding methods but I cannot seem to find something for 3D points. There are of course plane detection methods in open3d but these are exactly where I get my original non-perfect estimate from. I however want to find some method to refine the estimation based on what was found by those open3d methods.

1. Is it possible to fix the approach I am using?
2. Or if it is completely unsuitable, which approach should I pursue?

• python
• open3d

1 Answer 1

Your approach is hitting a pitfall did some research and...

you can do the following

Redefine the Cost Function a. Use a Continuous (Soft) Loss Function

Instead of counting inliers (a hard 0–1 indicator), define a loss that changes continuously with the box parameters. For example, for each point you might compute its “distance” to being inside the box. One strategy is to:

Compute Signed Distances: For each of the 6 planes of the OBB, calculate the signed distance of a point to that plane (positive if outside, negative if inside). Penalize Violations: Define a penalty for a point based on how far it is outside the box. For instance, you could use a function such as

 lossi=∑j=16max⁡{0,dij}2
 lossi​=j=1∑6​max{0,dij​}2

where dijdij​ is the distance of point ii from plane jj (adjusted so that points inside give negative or zero values). Aggregate Loss: Sum (or average) these penalties over all points. This loss will be zero when all points are inside and will increase smoothly as points lie outside.

2. Rotation Parameterization: Optimizing over rotations can be tricky. Rather than using three Euler angles (delta_angle_x, delta_angle_y, delta_angle_z) which can suffer from singularities and non-smooth behavior, consider:

Fix the cost function: Replace the discrete inlier count with a continuous loss (e.g. using distances from the box boundaries). Clean up the optimization loop: Remove any visualization or other side effects from the cost function. Consider a better rotation parameterization: Optimize over quaternions or rotation vectors to avoid pitfalls of Euler angles. Explore alternative strategies: Depending on your data’s characteristics (planar structure), project to 2D or use sampling/RANSAC methods to refine the bounding box.

a bunch o random links https://www.open3d./docs/latest/python_api/open3d.geometry.OrientedBoundingBox.html/

https://github/gabyx/ApproxMVBB

i cant control this auto crct it just indents me paragraphs appologies for unreadable content