Radius of curvature of a crater on a surface

I need to be able to use MATLAB to measure the radius of curvature of a "dimple" on a surface. I cannot post the images here, but imagine a nice rounded surface with a dimple or crater in it, as if something impacted it and left a forever dent.

I've seen some code on here regarding radius of curvature, but I cannot make it work for what I have. My images are imported as "STL" files, and as of now, I can only plot them using the scatter3 command, meaning I have a surface, but not a smooth one... one with scatter plot dots everywhere. You can easily visualize it, though.

Is there a way to do this? If ML can find the slope change where the surface stops and crater starts, it can give me the points for the "circle" to use, and I'm sure it can measure the depth, but as for the radius of curvature of that hole, I'm at a loss.

Any ideas?

[Sean de Wolski]

First you have to define the surface, i.e. the connectivity of the vertices. Do you have this information from the STL file?

[John]

When the STL file is imported, it dumps data into a two-part structure: faces and vertices.

Answer 1

[Sean de Wolski]

So you'll need to first determine all of the faces that are within the dimple. If the rest of the mesh is convex, you could do this by looking at vertices that are not part of the convex hull and all of the faces that contain these vertices. Once you have this information, calculation of the radius of curvature just becomes algebra.

For calculating the convex hull:

http://www.mathworks.com/help/releases/R2013b/matlab/ref/delaunaytriangulation.convexhull.html

[John]

That link won't open for me. Do you mind copying and pasting the instructions and/or code?

[Sean de Wolski]

Unfortunately I cannot. You can just type the following at the command prompt and then select the "convedHull" method:

doc delaunayTriangulation

[John]

Please bear with me here... I'm fairly new to MATLAB.

I understand what you saying, but I can't get it to work in ML. I know I need to get the points of the bound circle of the dimple, but I'm not connecting what you're suggesting to ML.

When importing the STL file into ML, I get a file structure that has the faces (roughly 2000x3 matrix) and the vertices (roughly 5000x3 matrix). I can use scatter3 to plot the vertices telling x, y, and z to be the first, second, and third column respectively of my vertices. It shows me the shape.

I found the delaunay Triangulation of the x, y, and z from above and used that to find the convex hull of that triangulation, and it returned a much smaller matrix (roughly 1000x3), but I don't know what those values mean, and I don't know what to do with them. I also don't know what to do going forward.

Any help is appreciated.

[Sean de Wolski]

Do you have a sample or dummy STL file you could give us?

[Sean de Wolski]

%Sample random data
fv = isosurface(randn(50,50,50),2);
faces = fv.faces;
verts = fv.vertices;

%triangulate
dt = delaunayTriangulation(verts);

%vertices on the hull
onHull = convexHull(dt);
k = dt.Points(onHull);

Now you have the points on the hull, remove them, and you'll have the remainder.

[Image Analyst]

What program measured your item and created the file? Can you have that program export the data into some kind of 2D image where the pixel value is the height information? If not a PNG file, then some kind of standard image file format or even a raw/proprietary format as long as we know what the format is.

[Sean de Wolski]

Unfortunately an internet search does not turn up any useful results for "STL file death star" :)

[John]

Sean and Image Analyst...

You are correct, no STL Death Star can be found. Since I cannot share the image I'm using (proprietary), I don't know exactly how to go any further. A google search for something similar pulls this up. My image is an image of a helmet with damage to it... imagine the helmet was impacted by some blunt object and it left a dent.

As for what program created the image I have, I don't know. I didn't create it, it was sent to me, so I can't recreate any 2D image for you.

[John]

// Now you have the points on the hull, remove them, and you'll have the remainder. //

May be a silly question, but what exactly does this comment mean? The remainder is the dimple I'm trying to find the radius of curvature of? And what exactly does "remove" them mean? There are almost 6000 rows of data. How do I remove the hull ones?

[Image Analyst]

Which axis? Radius in the z direction of the dimple? Radius of the base/substrate in the z direction (if it's on a curved substrate)? Radius in the lateral x/y direction?

[John]

I want to know the total diameter of the circle, which I can obviously find if I have the radius.

I want to know the depth of the dimple itself.

And I want to know the radius of curvature that goes from the outside edge of the dimple down the curve to the center of the dimple. Maybe that's what I need... now I'm not sure actually. :)

[Image Analyst]

OK, so you want both radii of curvatures and the minor (or major) axis of the ellipsoid. If you can get all points on the dimple surface, you can look at the FAQ: http://matlab.wikia.com/wiki/FAQ#How_can_I_fit_an_ellipse_or_other_shape_to_a_set_of_XY_data.3F

[Sean de Wolski]

My snippet above gives you the vertices on the hull, so use setdiff or setxor (with the 'rows' option) to remove them from the original list.

[John]

Sean, I followed your snippet, but instead of creating a random surface like you did to get faces and vertices, I used the faces and vertices I already had. Then I got the delaunayTriangulation and then the convexHull(dt) to the get the hull points. Out of curiousity, I did a scatter3 plot of the hull points, and it is a mess. Looks nothing like it should... almost a straight line with some random dots around it. I'm sure I did it right. I tried with the "onHull" points and with the "k" points, and neither looked correct.

[Sean de Wolski]

really hard to say...

[John]

Here's what I thought...

When I imported the STL and scatter3 plotted it, it wasn't convex, it was concave. The image was upside down, meaning the convexHull was trying to work with a concave image. I multiplied all my values by -1 to change the signs and flipped it over the correct way, but... it still doesn't work.

[John]

And I get this warning every time I do a triangulation. Could this be the cause of my issues?

Warning: Duplicate data points have been detected and removed.
 The Triangulation indices are defined with respect to the unique set of points in delaunayTriangulation property
 X.

[John]

Back after the Thanksgiving holiday...

Any ideas on why this still isn't working?

[Sean de Wolski]

Hi John,

I don't think duplicate values would describe the error you're seeing. It's really difficult for me to come up with an example without a file. But I do have some time on a plane Thursday and Friday so I could try then. Essentially, the goal would be to create an isosurface of a Death Star, probably by taking two partially overlapping spheres and removing the small sphere from both...

Answer 2

[Image Analyst]

I've never used STL files, but if you can get a topographic image (with pixel value being the height values) then you can threshold, and call regionprops(binaryImage, 'EquivDiameter') to get the equivalent circular diameter (ECD).

[Sean de Wolski]

There is a voxelize function on the FEX to convert the output from an STL file to a three-dimensional volumetric image. This might be a quick approach for you, John.

[John]

I'm having trouble getting this function to work. I keep getting errors.

[Image Analyst]

I haven't used it, and don't have plans to. I assume you contacted the author.