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!

[Mathieu NOE]

hmmm... there are a couple of side effects when you create an interpolated surface from a cloud of 3d scattered data - IMHO those points do represent a volume rather than a surface so I have 2 comments here

• some contour lines (red arrows) are not meaningfull as they are simply the effect of the interpolated surface having an uncontrolled shape outside the x,y area in which we have measurement data
• even inside the raw data volume, one can ask himself if the shape of the interpolated surface is meaning full , as this surface should permantenly jump from positive Z to negative Z data (of very close similar x,y points)

here a code that overlays the raw data points with the interpolation surface , NB the points and the surface may not match (i.e. the surface does not go exactly through the points , it's almost impossible in some areas where z make a lot of zig zags)

NB I don't apply axis equal because it tends here to squeeze the plot and does not help in my modest analysis of the situation

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);
hsc(1)
set(hsc(1),"EdgeColor","none","FaceAlpha",0.5);
% shading interp; % Smooths the surface
grid on
colormap(jet)
colorbar


%% adding my simple 3d dot plot / colored by slices 
hold on
zmin = min(z);
zmax = max(z);
levels = 25;
zt = linspace(zmin,zmax,levels+1);
colors = jet(levels+1);

for k = 1:levels
    ind = (z>=zt(k) & (z<zt(k+1)));

    plot3(x(ind),y(ind),z(ind),'*',"color",colors(k,:));
end
view(3)

[Mathieu NOE]

then I tried to better visualize the data - here's waht seems to work best using MyRobustCrust :

Surface Reconstruction From Scattered Points Cloud - File Exchange - MATLAB Central

code :

%% Load the data
load('workspace.mat');

[t,tnorm]=MyRobustCrust(xyzDATA);

%% plot of the output triangulation
figure(1)
title('Output Triangulation','fontsize',14)
p = trisurf(t,xyzDATA(:,1),xyzDATA(:,2),xyzDATA(:,3));
colormap("jet")
p.EdgeColor = 'interp';
p.FaceColor = 'interp';

view(-40,24)
box on
camlight(40,40)                                
camlight(-20,-10)

results : it looks like an irregular slice of pizza with a missing portion ....

seen from the top :

seen from the bottom :

[Mathieu NOE]

of course it's fairly simple to create N slices along the Z direction, but the resulting plot for each slice must also be interpreted with some precautions - what do we do next ? remove outliers ? take the boundary ?

example with 10 slices (demo) :

load('workspace.mat');
x = xyzDATA(:,1);
y = xyzDATA(:,2);
z = xyzDATA(:,3);

% plot individual levels 
zmin = min(z);
zmax = max(z);

levels = 10;
zt = linspace(zmin,zmax,levels+1);

% figure
% hold on
for k = 1:levels
    ind = (z>=zt(k) & (z<zt(k+1)));
    figure
    plot(x(ind),y(ind),'*');
end

Warning: Hardware-accelerated graphics is unavailable. Displaying fewer markers to preserve interactivity.

Warning: Hardware-accelerated graphics is unavailable. Displaying fewer markers to preserve interactivity.

Warning: Hardware-accelerated graphics is unavailable. Displaying fewer markers to preserve interactivity.

[Mathieu NOE]

you could also split your volume in two halves with positive and negative z surfaces , then do a contour plot - as we have scattered data , use tricontour : Contour Plot for Scattered Data - File Exchange - MATLAB Central

to extract the contour lines data , you can use C2xyz : C2xyz contour matrix to coordinates - File Exchange - MATLAB Central

example for positive z data - same principle is doable for negative z data

load('workspace.mat');
x = xyzDATA(:,1);
y = xyzDATA(:,2);
z = xyzDATA(:,3);

% contour plot for positive z data 
ind = (z>=0);

figure
tri=delaunay(x(ind),y(ind));           % triangulate scattered data
[C,h]=tricontour(tri,x(ind),y(ind),z(ind),10);
title Positive Z Contour
colormap("jet")
colorbar("vert")

% Given a contour plot, you want to know the (x,y) coordinates of the contours, 
% as well as the z value corresponding to each contour line.
[xc,yc,zc] = C2xyz(C);

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.