sparse to gridded 3-D interpolation

5 views (last 30 days)
B Huey
B Huey on 14 Nov 2019
Edited: John D'Errico on 14 Nov 2019
I regularly need to ingest 3-D datasets of an image Intensity signal that is acquired in such as way that the data is initially "sparse" (irregularly acquired) in the Z direction. To then analyze this data using conventional image analysis, stacking routines, etc., this data needs to be interpolated into gridded voxels which are regular along Z.
I have found that using the matlab command scatteredInterpolant provides excellent results, especially when using the method and extrapolation tags of 'natural' and 'none' respectively. Normally matlab calculates the 3DInterpolant in 1-10 seconds. No problem...
HOWEVER, when evaluating the 3Dinterpolant at gridded datapoints, this process takes hours to days.
Is there a way to accelerate the evaluation?
PLEASE HELP!
Details:
For instance, I need to convert data that is intially of the form (xregular,yregular,zarbitrary,Intensity) into output image stacks of (xregular,yregular,zregular,Intensityinterpolated). Typical input datasets include ~512regular * 512regular * 100sparse datapoints, with at least 8bit but ideally 16bit resolution for the Intensity channel at each of these 512*512*100 points.
My outputs need to have roughly (though not always) the same number of voxels (eg 512x512x150), but gridded in XY AND Z instead of sparse along Z.
Maybe it's helpful that the XY vectors are regular, though this may be irrelevant since Z is NOT regular.
Hardware: I have a brand new, pentium I9 desktop, with 16GB ram and a serious GPU, running Win10pro and matlab 2019b. But matlab seems to only use 10% of the CPU capability during this calculation, and only half of the overall memory, and none of the gpu. It ought to be MUCH faster.

Answers (1)

John D'Errico
John D'Errico on 14 Nov 2019
First, should MATLAB use all of the CPU capability? That is difficult in what you are doing, becaue it is not easily parallelized. The scatteredInterpolant is doing its work using a 3-d tessellation. Each point will lie in one simplex of the tessellation. So it needs to decide where a point lies, then interpolate inside that simplex. What happens is this is not necessarily easy to do in a way that uses all of your cores.
Some things in MATLAB that use linear algebra or some other operations are automatically parallelized. So it can be the case that you will see all of your cores running flat out. (Yes, I checked the other day on something that I was doing, and it was nice to see MATLAB running at 800%, thus all 8 cores on my machine, flat out. But then the easy part finished, and it dropped down to 100% of only one core, since the next step was not amenable to multi-core processing.) Perhaps a future release will see TMW include multi-thread processing in scatteredInterpolant, but that seems to be not the case today. They have only so many resources to apply to such enhancements.
(A quick check shows that scatteredInterpolant seems to use only one core on my 8 core hardware, even on a large problem. Done in R2019a.)
Next is the question of whether you are even doing this the correct way.
It seems from your question that you have a regular grid in the (x,y) plane already. So all points lie on that grid. But then you have points in z for each (x,y) pair, that are NOT regular. Your goal is to interpolate these points onto a 3-d lattice in (x,y,z), and the way you are doing it is to just throw it into scatteredInterpolant.
Effectively, it seems you are using a 3-d interpolation in a way that really needs only a 1-d interpolant.
My suggestion is to just reduce this into a 1-D interpolation, although many of them. At any (x,y) location, interpolate along z to grid the result. Each such interpolation is trivial and fast. It will run like blazes, even in a loop over x and y.
Don't use a high horsepower tool to solve a simple problem.
  2 Comments
B Huey
B Huey on 14 Nov 2019
I appreciate the speed enhancement if I can turn the problem into numerous 1D interpolations. But this is 3D data after all, and points which are near each other in xyz space should have similar values. Therefore this isn't Nx*Ny independent 1D problems--it really is a 3D problem.
Any other suggestions for how to better get to an answer?
And yes, leveraging all 16 cores would obviously make me much happier (you are right, matlab is maxing out 100% of just one of the cores). It would be great if mathworks prioritized this for a future update--there's a lot of horepower which is not being utilized.
-Bryan
John D'Errico
John D'Errico on 14 Nov 2019
Edited: John D'Errico on 14 Nov 2019
Some of what is involved can be easily parallelized. But not all. Actually building a tessellation in n-dimensions is probably somewhat difficult to parallelize. But then interpolating multiple points can be set up fairly easily, essentially a parfor loop, where each point to be interpolated uses the same global data, but doing different operations, picking out the appropriate simplex that contains the point, then interpolating inside that simplex.
It may SEEM like a 3-d problem to you. However, your intuition can easily lead you astray here.
If the triangulation scatteedInterpolant creates always have edges that are vertical in z (valid to some extent, but not always true) then the interpolant will purely be a function of z in those regions. Let me give an example. 2-d will suffice.
[x,z] = meshgrid(0:5,0:6);
z = z - z.*(6.25 - (x - 2.5).^2)/10;
plot(x,z,'ob')
So, in two dimensions, I've created a set that has fixed values of x, but the z spacing varies as a function of x.
Now, suppose I want to establish a nice regular grid in both x and z? Can I use scatteredInterpolant? What will happen if I do so?
tri = delaunay(x(:),z(:));
hold on
trimesh(tri,x(:),z(:),'color','k')
Now suppose I wanted to use scattered interpolant to predict at coordinates along a vertical line at x==3?
From z == 0, the interpolant will be exactly along an edge of one of those triangles. As such, a linear interpolant in such a scenario will be purely a function ONLY of the values at the end points of the line on that edge of the triangle. So, for x==3, when z is below roughly 2.5, the interpolant will be purely a function ONLY of the z value. You coule literally use a 1-d interpolant there, up to that point. A problem arises when z exceeds 2.5 though. Now it will use those very large triangles to interpolate for z>2.5. In fact for z larger then 4.5 or so, the triangle that will be used is huge. The triangle used there will have vertices at roughly (0,6), (4,3.5), and (5,6).
What I am trying tp point out are two problems with your intuition.
First, for some values of z, the interpolant will be purely a function of z. It will not remotely be a function of the nearby points. But at some point, the interpolant will need to start using other points, but the problem is? The other data points used will not be close at all! This is a common artifact of using linear interpolation in a tessellation to extrapolate outside of a non-convex region in multiple dimensions.
What happens from a linear interpolant is you get potential garbage in those regions. You may well have done as well to just extrapolate an interpolant that was purly a function of z, as I suggested.

Sign in to comment.

Categories

Find more on Parallel Computing Fundamentals in Help Center and File Exchange

Products


Release

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!