help with XYZ data file

I am trying to understand a XYZ datafile from a 3D file. Example file.

https://www.dropbox.com/scl/fi/9463jgix2t8bnz356s7ly/1-0331_06_mesh.xyz?rlkey=m73j8q1j9qqufe7xlqbockg68&st=vlsz2eoz&dl=0

I was able to google how to open the XYZ file in matlab and it gave me code that works (provided below) but I would like help manipulating the code. First question, I would like to be able to plot some of the planes of XY with fixed Z values/range. Can this be done direclty from the MATLAB variables XYZ used below? This seems to be making lines, it is some kind of vector data?

figure, imagesc(X(:,1),Y(:,1),Z(100:1000,1))

What is the scatterinterpolant doing? How can i get XY planes for each Z ?

clear
data = readmatrix('#1-0331_06_mesh.xyz', 'FileType', 'text');
X = data(:,1);
Y = data(:,2);
Z = data(:,3);

% Create a grid
[xq, yq] = meshgrid(linspace(min(X), max(X), 100), linspace(min(Y), max(Y), 100));
% Interpolate Z values onto the grid.  Different interploation methods

F = scatteredInterpolant(X, Y, Z);
zq = F(xq, yq);

% Plot the surface
figure, surf(xq, yq, zq);
shading interp; % Smooths the surface
axis equal

[Matt J]

Please attach the X,Y,Z data in a .mat file so we can demonstrate solutions.

[JM]

The data of the XYZ file seems to be a vector. The "scatterinterpolant" seems to be plotting points and then filling it in to make it look good in the figure, from what I can tell.

The main problem is that we need to get the data into a 3D image array as we want to overlay it onto a 3D image set.

The goal is to have a image(x,y) at each discrete step of height (z) in the image set the same way the other 3D image sets do this

[JM]

this is the .mat file of the data,. downsampled

I thought someone asked for it

data = readmatrix('0429_01_mesh rev3.xyz', 'FileType', 'text');

xyzDATA=data (1:2:end,1:2:end,1:2:end); % smaller XYZ file (1/8X)

clear data

[JM]

so we can do a contour of the surface

we can do it at each heigth level

Can i try to rephase the questions

can we create an array of the coutours at each height of the contour.?

so the highest point on the countor would be the largest values in the 3D image set (largest z value)?

the smallest countour value would basically have all teh coutour lines and be represented by the lowest z value?

[Matt J]

The main problem is that we need to get the data into a 3D image array as we want to overlay it onto a 3D image set.

A 3D image is defined by quadruplets (i,j,k,v) where (i,j,k) are voxel coordinates in the 3D image and v is the voxel value at (i,j,k).

Right now, you don't have quadruplets, but rather triplets (x,y,z). It is not clear what mapping

(x,y,z)--->(i,j,k,v)

you are trying to achieve.

Answer 1

[Star Strider]

You may only need to use reshape (with the same last two arguments) on each of the vectors, rather than using scatteredInterpolant. There's nothing wrong with using scatteredInterpolant, and it may have its uses here, depending on what you want to do, however that may not be the most efficient initial appropach.

Use the zip function to create a .zip file for the .xyz file so you can upload it here (providing it meets the 5MB size limit).

[JM]

the xyz file is uploaded in the first post (I posted it on dropbox). It is 12mb.

its really appreicated if you can show how to use reshape.

[JM]

its posted below, thank you

[Star Strider]

I've had no end of problems with Answers this afternoon.

On my laptop here. Desktop seems to need a reboot.

The reshape function requires gridded data for problems like this. Your data are scattered, so reshape isn't gong to work.

This works in MATLAB Online --

load('workspace.mat');
X = xyzDATA(:,1);
Y = xyzDATA(:,2);
Z = xyzDATA(:,3);

% Create a grid
gridsize = 1E+3;
[xq, yq] = meshgrid(linspace(min(X), max(X), gridsize), linspace(min(Y), max(Y), gridsize));
% Interpolate Z values onto the grid.  Different interploation methods
F = scatteredInterpolant(X, Y, Z);
zq = F(xq, yq);
% Plot the surface
figure
hsc = surfc(xq, yq, zq);
% get(hsc(2))
Levels = hsc(2).LevelList

Levels = 1×9

-20 -15 -10 -5 0 5 10 15 20

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

shading interp; % Smooths the surface
axis equal
grid on
colormap(turbo)
colorbar

M = hsc(2).ContourMatrix;
for k = 1:numel(Levels)
    idx = find(M(1,:) == Levels(k));
    ValidV = rem(M(2,idx),1) == 0;
    StartIdx{k,:} = idx(ValidV);
    VLen{k,:} = M(2,StartIdx{k});
end

cm = turbo(numel(Levels));

figure
hold on
for k1 = 1:numel(Levels)
    for k2 = 1:numel(StartIdx{k1})
        idxv = StartIdx{k1}(k2)+1 : StartIdx{k1}(k2)+VLen{k1}(k2);          % Index For Contour 'k1'
        xv = M(1,idxv);
        yv = M(2,idxv);
        hp{k1,k2} = plot(xv, yv, Color=cm(k1,:), DisplayName=sprintf('Level = %2d',Levels(k1)));
    end
end
hold off
axis('equal')
grid on
xlabel('X')
ylabel('Y')
title('Recovered Contours')
legend([hp{:,1}], Location='best')

Warning: Graphics acceleration hardware is unavailable. Graphics quality and performance might be diminished. See MATLAB System Requirements.

I'm still not certain what you want to do. This recovers the contours at each level, and plots them.

.

[JM]

wow, that is pretty amazing, I spent a few hours analyzing and using this.

This goal is combine laser TOF height measurements (ie lidar) with 3D xray image data created from a 3D imaging sytem that has a limted scanning angle. With limited angle systems, there are often reconstruction artifacts that protrude out of the rendering of the surface of the object. By having a set of contours in a cartesian coordinate systems from the laser scanning data, we can create a set of image masks to apply to the 3D xray data. This can be used to remove the artifacts in the 3D data.

It looks like there are 2 loops, one for the number of levels set by levellist (k1 = 1:numel(Levels)) and the other is for the number of coutours at each level (k2 = 1:numel(StartIdx{k1})). Is there a way to get all the coutours at each level into a single array and turn it into a mask? That way we have a mask to apply to each of the height levels in the 3D image.

[JM]

think need to implement something like

mask = poly2mask(x_coords, y_coords, size(Z, 1), size(Z, 2));

imshow(mask);

[Matt J]

Please do not post your remarks as an Answer when they are a response to somebody. Use Comments instead, otherwise we cannot tell who they are directed to. The Answers boxes are meant for actual answers to the posted question.

[JM]

ok, thank you

[Star Strider]

Thank you!

'Is there a way to get all the coutours at each level into a single array and turn it into a mask?'

The second loop (creating the last figure) gets all the contours at each level and puts them into the respective 'xv' and 'yv' vectors. (You can store those as sparate variables in 2D cell arrays to use later.) They can either be single lines, or closed, depending on how contour creates them. If they are single lines, you would have to artificially create closed contours by supplying the necessary edge regions. (So for example at level +20, you would have to create a 'x' line from 0 to the right end of it, and a separate 'y' line from 20 to the low end of it.) You can test that by comparing the first and last coordinates, If they're the same, the contour is closed, if they are not, the contour is open.

Here, the 0 level (pale green in the last figure) has several different contours (at least 10, probably many more than that since I didn't specifically have them counted here), so combining them all into one mask could be a challenge. You can change the number of levels and also specify their values, so those depicted here are not the only options.

Additionally, you can change the 'gridsize' value. I could increase it to 1E+5 in MATLAB Online, however that caused problems here, so I settled for 1E+3. Increasing it increases the resolution of the interpolation.

This is the best I can do at this point.

[JM]

amazing, thank you

[Star Strider]

As always, my pleasure!

Answer 2

[Matt J]

First question, I would like to be able to plot some of the planes of XY with fixed Z values/range.

The terminology "planes of XY" is not entirely clear to me, but I suspect that you want,

contour(xq,yq,zq)

In the documentation for contour, you will also see that there are further options to specify the precise isocontour z-levels that you want.

[JM]

rather than coutour at z levels, is there a way to get fixed data such as a imagesc(x,y @z)

[Matt J]

I'm afraid you're not going to be able to make your meaning understood by using fake Matlab command syntax like imagesc(x,y @z)

Right now you have a surface z(x,y). A surface and an image are not the same thing. If you take the intersection of a surface and a plane, you get a curve, not an image.

If x,y are integer-valued, and are supposed to represent pixel coordinates in some image, and you want to assign a value to each pixel (x,y), you must explain how the pixel value at each x,y is to be computed from (x,y,z).

Answer 3

[dpb]

"What is the scatterinterpolant doing?"

Creating a set of interopolated z values, zq, on a regular grid xq, yq over the range of X, Y (100 points in each direction).

"How can i get XY planes for each Z ?"

nZ=size(data,3);        % the number of planes in data array
for i=1:nZ              % iterate over all planes
  Zi=data(:,:,i);       % retrieve the ith plane
  % do whatever wanted here...
  
end

NOTA BENE:

MATLAB has regular arrays but the data in X, Y may not be regularly spaced (hence the use of the scatteredinterpolant() above, probably).. The above retrieves the actual data points only at whatever are the specific locations given.

[JM]

thank you

[JM]

if you dont mind, can you provide a more detail on the code example

should I call out coutour(x,y) at each z limit and then reshape it into 3D array?

Answer 4

[Matt J]

This might be what you're looking for:

zq=fspecial('gaussian',200,40);

Nlevels=5;

zRanges=linspace(min(zq(:)) ,max(zq(:)), Nlevels+1);

levelMap = discretize(zq,zRanges);

Masks = levelMap==reshape( 1:Nlevels,1,1,[]);

tiledlayout('h');
nexttile; imshow(zq,[]); 
nexttile; imshow(levelMap,[]); 
nexttile; imshow(Masks(:,:,3).*zq,[]); set(gcf(),'Position',[34 487 1529 463])

[JM]

This seems like it would work, Im was not able to get it to work though. I do appreicate the suggestion.