Get regular grid and points of a given stl file

Hello,

I have have a stl file, which has many points (ca. ten thousand) but they are not uniformly distributed. As you know, a stl object has Faces, Vertices and Normals. The Vertices are the points (they can be repeated a few times, depending on how many triangles you can form with them).

I would like to create points in between the vertices (I don't care if I lose these initial points) so that they are evenly distributed. Imagine I project the object (just the points) in a 2D plane. I would like the points to be regularly distributed (like a grid but without the lines).

[KSSV]

But why?

[Diego Hens]

I would like to do this for multiple images. I want to represent the height of the obejct in a two-dimensional plane with a grayscale value. Then I can compare these values between objects.

(Additionally, this two-dimensional image will look better than the object with its original points, as there are holes in this representation)

[KSSV]

Any pictorial example?

[Diego Hens]

sorry, I answered on the wrong spot. Can you take a look on my other comment?

Answer 1

[Bruno Luong]

I load the stl, cleanup the vertices, compute the face normal then I extract the top cap and interpolate on grid.

% Read mesh and compute normal vector
TR = stlread('Zahn_15_Bibliotheksmodell.stl'); % contains "Mesh" structure
[X,~,J] = unique(TR.Points,'rows');
F = J(TR.ConnectivityList).';
X = X.';
XF = reshape(X(:,F),3,3,[]);
N = cross(XF(:,2,:)-XF(:,1,:),XF(:,3,:)-XF(:,2,:),1);

% Extract the points on the top cap
Fup = F(:,N(3,:)>0);
Fup = sort(Fup,1);
n =  max(Fup(:));
A = accumarray([Fup([1 2],:), Fup([1 3],:), Fup([2 3],:)]', 1, [n,n], [], 0, true);
G = graph(A+A');
bins = G.conncomp;
n = accumarray(bins(:),1);
[~,icap] = max(n);
icap = find(bins==icap);
Xcap = X(:,icap);
x = Xcap(1,:);
y = Xcap(2,:);
z = Xcap(3,:);

% Interpolation on grid
xmin = min(x); xmax = max(x);
ymin = min(y); ymax = max(y);
nx = 513;  ny = 513; % define resolution of grids
xi = linspace(xmin,xmax,nx);
yi = linspace(ymin,ymax,ny);
[XI,YI] = meshgrid(xi,yi);
I = scatteredInterpolant(x(:),y(:),z(:),'linear','none');
ZI = I(XI,YI);

% Plot the result
figure(1);
ax=subplot(1,2,1);
imagesc(ax,xi,yi,ZI);
set(ax,'Ydir','normal','clim',[28 31])
colormap(ax,gray);
axis(ax,'equal');
axis(ax,'tight');
ax=subplot(1,2,2);
surf(ax,xi,yi,ZI,'linestyle','none');

There are some artefact, and did not try to see why, but it gives you an idea.

[Diego Hens]

Thank you so much for taking the time to help me.Really.

I'm going to take a careful look at it and understand what you did.

I'll come back then :)

[Diego Hens]

Hello Bruno,

I haven't understood the important parts of what you've done. I'm sorry, I tried.

I haven't understood how you created the mesh structure and thus I cannot replicate this code with other objects. The struct has 10 fields and I don't know what some of them are (for example Boundary, or Boundary.E).

I also haven't understood how you extract the points on the top cap.

[Bruno Luong]

I edit with code to create mesh F and X and N manually

[Diego Hens]

I'm a bit slow, but I think I've got it now (not really the code, but the problem in itself). Your answer was actually what I looked for (the part where you extract the points on the top cap I did not understand, but I could use the rest still). I think I can go on from here.

I have accepted it as an answer for this reason.

I will try to find out why there is noise and if I can avoid it. If not, I will look for other methods.

Thank you again, Bruno :D

Answer 2

[Diego Hens]

Ok, here's what I got.

Imagine the triangle on the left side was my object. It has values in X, Y and Z. I want to project the object into the XY plane (as in the right) and represent the Z value with a grayscale value. If I project the object I get the drawing ob the right. As I have only three vertices, I will only get three points with a grayscale value. That means (especially if the object is bigger) that I will have many holes, where no value Z is given, as there is no point.

I would like to be able to interpolate these values in between, so that I have no (or very little) holes.

Is the problem understandable now?

Thanks for the help :D

[Bruno Luong]

If you upload one of your typical STL file I might take a look.

[Diego Hens]

But the file should not matter, right?

Either way, I have uploaded one of them.

[Bruno Luong]

I take a quick look of your object; it doesn't seem one can represents such object as a 2D function z=f(x,y), so nothing like the triangle you have illustrated.

So sorry it's not clear to me what you want to do with projection and resampling.

[Diego Hens]

This is an example of what it would look like. This is a tooth seen from above (it's from a different tooth, but it should not matter). The "height" is represented by a grayscale value, dark being more in depth.

I'll go into much more detail to describe how I got it:

1) First transformed the STL into a pointcloud, removing the repeated points using the command unique (x, 'rows') and pointCloud(x)

2) I ran a pca (principal component analysis) to get the coefficients (principal vectors) to draw the axis of rotation for another problem, and the scores, which are my new values in X, Y and Z.

3) I name the first column X, the second Y adn the third Z and plot(X,Y) with this extra

scatter(x,y,25,z,'filled') ;
colormap(gray)
axis('image');
set(gca, 'color','k');

Now if you look at the image, you see the holes I was talking about. If I had more regular points they wouldn't be there. Of course I can fill the holes by making the markers bigger, but thats not the point. I would really like to know what the interpolated values at every point would be.

I'm sorry if I made you have the feeling like not explaining the problem detailed enough. I thought the more details I give, the more complex the question would be and the fewer answers I would get. I really thought a more simplified approach was better.

[Bruno Luong]

Have you tried

scatteredInterpolant

[Diego Hens]

I think that's exactly what I want. This is what I get:

F = scatteredInterpolant with properties:

Points: [25634×2 double]

Values: [25634×1 double]

Method: 'linear'

ExtrapolationMethod: 'linear'

How can I plot this now?

[Bruno Luong]

The difficult part is you need to filter the points and keep only the point on the cap of the tooth. Because - as I said - your data is not of the form z=f(x,y) and there are multiple z for the same (x,y) (close to the boundary if one looks from the above). This will polute the interpolation.