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

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')

[Jonas]

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]

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.

Answer 1

[Adam Danz]

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

[Kylen]

Thank you, Mr. Danz. This was extremely helpful.

If I could trouble you with one more issue. This script is being used to prepare data for a Simscape model. I am using a Grid Surface Block for contact modeling and want to input the shape you helped me solve above.

However, the input requirements to use the Grid Surface don't permit x,y,z scatter points. Rather, it follows the "ndgrid" format: "The grid surface has a similar format to the grid generated by the ndgrid function in MATLAB®. To model a grid surface, use the X Grid Vector and Y Grid Vector parameters to specify the coordinates in the x and y directions of the grid surface, then use the Z Heights parameter to specify the elevations of the grid points in the z direction."

Do you have any insight on how to adjust the x,y,z scatter points that you were able to generate to be useful using a Simscape Grid Surface?

More documentation here: https://www.mathworks.com/help/sm/ref/gridsurface.html

Thank you again.

[Adam Danz]

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]

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.

Answer 2

[Adam Danz]

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

[Kylen]

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!