how to project a contour along a surface plane?

How to plot the projection of a function f(x, y, z) onto a plane with normal vector N, passing through P0?

The goal is that the contour lines appear on top of the inclined plane. See code and plot below.

I guess I should use hgtransform on it, but I am getting nowhere.

N = 2:4;
P0 = 5:7;
d = -N*P0';

x = -5:5;
y = x;
[x, y] = meshgrid(x, y);
z = -(N(1).*x + N(2).*y + d)./N(3);

f = sqrt(x.^2 + y.^2 + z.^2)*0.1;

figure;
hold on;
surface(x, y, z);
xlabel("x");
ylabel("y");
zlabel("z");
shading interp;
cb = colorbar;
cb.Label.String = "z";
cb.Limits = [min(z, [], "all") max(z, [], "all")];

M = eye(4);
HGT = hgtransform();
HGT.Matrix = M;

contour(x, y, f, "k", "ShowText", "on", "Parent",HGT);

grid on;
box on;
view(120, 60)

Answer 1

[Geraldo Rebouças]

The challenge was more math than MATLAB itself.

Here is the solution:

close all;

% defining the plane:
N = 2:4; % normal vector to the plane
P0 = 5:7; % point on the plane
d = -N*P0';

x = -5:5;
y = x;
[x, y] = meshgrid(x, y);
z = -(N(1).*x + N(2).*y + d)./N(3);

% Coordinate transformation matrix from the plane z = 0, to
% the plane N(1)*x + N(2)*y +N(3)*z + d = 0.
M = eye(4);
M(3, 1) = -N(1)/N(3);
M(3, 2) = -N(2)/N(3);
M(3, 4) = -d/N(3);
M(3, 4) = M(3, 4)*1.01; % to ease viewing the contour lines and text

% function to be projected on the plane:
f = sqrt(x.^2 + y.^2 + z.^2)*0.1;

figure;
hold on;
surface(x, y, z);
xlabel("x");
ylabel("y");
zlabel("z");
shading interp;
cb = colorbar;
cb.Label.String = "z";
cb.Limits = [min(z, [], "all") max(z, [], "all")];

HGT = hgtransform();
HGT.Matrix = M;

Warning: The new value for the Matrix property may cause rendering problems.

hold on;
contour(x, y, f, "k", "ShowText", "on", "Parent", HGT);
grid on;
box on;
view(120, 60);

I don't know what's the problem with this transformation (i.e. the warning message).

Note: the above only works for projections over a plane.

Answer 2

[Matt J]

This looks applicable:

https://www.mathworks.com/matlabcentral/fileexchange/8589-contourz?s_tid=srchtitle

[Geraldo Rebouças]

Thanks for the suggestion! That is very interesting!

I did not understand what is the point of using the mouse, so I just commented it below.

I also had to update the input and output to the view() function.

I will stick to my solution because it is more concise, and do not rely on old code nor undocumented functions, e.g. contours().

%% contourz example:
figure
[x,y,v] = peaks;
z=-(x.^2+y.^2);
surf(x,y,z,'facecolor','none','edgealpha',.1)
hold on
contourz(x,y,z,v, 5);
view(2)
contourz('clabel');
view(3);

%% Question:
N = 2:4; % normal vector to the plane
P0 = 5:7; % point on the plane
d = -N*P0';

x = -5:5;
y = x;
[x, y] = meshgrid(x, y);
z = -(N(1).*x + N(2).*y + d)./N(3);

% function to be projected on the plane
f = sqrt(x.^2 + y.^2 + z.^2)*0.1;

figure;
hold on;
surface(x, y, z);
xlabel("x");
ylabel("y");
zlabel("z");

contourz(x,y,z,f,10);
view(120, 60);

%%
function handles = contourz(xi,yi,zi,vi,vv,theColor)
%CONTOURZ   3D contour plot on a surface
%   Similar to contour3 but contours are drawn over any surface
%   This function is also used to create clabels.
%
%   Syntax:
%      H = CONTOURZ(X,Y,Z,V,VV,COLOR)
%      H = CONTOURZ('clabel')
%
%   Inputs:
%      XI, YI   2-D arrays
%      ZI       Surface where contours will be drawn, defined at XI,YI
%      VI       2-D array to contour, defined at XI,YI
%      VV       Contour levels or number of contours [ 10 ]
%      COLOR    Color of contours, otherwise a patch is created
%      'clabel' Clabels are created manually as by clabel(cs,'manual')
%               Clabels are done clicking with left mouse button
%               until click with other button
%
%   Outputs:
%      H   handles of lines, patches for contours or handles for labels
%          (marker and text) for clabels
%
%   Comments:
%      This may be useful if you wanna add contours to a
%      surface created by surf(x,y,z,v).
%      Notice that if v = z, is the same as use contour3
%
%   Example:
%      figure
%      [x,y,v] = peaks;
%      z=-(x.^2+y.^2);
%      surf(x,y,z,'facecolor','none','edgealpha',.1)
%      hold on
%      contourz(x,y,z,v);
%      view(2)
%      contourz('clabel');
%      view(3)
%
%   MMA 8-2004, [e-mail]
%   Department of physics
%   University of Aveiro
handles = [];
% --------------------------------------------------------------------
% clabel
% --------------------------------------------------------------------
if isequal(xi,'clabel')
  is_hold = ishold;
  hold on
  % vw = view;
  [az, el] = view;
  %view(2);
  h=[];
  % while 1
  %   [x,y,mouse] = ginput(1);
  %   if mouse ~= 1
  %     break
  %   end
  %   tag = get(gco,'tag');
  %   if isequal(tag,'contourz')
  %     % get closest point:
  %     xdata = get(gco,'xdata');
  %     ydata = get(gco,'ydata');
  %     zdata = get(gco,'zdata');
  %     dist = (xdata-x).^2 + (ydata-y).^2;
  %     i=find(dist == min(dist));
  %     level = get(gco,'userdata');  
  %     p=plot3(xdata(i),ydata(i),zdata(i),'r+');
  %     t=text(xdata(i),ydata(i),zdata(i),num2str(level), ...
  %       'HorizontalAlignment','left');
  %     h=[h;p;t];
  %   end
  % end
  % view(vw)
  view(az, el)
  handles = h;
  if ~is_hold
    hold off
  end
  return
end
% --------------------------------------------------------------------
% contour
% --------------------------------------------------------------------
if nargin < 4
  disp('» contourz: arguments required');
  return
end
if nargin == 4
  vv    = 10;
end
if nargin == 5 & isstr(vv)
  theColor=vv;
  vv = 10;
end
do_line = 1;
eval('lineColor = theColor;','lineColor = [];')
if isempty(lineColor)
  do_line  = 0;
else
  if ~isempty(str2num(lineColor)); % allow colors as ['r g b'] (as string)
    lineColor = str2num(lineColor);
  end
end
c_xz = contours(xi,zi,vi,vv);
c_yz = contours(yi,zi,vi,vv);
i = 1;
n = 1;
if isempty(c_xz) | isempty(c_yz)
  return
end
while 1
  lev(n) = c_xz(1,i);
  nvals  = c_xz(2,i);
  x{n}   = c_xz(1,i+1:i+nvals);
  y{n}   = c_yz(1,i+1:i+nvals);
  z{n}   = c_xz(2,i+1:i+nvals);
  i = i+nvals+1;
  if i > size(c_xz,2)
    break
  else
    n=n+1;
  end
end
h=[];
for n=1:length(x)
  if do_line
    p=plot3(x{n},y{n},z{n},'color',lineColor,'userdata',lev(n),'tag','contourz'); hold on
  else
    xdata = [x{n} nan];
    ydata = [y{n} nan];
    zdata = [z{n} nan];
    level = lev(n);
    vdata = [repmat(level,size(x{n})) nan];
    p = patch('XData',xdata,'YData',ydata,'ZData',zdata, ...
              'CData',vdata,'userdata',level,'tag','contourz', ...
              'facecolor','none','edgecolor','flat');
  end
  h=[h;p];
end
handles = h;
end