Find all integer values contained within an irregular volume

I'm setting up a 3D volume within which points will be distributed, and want to mark some areas as "inaccessible". The areas are generated using e.g. a parametric description of ellipsoids:

% starting parameters
axes = randi(10,1,3) ; % generate axes lengths
theta = linspace(-0.5*pi,0.5*pi) ; % eccentric anomaly
lambda = linspace(0,2*pi) ; % azimuthal angle
angles = combvec(theta,lambda) ; % all angles for evaluation
% evaluate ellipsoid points
points = axes .* [cos(angles(1,:)) .* cos(angles(2,:)) ; ...
    cos(angles(1,:)) .* sin(angles(2,:)) ; ...
    sin(angles(1,:))] ;
% convert to integers, find unique rows
unique_points = unique(floor(points), 'rows') ;

At this point, I have a set of integer points for an ellipsoid centred on the origin. But I would also like all integers contained within these points. Additionally, in the real version, these ellipsoids can be asymmetric across the origin, are rotated, and offset. They might also be cylinders. The shape is determined by the user at initialisation.

Lacking a precise equation to describe such an irregular shape, how can I find all integer values contained within the ellipsoid? The only way I can think of at the moment requires nested for-loops, which takes TIME. I suspect there is a way to vectorise this code, but I can't see it at the moment.

EDIT: clarified use of "irregular shape/volume" whilst giving an example ellipsoid.

[John D'Errico]

"I'm setting up a 3D volume within which points will be distributed"

That seems to imply you want to eventually generate random numbers, or something like that. What are you really looking to do in the end with this?

Will the volume ALWAYS be an ellipsoid? Or will it be irregular? Thus anything? It seems like you are using an ellipsoid, so admit it is so, instead of just calling it irregular. Something that is truly irregular is vastly different from an ellipsoid, and you may be able to use the information about it being an ellipsoid if it truly is one.

So... Are you looking to generate the set of all points on an integer lattice that lie within a fully general ellipsoid in 3-d?

If so, then what you are doing does not seem to achieve that, thus working with floating point grids in spherical coordinates, then converting to cartesian, then rounding. If you don't want to miss points inside that volume, then you need to work in cartesian coordinates from the start.

And so, my next question is, how large is the ellipsoid? That is, unless it has a very large volume, the simplest and probably fastest solution really is to just use meshgrid/ndgrid, and then test the set of all lattice points in that volume, accepting that you will be testing many points. Since the test will be fully vectorized, it will still be fast even if the sampling is overly done. 3-d is not that huge of a space to search, and 10-1 is not that big of an aspect ratio, if that is as large as it gets.

[Paul Gratrex]

To explain the reason for using both floating-point and integers: the integers are for a grid. Using a grid means that when floating points are placed within the volume, and I check how far this point is from earlier points, only a finite area of the grid needs to be searched. This is necessary for issues of scale, as, ultimately, each grid edge will be several hundred units long. A further reason for switching between Cartesian and integers is that the inaccessible areas need not be centred on an integer.

The shape of the grid is defined by user input, and is (at the moment) anything from a sphere to a cylinder to an ovoid within which all symmetries are broken. I included an ellipsoid in the example as an illustration. Perhaps my post was poorly worded.

The ellipsoids (or general shapes) use a rng for their dimensions, typically using e.g.

radius = random('Normal',25,5) ;

The entire volume generated, which as I said is arbitrary NxN where N>>10, has a grid to simplify searches. There will be thousands of points inside, with constraints on the generation of new points. One such constraint is that it should not be within the e.g. ellipsoid areas generated at the initialisation of the volume. Thus, I want a fast check of whether or not a given point is inside one of the e.g. ellipsoid areas. Additionally, I want a fast check of how much of the total NxN volume is contained within the e.g. ellipsoid areas.

Answer 1

[Cris LaPierre]

I wonder if the functions inpolygon and inpolyhedron might be helpful here. See this post about finding the points inside a 3D volume as well.

[Paul Gratrex]

Perhaps you can help me understand how to use inpolyhedron?

It seems I need to feed my Nx3 matrix into an isosurface, but I'm not sure how to do that. isosurface asks for a 3D array. Converting an Nx3 matrix into a 3D array seems to require a second array converting e.g. the z-indices into z-values, since the original matrix can contain negative integers. It also might require a prohibitively large matrix, multiple GB.

If you could perhaps provide a MWE of how to use inpolyhedron on an ellipsoid, as I described it above?

Thanks.

EDIT: I've got it working for inpolygon, which is a significant speed boost compared with my original nested for-loop. So, thanks for that suggestion. I'd still be interested in using inpolyhedron if you can demonstrate its usage in this situation. Thanks!

[Cris LaPierre]

MATLAB has functions for generating ellipsoids and cylinders. Perhaps those will be useful as well.

I found the following pages describing inpolyhedron. The challenge would be getting your shapes to be triangulated meshes.

https://blogs.mathworks.com/pick/2013/09/06/inpolyhedron/

https://blogs.mathworks.com/community/2013/06/20/paul-prints-the-l-shaped-membrane/

[Paul Gratrex]

Thanks for the help, I think I'll probably end up using MATLAB's built-in functions at a later date. In the meantime, inpolygon works, just about.