1

I have two sets of discrete points in spherical coordinates, each representing top and bottom surfaces of an object.

illustration for the point cloud

I am trying to create volume from these points to separate points which lies inside and outside the object. Any suggestions where to look or which library to use?

Blue and red points represents top and bottom surfaces. Red points are generated by shifting top surface radially downwards with some constant radius.

xdze2
  • 3,986
  • 2
  • 12
  • 29
GTR
  • 13
  • 5
  • didnt convex hull work? – DrBwts Aug 27 '18 at 10:26
  • @DrBwts I created volume using convex hull but it is not properly distinguishing the points inside and outside the volume. Here [outcome] (https://imgur.com/Z90Cqff) of the convex-hull algorithm. In image fig: 3 as points which lies outside the two surfaces. – GTR Aug 27 '18 at 11:13
  • @GTR could you explain what we see on the graphs? The blue and red surfaces? The mesh? It's look like the top (blue?) and bottom (red?) surfaces are actually spherical but with different radius? Is a test on the radius not sufficient to test if the point is in-between 'blue' and 'red'? – xdze2 Aug 27 '18 at 13:51
  • it would be helpful if you put the images and explanation in the OP – DrBwts Aug 27 '18 at 14:47
  • @xdze2 Blue and red points represents top and bottom surfaces. Red points are generated by shifting top surface radially downwards with some constant radius. Mesh is just for representation or visualisation purpose. At current resolution, surfaces looks like spherical but it is flat in some places. Since surfaces are not entirely spherical, radius testing is futile. – GTR Aug 27 '18 at 23:56
  • @xdze2 Is it possible to create top and bottom surfaces and vertical walls then join them to form volume of the object to check the position the points? Scipy has inbuilt spatial function to generate surfaces in spherical coordinates. After generating the surfaces I am not sure how to proceed to create volume. – GTR Aug 28 '18 at 00:08
  • I will go this way: First, find how to mesh the blue surface (it is non-convex). Then, extrude this surface to create the volume (copy, translate, create new volume elements between old and new elements). I would look at, for instance, the CGAL library instead of scipy... I think the first step is the more complex, because the area on the spherical surface is non-convex – xdze2 Aug 28 '18 at 08:08
  • and the volume you want to create is non-convex, so convex hull is not going to work – xdze2 Aug 28 '18 at 08:09
  • @xdze2 a surface patch on a sphere with the normals pointing outwards is convex – DrBwts Aug 28 '18 at 10:22
  • 1
    @DrBwts I mean in the sense of a [convex set](https://en.wikipedia.org/wiki/Convex_set), that "for every pair of points within the region, every point on the straight line segment that joins the pair of points is also within the region". For instance, the boundary along the spherical surface, which looks like a coast line, is not convex: a convex-hull algorithm on this will short-cut all the re-entering parts? – xdze2 Aug 28 '18 at 10:34
  • @xdze2, DrBwts Thanks for the suggestions and discussion. It is difficult for me use CGAL library to solve the problem because I am using python based FEM code developed by small community through docker container. It is bit complicated to install a new library and its dependent packages. I found this [paper] (http://www.iis.sinica.edu.tw/page/jise/2012/201205_10.pdf) online with 3D concave hull algorithm. I will start with this paper and update you soon. – GTR Aug 29 '18 at 05:03
  • What do all the dots represent ? And if red points are generated from the blue points, why do you need to separate them later rather than keeping them apart ? –  Aug 29 '18 at 06:24
  • I don't see a need to stay in spherical coordinate. Seeing the angles as Cartesian coordinates wouldn't change the topology of the surface it seems. –  Aug 29 '18 at 06:28
  • @YvesDaoust dots are the points I am trying separate based on their position whether it lies inside or outside the volume created by the surfaces (blue + red points). Please have a look at the image in second comment. I will agree with your second comment. I just mentioned the coordinate system in which points and surfaces are present. – GTR Aug 29 '18 at 10:25

2 Answers2

0

If I am right, the blue and red surfaces are meshed (and watertight). So for every point you can draw the line from the sphere center and look for intersections with the mesh. This is done by finding the two triangles such that the line pierces them (this can be done by looking at the angular coordinates only, using a point-in-triangle formula), then finding the intersection points. Then it is an easy matter to classify the point as before the red surface, after the blue or in between.

Exhaustive search for the triangles can be costly. You can speed it up for instance using a hierarchy of bounding boxes or similar device.

0

Here is a custom tinkered method which may works at the condition that the average distance between points in the original surface is much smaller than the thickness of the volume and than the irregularities on the surface contour. In other words, that there are a lot of points describing the blue surfaces.

import matplotlib.pylab as plt

import numpy as np
from scipy.spatial import KDTree

# Generate a test surface:
theta = np.linspace(3, 1, 38)
phi = np.zeros_like(theta)
r = 1 + 0.1*np.sin(8*theta)

surface_points = np.stack((r, theta, phi), axis=1)  #  n x 3 array

# Generate test points:
x_span, y_span = np.linspace(-1, 0.7, 26), np.linspace(0.1, 1.2, 22)
x_grid, y_grid = np.meshgrid(x_span, y_span)

r_test = np.sqrt(x_grid**2 + y_grid**2).ravel()
theta_test = np.arctan2(y_grid, x_grid).ravel()
phi_test = np.zeros_like(theta_test)

test_points = np.stack((r_test, theta_test, phi_test), axis=1)  #  n x 3 array

# Determine if the test points are in the volume:

volume_thickness = 0.2  # Distance between the two surfaces
angle_threshold = 0.05  # Angular threshold to determine for a point 
                        # if the line from the origin to the point 
                        # go through the surface

# Get the nearest point: (replace the interpolation)
get_nearest_points = KDTree(surface_points[:, 1:]) # keep only the angles
# This is based on the cartesian distance,
# and therefore not enterily valid for the angle between points on a sphere
# It could be better to project the points on a unit shpere, and convert 
# all coordinates in cartesian frame in order to do the nearest point seach...

distance, idx = get_nearest_points.query(test_points[:, 1:])

go_through = distance < angle_threshold

nearest_surface_radius = surface_points[idx, 0]

is_in_volume = (go_through) & (nearest_surface_radius > test_points[:, 0]) \
                & (nearest_surface_radius - volume_thickness < test_points[:, 0])

not_in_volume = np.logical_not(is_in_volume)

# Graph;
plt.figure(figsize=(10, 7))
plt.polar(test_points[is_in_volume, 1], test_points[is_in_volume, 0], '.r',
          label='in volume');
plt.polar(test_points[not_in_volume, 1], test_points[not_in_volume, 0], '.k',
          label='not in volume', alpha=0.2);
plt.polar(test_points[go_through, 1], test_points[go_through, 0], '.g',
          label='go through', alpha=0.2);
plt.polar(surface_points[:, 1], surface_points[:, 0], '.b',
          label='surface');
plt.xlim([0, np.pi]); plt.grid(False);plt.legend();

The result graph, for 2D case, is:

example test

The idea is to look for each test point the nearest point in the surface, by considering only the direction and not the radius. Once this "same direction" point is found, it's possible to test both if the point is inside the volume along the radial direction (volume_thickness), and close enough to the surface using the parameter angle_threshold.

I think it would be better to mesh (non-convex) the blue surface and perform a proper interpolation, but I don't know Scipy method for this.

xdze2
  • 3,986
  • 2
  • 12
  • 29