How to plot contour on a curved surface for the below problem?

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.

[Nathan Hardenberg]

Where/How do you calculate the surface? I can't seem to find it. You have your two functions A(x) and zc(x), but it is not commented on what they should be.

In the best case you should have a function like

x_new = surface_x_value(y, z)

where you give the current point(s) of your circle in the yz-plane and you get your new x result. Then you

plot3(x_new, y, z)

[Nathan Hardenberg]

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:

1. in "distance_from_pith" you use y^2. But you probably need y.^2
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)

Warning: Imaginary parts of complex X, Y, and/or Z arguments ignored

Warning: Imaginary parts of complex X, Y, and/or Z arguments ignored

    x_new = sqrt(Ri0^2-y.^2);  % calculate points on cylinder (Ri0 = radius of cylinder)
end

[Nathan Hardenberg]

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

Warning: Imaginary parts of complex X, Y, and/or Z arguments ignored

Warning: Imaginary parts of complex X, Y, and/or Z arguments ignored

% 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

Answer 1

[J. Alex Lee]

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