How to plot contour on a curved surface for the below problem?
4 次查看(过去 30 天)
显示 更早的评论
Hello,
Can anyone here help me with the plotting contours on a curved surfaceThe code should be able to generate the plot as in ooutput image.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1478986/image.png)
3 个评论
Nathan Hardenberg
2023-9-11
Yeah, works as expected! But your "surface_x_value" function just does not give correct results. Since I don't know how you calculate those I can't help you further. I see that you do an optimization, but somewhere lies an error.
Two small things I found:
- in "distance_from_pith" you use y^2. But you probably need y.^2
- Your inital guess seems to be very weird. Your distance to pith should not be a good guess for your x-value. Maybe a better guess would be the surface of the cylinder.
As you can see below, if I draw the red circles on the cylinder it works. (except for the imaginary parts, since they are not on the cylinder)
... % Only changed the surface_x_value function (see below)
function x_new = surface_x_value(y, z, k, zc0, Ri0)
x_new = sqrt(Ri0^2-y.^2); % calculate points on cylinder (Ri0 = radius of cylinder)
end
Nathan Hardenberg
2023-9-14
Here is all the code that is needed to function. Again only your code except the function at the end.
clear all
close all
z_root=0;
z_top=150;
dk_log=40; % extension of drawn knot, i.e. length in addition to Ri0
Ri0=80;
% Make a drawing to illustrate a part of a log (two circles in 3D space representing the circumference of the log
% and dashed lines representing the pith and a plane through the log.)
figure(1)
hold on
teta_mesh=[0:5:360]/180*2*pi;
plot3([-Ri0 Ri0 Ri0 -Ri0 -Ri0]',[0 0 0 0 0]',[z_root z_root z_top z_top z_root],':k')
plot3([Ri0*cos(teta_mesh)]',[Ri0*sin(teta_mesh)]',[teta_mesh*0+z_root]','-k')
plot3([Ri0*cos(teta_mesh)]',[Ri0*sin(teta_mesh)]',[teta_mesh*0+z_top]','-k')
plot3([0 0]',[0 0]',[z_root z_top],':r')
box on
axis equal
view(3)
% Define parameters for the geometry of a knot according to Foley's
% model/thesis.
zc0=40; % z-position of the knot apex at the pith
Ax100=41.5; %page 88
L100=28.6; %page 88
Dl=1.22e-2*L100+0.1252; %page 91
Cax=9.7e-3*Ax100+0.1725; %page 91
kappa=0.95; %page 93
c1=-1.458e-3; %page 91
c2=0.01*Ax100+0.1458; %page 92,
x_mesh=[0:4:(Ri0+dk_log)];
lx_mesh=-1e-7*x_mesh.^4+3e-5*x_mesh.^3-4e-3*x_mesh.^2+Dl*x_mesh; %page 91
bx_mesh=lx_mesh*kappa; %page 93
zc_mesh=c1*x_mesh.^2+Cax*x_mesh; %page 92
knot_rad_long_xpdk=(-1e-7*(Ri0+dk_log)^4+3e-5*(Ri0+dk_log)^3-4e-3*(Ri0+dk_log)^2+Dl*Ri0)/2;
k=1/kappa; %page 94, eq. 5.6
% Add to the drawing lines representing the geometry of a knot
XYZp=zeros(181,3);
XYZp(:,1)=Ri0+dk_log;
zc=c1*(Ri0+dk_log)^2+Cax*(Ri0+dk_log);
plot3(x_mesh,x_mesh*0,zc_mesh+zc0,'--k')
plot3(x_mesh,x_mesh*0,zc_mesh+zc0-lx_mesh/2,'-k')
plot3(x_mesh,x_mesh*0,zc_mesh+zc0+lx_mesh/2,'-k')
knot_ovals_x_pos=[30:30:120];
for xov=knot_ovals_x_pos
knot_rad_long_xpdk=(-1e-7*xov^4+3e-5*xov^3-4e-3*xov^2+Dl*xov)/2;
dp=knot_rad_long_xpdk;
XYZp=zeros(181,3);
XYZp(:,1)=xov;
pos=0;
zc=c1*xov^2+Cax*xov;
for v=0:2:360
pos=pos+1;
vr=v/180*pi;
ypos=dp*cos(vr)/k;
zpos=dp*sin(vr);
XYZp(pos,2:3)=[ypos zpos];
if pos>1
plot3(XYZp((pos-1):pos,1),XYZp((pos-1):pos,2),XYZp((pos-1):pos,3)+zc0+zc,':k');
end
end
end
% Add to the figure a number of contour curves (red color) in the 'bulged part' of a 'growth surface' with
% Ri0 = 80 mm, see pages 62-67 in Foley's thesis.
Bbump=2.0; %page 93
Abump=(0.0217*L100-0.2); %page 93
Aexp=0.0056*L100+1.96; %page 93
kappa_limit=0.99; %page 94
Rik=80;
%--- Add the missing code here ---
% Define the function A(x) using equation (5.3.a)
A = @(x) -1e-7*x.^4 + 3e-5*x.^3 - 4e-3*x.^2 + (1.22e-2*L100 + 0.1252).*x;
% Define the function zc(x) using equation (5.3.b)
zc = @(x) -1.458e-3*x.^2 + (9.7e-3*Ax100 + 0.1725).*x;
% Plotting the growth surface contours
p_vals = [15:5:75];
for zp = p_vals % page 65, fig 4.7
theta_vals = linspace(0, 2 * pi, 100);
% Calculate the distance from the pith to the contour curve
distance_to_pith = sqrt(Ri0^2 - zp^2);
% Adjust the radius based on the distance to the pith and aspect ratio k
r_vals = distance_to_pith * k; % k is the aspect ratio
y_vals = r_vals .* cos(theta_vals); % Adjusted for YZ plane
z_vals = 35 + zc0 + r_vals .* sin(theta_vals); % Adjusted for YZ plane
% Get the x values for the growth surface
x_new = arrayfun(@(y, z) surface_x_value(y, z, k, zc0, Ri0), y_vals, z_vals);
plot3(x_new, y_vals, z_vals, 'r'); % plotting the growth surface
end
% Adding labels and title
xlabel('X-axis')
ylabel('Y-axis')
zlabel('Z-axis')
% Function to calculate the x-coordinate of the growth surface for a given y and z
%-----------------------------------
grid on
axis equal
function x_new = surface_x_value(y, z, k, zc0, Ri0)
x_new = sqrt(Ri0^2-y.^2);
end
回答(1 个)
J. Alex Lee
2023-9-16
This is how I would do it with the bump.
% parameters
z_root=0;
z_top=150;
Ri0 = 80;
trunk.TRes = 360/5;
trunk.ZRes = 30;
dk_log=40; % extension of drawn knot, i.e. length in addition to Ri0
zc0=40; % z-position of the knot apex at the pith
Ax100=41.5; %page 88
L100=28.6; %page 88
Dl=1.22e-2*L100+0.1252; %page 91
Cax=9.7e-3*Ax100+0.1725; %page 91
c1=-1.458e-3; %page 91
c2=0.01*Ax100+0.1458; %page 92,
kappa = .7; % free choice how oval are the rings
% 1. the symmetrical gaussian bump
bumpH = 10; % free choice approx height of bump relative to Ri0
bumpVar = 200; % free choice how wide/gradual bump is
numRings = 15; % free choice number of rings to draw
ringRes = 200;
% knot geometry
knotCenterFcn = @(r)(polyval([c1,Cax,zc0],r));
knotRadiusFcn = @(r)(polyval([-1e-7,3e-5,-4e-3,Dl,0],r)/2);
knot.centerline.n = 50;
knot.centerline.r = linspace(0,Ri0 + dk_log,knot.centerline.n).';
knot.centerline.theta = 0 * ones(knot.centerline.n,1);
[knot.centerline.x,knot.centerline.y] = pol2cart(knot.centerline.theta,knot.centerline.r);
knot.centerline.z = knotCenterFcn(knot.centerline.r);
% create gridded surface coordinates (easy way)
[trunk.X,trunk.Y,trunk.Z] = cylinder(Ri0*ones(trunk.ZRes,1),trunk.TRes);
trunk.Z = trunk.Z*(z_top-z_root) + z_root;
% warning: this will break for knot.centerline.theta .ne. 0
[knot.surf.Z,knot.surf.Y,knot.surf.X] = cylinder(knotRadiusFcn(knot.centerline.r));
knot.surf.Y = knot.surf.Y * kappa;
%--- Add the missing code here ---
% 1. model the bump as a gaussian
% 2. generate the contours using built in contourc
% 3. use Adam Danz's contour line extractor
% 4. wrap the contour line cooridnates on the trunk cylinder with trig
%
% steps 2 and 3 can be done manually, but why not use built-in
zPith = knotCenterFcn(Ri0 + bumpH);
rPithAtTrunk = knotRadiusFcn(Ri0 + bumpH);
% need to calculate the gaussian exponential factor A s.t. at the height of
% the bump for rPothAtTrunk = bumpH
% bumpH = A*exp(... + kappa*y^2/A)
A = bumpH/exp(-(rPithAtTrunk^2)/bumpVar);
% the double gaussian function
bumpFcn = @(z,y)(A*exp(-((z-zPith).^2+(y/kappa).^2)/bumpVar) + Ri0);
z = linspace(z_root,z_top,ringRes);
y = linspace(-Ri0,Ri0,ringRes).';
XFlat = bumpFcn(z,y);
% some logic to space the rings
skewedlevels = linspace(0,1,numRings).^.3;
levelLims = log10([1e-8,bumpH]);
loglevels = skewedlevels*diff(levelLims) + levelLims(1);
levels = Ri0 + 10.^loglevels;
% get the contour coordinates using built-ins plus FEX
[cm] = contourc(z,y,XFlat,levels);
[ct,~] = getContourLineCoordinates(cm);
grps = unique(ct.Group);
for i = numel(grps):-1:1
grp = ct(ct.Group==grps(i),:);
% contour assumed z-axis is height map
% so need to switch some coordinates up
x = grp.Level;
y = grp.Y;
% use trig to wrap around trunk
rings(i).x = x .* cos(y./x);
rings(i).y = x .* sin(y./x);
rings(i).z = grp.X;
end
% draw
figure(1);
tl = tiledlayout(2,2);
ax(1) = nexttile([1,2]);
ax(2) = nexttile(tl);
ax(3) = nexttile(tl);
set(ax,"DataAspectRatio",[1 1 1],"NextPlot","add","Box","on")
mesh(ax(1),trunk.X,trunk.Y,trunk.Z,"EdgeAlpha",0.1,"FaceAlpha",0,"EdgeColor","k")
mesh(ax(1),knot.surf.X+knot.centerline.x,knot.surf.Y+knot.centerline.y,knot.surf.Z+knot.centerline.z,"EdgeAlpha",0.1,"FaceAlpha",0,"EdgeColor","k")
plot3(ax(1),knot.centerline.x,knot.centerline.y,knot.centerline.z,"--b")
for i = 1:numel(grps)
plot3(ax(1),rings(i).x,rings(i).y,rings(i).z,"-","Color",[0.8,0,0])
end
copyobj(allchild(ax(1)),ax(2))
copyobj(allchild(ax(1)),ax(3))
view(ax(1),[ 39.0438 17.8241])
view(ax(2),[90,0])
view(ax(3),[0,0])
xlabel(ax,"X-axis")
ylabel(ax,"Y-axis")
zlabel(ax,"Z-axis")
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Lighting, Transparency, and Shading 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!