Voxel binary dilation solves sampling space binary voxel dilation in 3D space

binary voxel dilation in 3D space

  • foreword
  • 1 Binary inflation
  • 2 voxel binary dilation
    • 2.1 Difficulties
    • 2.2 Write a voxel expansion function based on ndimage.binary_dilation()
  • 3 Voxel expansion test code and effect

Foreword

In the process of doing the project, I encountered a problem: it is known that there is a 3D point cloud model of the target object. It is necessary to expand the model, expand its space, and then sample the three-dimensional coordinate points in the newly expanded space. The illustration is as follows, the Target needs to be expanded to the purple and green boundaries respectively, and the sampling space between green and purple is the sampling space I want.
Schematic diagram of space expansion


Schematic diagram of space expansion

To solve this problem, the method of voxel binary expansion is used. The general process is: first convert the point cloud into voxels, extract the index of the voxels, and use the expansion operation in morphology to solve the voxels after expansion. Index, and then calculate the coordinates of each voxel to obtain the expansion space. However, since the actual voxel is different from the binary expansion in the image, the solution process is cumbersome, and some conversions are required in some places to perform the expansion operation.

1 binary expansion

First understand what binary expansion is, and the library function scipy.ndimage.binary_dilation in python.
Dilation is the process of merging all background pixels (black areas, 0 in binary) that touch an object (white area, 1 in binary) into that object.
Image expansion diagram


Image Expansion Schematic

Function of expansion operation: scipy.ndimage.binary_dilation(input, structure=None, iterations=1, mask=None, output=None, border_value=0, origin=0, brute_force=False), where input is the binary value to be expanded Matrix, structure is an expansion structure, and other parameters can be consulted by yourself.
The following example illustrates the usage of the function.
First import the library, create a 5*5 matrix a at will, and set a[2, 2] = 1, and the remaining elements are 0, and expand.

from scipy import ndimage
import numpy as np

a = np. zeros((5, 5))
a[2, 2] = 1
b = ndimage. binary_dilation(a)
c = ndimage.binary_dilation(a) + 0 # + 0 makes the output bool matrix into 1 and 0

a, b, and c are as follows:

After that, specify the structure to generate the expansion, and then see the effect. In the code below, the first parameter in ndimage.generate_binary_structure(2, 1) specifies a matrix of several dimensions, and the second parameter is the shape to be expanded. If you want to inflate a three-dimensional matrix, then ndimage.generate_binary_structure(3, 1).

struct1 = ndimage.generate_binary_structure(2, 1)
d = ndimage.binary_dilation(a, structure=struct1) + 0

struct2 = ndimage.generate_binary_structure(2, 2)
e = ndimage.binary_dilation(a, structure=struct2) + 0

2 voxel binary dilation

2.1 Difficulties

To perform binary dilation on voxels, the process is a little more complicated. The main reasons are as follows: (1) In the binary expansion mentioned above, the expanded elements of the matrix a to be expanded do not exceed the index boundary, but the indices in the voxels are closely arranged, so a larger three-dimensional matrix needs to be created first, and the The voxel index corresponds to it. (2) It is also necessary to solve the voxel coordinates, because binary values are used when expanding, which results in the need to convert the coordinates after calculating the index. It is necessary to understand the three-dimensional matrix, two-dimensional matrix, matrix index, and matrix values. difference, and have a certain spatial imagination.

2.2 Write voxel expansion function based on ndimage.binary_dilation()

Assuming that open3d has been used to convert point cloud data into voxel data voxel_grid, it is necessary to solve voxel coordinates voxels_xyz after expansion.
Function inputs include:
voxel_grid, voxel data.
expand_dis, the space distance to be expanded outwards, in meters, the number of iterations during expansion in the function will be determined by both expand_dis and the voxel size.
scale_factor, the reciprocal of the coordinate scaling factor when the point cloud is converted to a voxel, that is, if the scaling factor is 1/7, the coordinate values of all voxels will be reduced by 7 times to the center of the point cloud, and scale_factor=7 is sufficient when calling the function.
limit_z, add a limit condition, when the meaning of limit_z in this function, the voxel obtained after expansion will be deleted if the z coordinate is less than limit_z.
Function output:
voxels_xyz, voxel coordinates after dilation.
The complete code of the voxel dilation function is as follows:

import numpy as np
from scipy import ndimage
import math
import matplotlib.pyplot as plt

def generate_voxel_dilate_xyz(voxel_grid, expand_dis, limit_z):
    # get voxel index
    voxels = voxel_grid. get_voxels()
    voxel_coords = np.asarray([voxel.grid_index for voxel in voxels])
    # print("voxel_coords--------------------:\
", voxel_coords)
      
    voxel_size = voxel_grid.voxel_size
    dilate_iter = math.ceil(expan_dis/voxel_size) #Expansion iterations
    expand_coord = dilate_iter * 2 #The index of the expanded space array should expand to both sides of the coordinate axis
    
    # The array of pre-expansion space, first expand the voxel boundary, assign 1 to the place with voxel, and 0 to the place without voxel, so as to facilitate subsequent expansion
    # (np.max(voxel_coords, axis=0) returns the maximum value of each column of voxel_coords, that is, the maximum index value in the three directions of xyz
    expan_space = np.zeros((np.max(voxel_coords, axis=0) + 1 + expan_coord)) #np.max calculates the maximum value of the voxel index, and it should be + 1 when specifying an array
    for i in range(len(voxel_coords)):
        coord_x = voxel_coords[i, 0] + expan_coord//2
        coord_y = voxel_coords[i, 1] + expan_coord//2
        coord_z = voxel_coords[i, 2] + expan_coord//2
        expand_space[coord_x, coord_y, coord_z] = 1
    # print("expan_space-----------------------------\
",expan_space)
        
    # Generate spherical structure elements before expansion, this is the structure after expansion, and then expand
    se = ndimage.generate_binary_structure(3, 3)
    for i in range(dilate_iter):
        expand_space = ndimage.binary_dilation(expan_space, structure=se)
    # print("expan_space-----------------------------\
", expan_space + 0)
      
    # Set the benchmark index and coordinate value, then traverse, calculate the coordinates of each voxel according to the index difference and save it to the two-dimensional array voxels_xyz
    voxel0 = voxels[0].grid_index
    voxel0_xyz = voxel_grid.get_voxel_center_coordinate(voxel0)
    base_coord_x, base_coord_y, base_coord_z = voxel0 + expan_coord//2 #(0, 0, 0) in the original voxel, the index in the current voxel is 0 + expan_coord//2
    base_x, base_y, base_z = voxel0_xyz[0], voxel0_xyz[1], voxel0_xyz[2]

    voxels_xyz = np.empty((np.count_nonzero(expan_space), 3))
    point_order = 0
    for index, value in np.ndenumerate(expan_space):
        if value != 0:
            coord_x, coord_y, coord_z = index
            voxels_xyz[point_order, 0] = base_x + (coord_x - base_coord_x) * voxel_size
            voxels_xyz[point_order, 1] = base_y + (coord_y - base_coord_y) * voxel_size
            voxels_xyz[point_order, 2] = base_z + (coord_z - base_coord_z) * voxel_size
            point_order = point_order + 1
    
    
    # Delete coordinates where z<limit_z. iterate over each row
    voxels_xyz_original = np.copy(voxels_xyz) # backup original array
    condition = voxels_xyz[:, 2] >= limit_z # Use the NumPy condition index to directly filter out the rows that meet the condition
    voxels_xyz = voxels_xyz[condition]
    
    return voxels_xyz


# Call this function to draw the voxel coordinate index before and after expansion at the same time
def plot_3d_points(points1, points2, colors=None):
    # Create a 3D coordinate axis
    fig = plt. figure()
    ax = fig.add_subplot(111, projection='3d')

    # draw each point
    ax.scatter(points1[:, 0], points1[:, 1], points1[:, 2], c='red', label='Points 1', alpha=0.1)
    ax.scatter(points2[:, 0], points2[:, 1], points2[:, 2], c='green', label='Points 2', alpha=0.1)

    # set axis labels
    ax. set_xlabel('X')
    ax. set_ylabel('Y')
    ax. set_zlabel('Z')

    # add legend
    ax. legend()

    # display graphics
    plt. show()

3 Voxel expansion test code and effect

Based on open3d, the point cloud is filtered and converted into voxels, and the voxel coordinates after solving the expansion are saved as a .mat file, and then the voxel coordinates before and after expansion are used for three-dimensional drawing.
The code written in the previous section is saved as voxel_dilate.py, which is called in the next code, code:

import open3d as o3d
import numpy as np
import voxel_dilate
import copy
import scipy.io

pcd = o3d.io.read_point_cloud("E:\learn_test\MyMatlab\models\plane\plane_expans.ply")
print("start pcd:", pcd)
print(pcd. get_center(), pcd. get_min_bound(), pcd. get_max_bound())

# Point cloud to voxel
voxel_pcd = copy.deepcopy(pcd)
print(voxel_pcd. get_center(), voxel_pcd. get_min_bound(), voxel_pcd. get_max_bound())
voxel_pcd.colors = o3d.utility.Vector3dVector(np.random.uniform(0, 1, size=(len(voxel_pcd.points), 3)))
voxel_grid = o3d.geometry.VoxelGrid.create_from_point_cloud(input = voxel_pcd, voxel_size=0.4)
print("voxel_grid:", voxel_grid)
print(voxel_grid. get_center(), voxel_grid. get_min_bound(), voxel_grid. get_max_bound())
o3d.visualization.draw_geometries([voxel_grid], window_name="PointCloud Voxelization", width=800, height=600, mesh_show_back_face=True)


# Voxel expansion, call the self-written function, return the expanded coordinate value, save it as a MAT file
voxels_xyz = voxel_dilate.generate_voxel_dilate_xyz(voxel_grid, expan_dis=0.5,
                                                    limit_z=-2.5)
# scipy.io.savemat('.\sample_space\voxels_xyz_min.mat', {'array_2d': voxels_xyz})
# voxels_xyz = voxel_dilate.generate_voxel_dilate_xyz(voxel_grid, expan_dis=5.7, #(5.5^2 * 3)^(1/2) = 10
# limit_z=0.38)
# scipy.io.savemat('.\sample_space\voxels_xyz_max.mat', {'array_2d': voxels_xyz})


# Get the voxel coordinates, and then draw the old and new coordinates together
voxels = voxel_grid. get_voxels()
voxel_coords = np.asarray([voxel.grid_index for voxel in voxels])
voxel_xyz_ori = np.empty((len(voxel_coords), 3))
i = 0
for voxels in voxels:
    voxel_xyz_ori[i] = voxel_grid.get_voxel_center_coordinate(voxel.grid_index)
    i = i + 1

voxel_dilate.plot_3d_points(voxel_xyz_ori, voxels_xyz)

The renderings are as follows: