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:

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.