Clear Filters
Clear Filters

How to get "Clean" edges on a surface plot?

35 views (last 30 days)
Hello All,
I am trying to plot a surface using the surf() command with some imported shape data from a CSV. The shape in question is created by revolving our data from 0 to 180 degrees. The data starts out as axial position and a corresponding radius. Essentially, it's a half-pipe with asymmetrically tapered ends.
I've been able to import my data, revolve it, and use meshgrid() and griddata() to create the Z coordinate matrix. However, as you can see from my attached screenshot, there are "flats" at the ends of the revolved shape. I have a feeling this has something to do with the interpolation, but I'm not very familiar with this type of stuff. I can't seem to get rid of the flat protrusions.
In my attempts to rectify this, the closest I could get was ending up with jagged edges instead. But this won't work either since this shape data is going to be imported into a Simscape model, and we want it to be as true to life as possible (no flat protrusions and no jagged edges).
Any ideas? Data attached as a text file (wouldn't let me attach a CSV) and code pasted below.
Thanks Everyone!
GS = readtable('GSDataEdited.txt');
GS = table2array(GS);
GSx = GS(:, 1);
GSrad = GS(:, 2);
% Define the angle step (in degrees)
angle_step = 7.5;
all_points = [];
% Loop through angles from 0 to 180 degrees
for angle = 0:angle_step:180
% Convert angle to radians
theta = deg2rad(angle);
% Create rotation matrix for current angle
R = [1 0 0; 0 cos(theta) -sin(theta); 0 sin(theta) cos(theta)];
% Rotate the coordinates
rotated_coords = R * [GSx'; GSrad'; zeros(size(GSx'))];
% Append the rotated coordinates to the matrix
all_points = [all_points, rotated_coords];
end
all_points = all_points';
GSx = all_points(:, 1);
GSy = all_points(:, 2);
GSz = all_points(:, 3);
GSall = [GSx,GSy,GSz;GSx,GSy,-GSz];
GSall = unique(GSall, 'rows');
% Define grid for interpolation
[GSX,GSY] = meshgrid(min(GSx):1:max(GSx), min(GSy):1:max(GSy)); % Define the grid spacing
% Interpolate data onto grid
ZVect = griddata(GSx, GSy, GSz, GSX, GSY, 'linear');
Warning: Duplicate data points have been detected and removed - corresponding values have been averaged.
XVect = GSX(1, :);
YVect = GSY(:, 1);
ZVect(isnan(ZVect))=0;
surf(XVect, YVect, ZVect, 'EdgeColor', 'none')
  2 Comments
Jonas
Jonas on 5 Aug 2024
i am not exactly sure what you want not to appear, amybe one of the two versions below is what you are searching for
GS = readtable('GSDataEdited.txt');
GS = table2array(GS);
GSx = GS(:, 1);
GSrad = GS(:, 2);
% Define the angle step (in degrees)
angle_step = 7.5;
all_points = [];
% Loop through angles from 0 to 180 degrees
for angle = 0:angle_step:180
% Convert angle to radians
theta = deg2rad(angle);
% Create rotation matrix for current angle
R = [1 0 0; 0 cos(theta) -sin(theta); 0 sin(theta) cos(theta)];
% Rotate the coordinates
rotated_coords = R * [GSx'; GSrad'; zeros(size(GSx'))];
% Append the rotated coordinates to the matrix
all_points = [all_points, rotated_coords];
end
all_points = all_points';
GSx = all_points(:, 1);
GSy = all_points(:, 2);
GSz = all_points(:, 3);
GSall = [GSx,GSy,GSz;GSx,GSy,-GSz];
GSall = unique(GSall, 'rows');
% Define grid for interpolation
[GSX,GSY] = meshgrid(min(GSx):1:max(GSx), min(GSy):1:max(GSy)); % Define the grid spacing
% Interpolate data onto grid
ZVect = griddata(GSx, GSy, GSz, GSX, GSY, 'linear');
Warning: Duplicate data points have been detected and removed - corresponding values have been averaged.
XVect = GSX(1, :);
YVect = GSY(:, 1);
ZVect(isnan(ZVect))=0;
% a more technical view, such that you can see the mesh only
figure;
surf(XVect, YVect, ZVect, 'FaceColor','none');
% a more smooth experience
figure;
surf(XVect, YVect, ZVect, 'EdgeColor', 'none','FaceColor','interp')
Kylen
Kylen on 5 Aug 2024
Hi and thank you for your answer.
Perhaps the attached images will make my intent more clear. The areas highlighted in red are what I'm trying to omit from the surface plot. I have also included a screenshot of the revolved profile from Solidworks.
Thanks again.

Sign in to comment.

Accepted Answer

Adam Danz
Adam Danz on 6 Aug 2024
Edited: Adam Danz on 6 Aug 2024
Conversion from 3D scatter points to a surface is not straightforward. Here, I'm using Delaunay triangulation to create a mesh of triangles that approximates the surface defined by the scatter points.
Code from OP
GS = [0,0
0.225,2.35
1.325,5.525
2.405,7.3375
3.805,9.15
5.76,11.3125
7.45,12.95
9.35,14.5
11.855,16.0875
14.375,17.325
16.035,18.025
18.495,18.9125
21.22,19.6875
24.095,20.2375
27.14,20.5375
29.575,20.6125
32.495,20.625
32.5,20.625
101,20.625
106,20.625
117.225,20.4125
125.15,20.025
133.785,19.425
138.675,18.9875
145.53,18.1375
152.38,16.7
157.47,15.1
160.775,13.8375
164.645,12.1875
169.765,9.8125
174.14,7.6
177.11,5.9125
178.795,4.8375
180.225,3.85
181.895,2.35];
% GS = table2array(GS);
GSx = GS(:, 1);
GSrad = GS(:, 2);
% Define the angle step (in degrees)
angle_step = 7.5;
all_points = [];
% Loop through angles from 0 to 180 degrees
for angle = 0:angle_step:180
% Convert angle to radians
theta = deg2rad(angle);
% Create rotation matrix for current angle
R = [1 0 0; 0 cos(theta) -sin(theta); 0 sin(theta) cos(theta)];
% Rotate the coordinates
rotated_coords = R * [GSx'; GSrad'; zeros(size(GSx'))];
% Append the rotated coordinates to the matrix
all_points = [all_points, rotated_coords];
end
Show the scatter points
plot3(all_points(1,:), all_points(2,:), all_points(3,:),'k*')
axis equal
Triangulate and plot the surface
x = all_points(1,:);
y = all_points(2,:);
z = all_points(3,:);
figure()
tri = delaunay(x,y);
Warning: Duplicate data points have been detected and removed.
Some point indices will not be referenced by the triangulation.
h = trisurf(tri, x, y, z, 'EdgeColor','none','FaceColor','interp');
axis equal
You can play around with it to get different coloration,
figure
h = trisurf(tri, x, y, z,'EdgeColor','none','FaceColor',[.5 .5 .5]);
axis equal
camlight
material dull
  3 Comments
Adam Danz
Adam Danz on 6 Aug 2024
Now I see why you were using griddata.
I don't have experience with simscape's modeling graphics. Continuing with your griddata work, the parts of the grid you want to remove are on a single plane at the bottom. If simscape's modeling graphics accepts NaN values, you could fill in those corner values with NaNs. Some MATLAB graphics functions ignore NaN values. To identify those values, I'd start by using inpolygon which determines what coordinates are within a polygon. The polygon boundary would be defined by your GS coordinates.
Kylen
Kylen on 6 Aug 2024
Unfortunately, the NaN's can't be ignored. I can replace them with zeros easily enough, but I end up with the "flat corners" like I have in my original post.
Thanks again, your help was still very much appreciated.

Sign in to comment.

More Answers (1)

Adam Danz
Adam Danz on 5 Aug 2024
Use convhull to enclose the polyton and trisurf to plot the results.
Code from OP
GS = [0,0
0.225,2.35
1.325,5.525
2.405,7.3375
3.805,9.15
5.76,11.3125
7.45,12.95
9.35,14.5
11.855,16.0875
14.375,17.325
16.035,18.025
18.495,18.9125
21.22,19.6875
24.095,20.2375
27.14,20.5375
29.575,20.6125
32.495,20.625
32.5,20.625
101,20.625
106,20.625
117.225,20.4125
125.15,20.025
133.785,19.425
138.675,18.9875
145.53,18.1375
152.38,16.7
157.47,15.1
160.775,13.8375
164.645,12.1875
169.765,9.8125
174.14,7.6
177.11,5.9125
178.795,4.8375
180.225,3.85
181.895,2.35];
% GS = table2array(GS);
GSx = GS(:, 1);
GSrad = GS(:, 2);
% Define the angle step (in degrees)
angle_step = 7.5;
all_points = [];
% Loop through angles from 0 to 180 degrees
for angle = 0:angle_step:180
% Convert angle to radians
theta = deg2rad(angle);
% Create rotation matrix for current angle
R = [1 0 0; 0 cos(theta) -sin(theta); 0 sin(theta) cos(theta)];
% Rotate the coordinates
rotated_coords = R * [GSx'; GSrad'; zeros(size(GSx'))];
% Append the rotated coordinates to the matrix
all_points = [all_points, rotated_coords];
end
all_points = all_points';
GSx = all_points(:, 1);
GSy = all_points(:, 2);
GSz = all_points(:, 3);
GSall = [GSx,GSy,GSz;GSx,GSy,-GSz];
GSall = unique(GSall, 'rows');
% Define grid for interpolation
[GSX,GSY] = meshgrid(min(GSx):1:max(GSx), min(GSy):1:max(GSy)); % Define the grid spacing
% Interpolate data onto grid
ZVect = griddata(GSx, GSy, GSz, GSX, GSY, 'linear');
Warning: Duplicate data points have been detected and removed - corresponding values have been averaged.
XVect = GSX(1, :);
YVect = GSY(:, 1);
Zorig = ZVect;
ZVect(isnan(ZVect))=0;
Perform convex hull and plot results
k = convhull(GSX,GSY,ZVect,'Simplify',true);
trisurf(k,GSX,GSY,ZVect,'EdgeColor','none','FaceColor', 'interp')
axis padded
  1 Comment
Kylen
Kylen on 5 Aug 2024
Hi and thank you for the answer.
Unfotunately, I don't think I was descriptive enough in my original post. I'm not try to include the squared off edges, in fact I am trying to omit them. Please see the attached images. One is a Solidworks profile revololution that shows what the "real" geometry looks like. The other highlights in red the area that I'm trying to omit.
Thanks again for your response and apologies for not being clear enough in my original post!

Sign in to comment.

Products


Release

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!