Efficient display of 2D images in 3D space?

How can I plot a set of 2D ground penetrating radar sections in 3D, as rendered 2D images that have the appearance of 3D slices? [An analogy might be that I have a set of paintings, which I want to stack in a rack. The rack creates a 3D geometry for viewing the paintings, but the canvases themselves are still 2D.]

I have data from a series of parallel radar scans through a highway. Each scan generated a time series that can be arranged as xyzA where in a simple case: x is distance along the line, y is a constant for each line, z is depth, and A is the radar signal's amplitude. A typical scan comprises ~3 million data points.

A single scan can be displayed in 2D using imagesc (using the x,z, and A components). My task is to visualize the data in 3D. In theory, I could interpolate all of the scanned data to make a 3D volume, and slice it, but that will not differentiate between noise and signal, and is likely to create persistent artifacts. Instead, I'd like to just display a series of 2D slices in 3D space, where each slice is an individual scan.

This sounds like a trivial problem (after all, a 3 million point slice is equivalent to a low-resolution photo, and 3D perspective computer games handle this kind of problem in real time). So far, though I've had no success. The prospective solutions that I've found online can't handle the number of points involved, and produce memory errors. I've tried a couple of methods so far:

- 3D scatter plotting where each plotted point represents one data point from the file, and its colour is controlled by the signal amplitude; and - A slicing method that I found via Matlab code sharing. This used meshgrid to generate a distribution of data across a 3D plane, but could not cope with the number of points to be gridded.

I'm wondering whether the solution lies in the way the data are being managed. Would some approach involving projection of a jpeg onto a 3D surface be more efficient than trying to plot each value from a matrix of 3D data points?

An alternative is to go back to the start, and say "OK. If my data were in a commercial radar format, I could just use dedicated radar processing software. So then how can I flip the file from its current proprietary format to be readable by such a package?" I'll do that if this avenue turns out to be a dead end, but for now I'm looking for a Matlab visualization solution.

Answer 1

[Mike Garrity]

I wouldn't think that going to a scatter plot would be a good approach here.

What about something like this?

img = imread('ngc6543a.jpg');
for z=-100:10:100
  g = hgtransform('Matrix',makehgtform('translate',[0 0 z]));
  image(g,img)
end
view(3)

[Prairie geophysics]

Thanks Mike for your rapid reply! I'll have a go at adapting that approach and let you know how it works out.

[Prairie geophysics]

I'm having some trouble with that code, although I'm trying to work through it. Unchanged apart from the filename, it returns:

Error using image Incorrect number of arguments.

Error in SANDBOX (line 4) image(g,img)

I can see that the code creates a transform that would step through a series of xyz coordinates for the image, but it does not appear to apply that transform for each value of z to the image itself. Passing the transform to image() causes an exception; I think image() is looking for coordinates rather than the transform value. Nevertheless, the stacked images are a rotated version of what I'm after, so thanks - I'm just not sure how the code as it is posted could have generated them. Maybe I'm not implementing it correctly?

[Mike Garrity]

Sorry about that. I think that syntax was introduced in R2016a. An equivalent which has been around forever is this:

image(img,'Parent',g)

[Prairie geophysics]

Hey Mike! I must admit that I am still trying to understand the syntax, but like a monkey with a typewriter and plenty of time, I got something basic to work, and can expand it to do everything I want it to. I thought I'd share it with you, and anyone else who might be trying this in the future.

I am using photos for now instead of data - it is much easier to tell whether the translation and rotation are behaving properly this way. One radargram looks rather like another...

The code reads two photos in, then creates transforms to rotate them both, and translate one of them. It runs in less than 0.2 seconds, and the rotation/zoom tools in the viewer are real-time. That's exactly what I was looking for! Thanks for your help. Warning: ropey novice code follows...

if true
  % code
tic
imagepath = 'C:\\Matlab\\Radar\\Sandbox';
cd (imagepath);
ax = axes('XLim',[0 2000],'YLim',[-1 1],'ZLim',[-1000 0]);
ax.XLabel.String = 'X';
ax.XLabel.FontSize = 12;
ax.YLabel.String = 'Y';
ax.YLabel.FontSize = 12;
ax.ZLabel.String = 'Z';
ax.ZLabel.FontSize = 12;
view(3)
img = imread('towers.jpg'); %read in photo1 - Porcupine Ridge
img2 = imread('MattDM.jpg'); %read in photo2 - Prairie Mountain
t = hgtransform('Parent',ax);

   q = pi()/2;
   r = [0 0 0.5];
   h = hgtransform('Matrix',makehgtform('xrotate',-q));
   h2 = hgtransform('Matrix',makehgtform('xrotate',-q,'translate',r));
   set(h,'Parent',t)
   set(h2,'Parent',t)
image(img,'Parent',h)
  image(img2,'Parent',h2)

toc end

[Mike Garrity]

That's great.

Yeah, the syntax is a bit weird. The basic idea is that you're building a tree structure of objects. The Axes has two children which are HGTransforms, and each of those has one child which is an Image. The point of this is that the HGTransform object has a Matrix property which applies to all of its children. This allows you to move things around, just like changing the Limits of the Axes would. It's just that adding a couple of HGTransform objects gives you more knobs that you can play with.

[Walter Roberson]

Guojin Feng comments to Prairie Geophysics:

Good solution

[Con Ahern]

I combined this solution with imagesc to scale the image I have to "paint" walls (using multiple rotations and translations) with a certain marble described by a JPEG file. Thanks very much!

Answer 2

[Mohamed Amine Henchir]

Hey Mike, Thanks for your solution! I was able to use it to show a stack of images updating every iteration. however I still have another challenge which is to plot detected circles in these images along with the original images on the same stack.

here is my code

    [centersBright1, radiiBright1, metric1] = imfindcircles(I,[Rmin Rmid ],'ObjectPolarity',...
        'bright','Method','PhaseCode','Sensitivity',slider2Value,'EdgeThreshold',slider1Value);
    [centersBright2, radiiBright2, metric2] = imfindcircles(I,[Rmid Rmax],'ObjectPolarity',...
        'bright','Method','PhaseCode','Sensitivity',slider2Value,'EdgeThreshold',slider1Value);
    [x,y,r] = sort_cirlces(centersBright1,radiiBright1,metric1,centersBright2,radiiBright2,metric2,Rmin,Rmid);
    center = [x y];
    viscircles(center,r,'EdgeColor','g');
    h = hgtransform('Matrix',makehgtform('translate',[0 0 p*10]));   
    imagesc(h,I)
    colormap(gray(65536));

here is a picture of the current output, I would like to have the circles following its corresponding picture, but for some reason they keep being updated at the bottom.

any Idea how to fix this? thank you very much!