Is it possible to speed up finding points inside a 3D shape (ie inShape)?

Hi-

I am hoping to figure out a way (and it may not be possible) to speed up finding points inside a 3D shape, ie by inShape as below. This is for a research lab analysis of some microscopy images. A snippet of the code from the larger analysis block is below, along with an image of the graphs that are output from this snippet. The general idea is

--Start with I1, which is an m x n x p array of the points shown in the leftmost graph, pixellated (so values = number of points in that "pixel").

--Goal is to make mask1 that is an evenly spaced m x n x p array of points with values = 1 that are in the shp, and values = 0 that are outside the shp

--The slow step is finding which points inside the mask are within the alphaShape (ie, inShape). It doesn't take so long on a single run (1-3 seconds depending on size of the original image), but this bit of code gets run frequently in the overall analysis workflow. I haven't really found any solutions to speed it up beyond cropping the set of test points down to the limits of the image like I have already. inShape doesn't run on the GPU and it's possible the overhead of transfer wouldn't help anything anyway. I recognize that this just might be the limit of working with relatively large data.

--I can upload a sample I1 if needed, but I'm wondering first if there is an obvious conceptual coding change I can make.

Thanks-

So :

% Example code
% I1 is set previously, it's an m x n x p array of 3d space where the values are pixel intensity at each point
% x1 y1 and z1 are m x 1 arrays that contain all the xy or z points of the original image acquisition after pixellation
% 
% Copy the image array and set points that are outside the max and min edges of the real image to nan. the +/-1 is a buffer to include all points
mask1 = I1; 
mask1(1:min(x1)-1,:,:) = nan;                  
mask1(:,1:min(y1)-1,:) = nan;
mask1(:,:,1:min(z1)-1) = nan;
mask1(max(x1)+1:end,:,:) = nan;
mask1(:,max(y1)+1:end,:) = nan;
mask1(:,:,max(z1)+1:end) = nan;
mask1(~isnan(mask1)) = 0;   % set those points that are inside the image to 0 (in image = 0, out image = nan)

% convert the mask points that are in the image to a grid, which is a rectangular prism at this point
[cx, cy, cz] = ind2sub(size(mask1), find(mask1 == 0));
marr = [cx cy cz];

% Create 3D alpha shape for the image data and find points in marr that are within the shape
shp = alphaShape(x1,y1,z1,48);           
tf = inShape(shp,marr(:,1),marr(:,2),marr(:,3));          % THIS IS THE SLOW STEP

% Crop marr to the positions that are inShape and populate mask1 with 1s.  Yield is mask1, and m x n x p array with 1s where the image is.
marr = marr(tf == 1,:);
for i=1:sum(tf)
    mask1(marr(i,1),marr(i,2),marr(i,3)) = 1;
end
mask1(isnan(mask1)) = 0;
figure

hold on
scatter3(x1,y1,z1,'.');
plot(shp,'FaceAlpha',0.2)
hold off
title('Original points and shp')

figure
hold on
scatter3(cx,cy,cz,1,'.');
plot(shp,'FaceAlpha',0.2)
hold off
title('Cropped rectangular mask')

figure
hold on
scatter3(marr(:,1),marr(:,2),marr(:,3),'.')
plot(shp,'FaceAlpha',0.2)
hold off
title('End product after inshape')

Here's the graph outputs (this is the xy view)

Answer 1

[raymeout]

Have you looked at the function pointLocation? But to use this function, you'll need to create a tetrahedral triangulation first (e.g. tri = alphaTriangulation(shp)).

pointLocation will find the tetrahedron, which encloses the specific points → if a point is in none (→ outside), it will get NaN...

I did not do any speed benchmark, but in the cases I used the function, it seemed to be pretty fast, even for bigger point clouds and tetrahedral meshes.

[GolgiWhillikers]

I tried this out and a few things-

1) alphaTriangulation doesn't work for this, since pointLocation requires a traingulation object and won't just take the vertices that come from the alphaTriangulation function. But...

2) Either way, it seems to be the same speed, at least with the point cloud I tested. In one specific instance of my code above, I'm getting about 0.4-0.41 vs 0.38-0.39 seconds for alphaShape/inShape vs delaunayTriangulation/pointLocation. So it's a very similar speed. I suspect that any performance increase is that delaunayTraingulation is a bit faster than alphaShape, while pointLocation and inShape take a similar amount of time. I'm not a code person but I'd guess this is just running into how fast search algorithms are, since it's essentially the same ask either way.

Thanks for the suggestion.

[raymeout]

alphaTriangulation returns the full triangulation of the alphaShape... You can create the required object by e.g.:

tri_object = triangulation(tri, 3DPoints);

But either way: if it isn't faster, it won't be of any use in this case...

(A possible reason, why your code with inShape might be a bit slower compared to the pointLocation-version could also be, that inShape might use pointLocation and alphaTriangulation internally anyways (since pointLocation was introduced earlier). But I didn'd check this statement any thurther...)

[GolgiWhillikers]

Ah, thanks for the object conversion code, I had tried to find something like that but couldn't.

And yeah, I figure that's what's happening for the speed. Thanks.