How to identity all surface voxels from a VoxelGrid object effectively?

### Checklist

- [X] I have searched for [similar issues](https://github.com/isl-org/Open3D/issues).
- [X] For Python issues, I have tested with the [latest development wheel](http://www.open3d.org/docs/latest/getting_started.html#development-version-pip).
- [X] I have checked the [release documentation](http://www.open3d.org/docs/release/) and the [latest documentation](http://www.open3d.org/docs/latest/) (for `master` branch).

### My Question

I'm going to do a project by voxelizing a mesh object. One project step is to evaluate the voxel grid quality by surface voxels normals consistency.

Actually, I created a function to filter surface voxels. But, the problem is that as the voxel size decreases, there is a sudden increase in the time spent on the procedure. I tried to optimize the function, but it didn't make significant work.

My old solution:
1. load the .stl file and convert it to a mesh object
2. create a voxel grid object based on the mesh with a size (0.1)
3. for each voxel in the voxel grid:
3.1 locating the voxel index
3.2 finding the voxel's neighbors index by a custom function (at first it would find all 26 neighbors, but now, it only find 6 neighbors)
3.3 check if all neighbors' voxel index in the voxel grid, if any neighbor is in a "no record" state, then the voxel is a surface voxel

my [sample file](https://drive.google.com/file/d/1-x90cyibmVRAgQ_yId3emHnnCZo6J0Vo/view?usp=sharing)

Below is my code:
```
import open3d as o3d
import numpy as np

# settings
sample_path = 'data/sample.stl'
voxel_size = 0.1

# custom function to generate a voxel's neighbors' index
def get_voxel_neighbors(index):
x_neighbors_1 = index - [1, 0, 0]
x_neighbors_2 = index + [1, 0, 0]
y_neighbors_1 = index - [0, 1, 0]
y_neighbors_2 = index + [0, 1, 0]
z_neighbors_1 = index - [0, 0, 1]
z_neighbors_2 = index + [0, 0, 1]
return [x_neighbors_1, x_neighbors_2, y_neighbors_1, y_neighbors_2, z_neighbors_1, z_neighbors_2]

# get mesh
mesh = o3d.io.read_triangle_mesh(sample_path)
mesh = mesh.compute_triangle_normals()
total_triangles = len(mesh.triangles)
print('the mesh contains {} triangles'.format(total_triangles))

# create voxel grid
voxel_grid = o3d.geometry.VoxelGrid.create_from_triangle_mesh(mesh, voxel_size=voxel_size)
voxel_indices = get_voxel_indices(voxel_grid)
print(f'the voxel grid contains {len(voxel_indices)} voxels')

# get surface voxels
surface_voxels = []
for voxel in voxel_grid.get_voxels():
score = 0
index = voxel.grid_index
surface_voxel = SurfaceVoxel(voxel_size=voxel_size, voxel_index=index, mesh=mesh)
if surface_voxel.quick_check():
neighbors = get_voxel_neighbors(index)
for neighbor in neighbors:
if isin_check(voxel_indices, neighbor):
score += 1
else:
surface_voxel.add_surface_id(neighbor)
if score != 6:
surface_voxels.append(surface_voxel)
total_surface_voxel_count = len(surface_voxels)
print('surface voxel count:', total_surface_voxel_count)
```