How to project a line on a surface?

Hi,

I am plotting a 3D surface and a straight line together. The line is not on the surface nor intersecting the surface. Now I want to project that line onto the surface.

Anyone can help me how to do this?

Thank you!

[Sargondjani]

What is the shape of your surface? Is it a function? And how do you want to project it?

I assume your surface consists of a finite number of data points, and you use vertical projection. You could use interpolation to let the line closely follow the surface. Assume your surface consists of vectors X,Y,Z

If the line consists of a vectors with x,y coordinates in pairs xp, yp, you could do:

F=scatteredInterpolant(X,Y,Z);

Proj_line = F(xp,yp);

Otherwise you could try to change your data to something like this.

[Sachal Hussain]

It's an anatomical shape, not a function. The line is above the surface so I want to project it straight down on the surface.

The line is plotted in a space so it has X,Y,Z coordinates.

Answer 1

[Star Strider]

Perhaps something like this using the scatteredInterpolant function --

x = linspace(-3, 3);
y = linspace(-4, 4);
[X,Y] = ndgrid(x, y);
zfcn = @(x,y) exp(-(x-mean(x(:))).^2 - (y-mean(y(:))).^2);
Z = zfcn(X,Y);

yline = 0.5*x - 0.5;
zline = 2*ones(size(x));

figure
surf(X, Y, Z, EdgeColor='none', DisplayName='Surface')
hold on
plot3(x, yline, zline, '-g', DisplayName='Offset Line')
hold off
colormap(turbo)
legend(Location='bestoutside')
title('Original')

F = scatteredInterpolant(X(:), Y(:), Z(:));

sline = F(x(:), yline(:));

figure
surf(X, Y, Z, EdgeColor='none', DisplayName='Surface')
hold on
plot3(x, yline, zline, '-g', DisplayName='Offset Line')
plot3(x, yline, sline, '-r', LineWidth=2, DisplayName='Projected Line')
hold off
colormap(turbo)
legend(Location='bestoutside')
title('Original With Projected Offset Line')

●

.

[DGM]

OP's data was a triangulated surface from an STL file, not gridded data. How would you apply this to such a surface?

Answer 2

[KSSV]

Let X, Y, Z be your surface data. And (x,y) be the coordinates of the line data. Use interpolation to get the z values of line from surface and then plot.

z = interp2(X,Y,Z,x,y) ;
surf(X,Y,Z)
shading interp
plot3(x,y,z,'r')

[Sachal Hussain]

But in my case X Y Z are column vectors, how to grid them?

[KSSV]

Let X, Y, Z be your column data.

F=scatteredInterpolant(X,Y,Z) ;
z=F(x, y) ;

[Sachal Hussain]

I tried to follow your suggested way but the result is still not correct. Below is the code I wrote, Please have a look and point out where I am doing wrong. Thak you!

import_stl = stlread('mesh_40_tagli.stl');

figure,hold on,trisurf(import_stl, 'FaceColor','flat', 'LineStyle', 'none');view(3),

% Adding noise to x,y,z coordinates of surface to make them 'unique'

a = import_stl.Points(:,1); t = a == unique(a)'; out1 = sum(cumsum(.0001*t).*t,2)-.0001 + a; % x-coordinates

a = import_stl.Points(:,2); t = a == unique(a)'; out2 = sum(cumsum(.0001*t).*t,2)-.0001 + a; % y-coordinates

a = import_stl.Points(:,3); t = a == unique(a)'; out3 = sum(cumsum(.0001*t).*t,2)-.0001 + a; % z-coordinates

X = [133.4843 131.3938]; Y = [70.1121 66.5661]; Z = [51.3541 27.6482]; % query points

xq = linspace(X(1),X(2),50); yq = linspace(Y(1),Y(2),50); % sampling of query points

F=scatteredInterpolant(out1,out2,out3) ;

zq=F(xq,yq) ;

plot3(xq,yq,zq);

Answer 3

[DGM]

This isn't quite the same, but it's close enough to be refined if the goals required it. We are given a triangulated surface in 3D and any curve in x,y or in x,y,z. Rather than a simple projection atop the model, this finds the intersection between the surface and the projection of the curve through the model.

That means that the projection is not occluded by the model, and the z-position of the curve is completely ignored. If we want only the part of the projection that's not occluded by the model, we'd have to do further work to remove the occluded part.

unzip utahteapot.stl.zip % for the forum

% the triangulated surface
T = stlread('utahteapot.stl');
[F V] = t2fv(T);

% a curve in 3D
nline = 100;
x = linspace(-4,4,nline);
y = 4*sin(1.5*x);
z = linspace(8,10,nline); % this gets ignored anyway

% project the curve through the model along z
[curves Fc Vc] = getprojection(F,V,x,y);

% preview display
plot3(x,y,z,'linewidth',3); hold on; grid on 
for k = 1:numel(curves)
    plot3(curves{k}(:,1),curves{k}(:,2),curves{k}(:,3),'linewidth',3)
end
patch('faces',F,'vertices',V,'facecolor','w','edgecolor','none','facealpha',1);
patch('faces',Fc,'vertices',Vc,'facecolor','r','edgecolor','none','facealpha',0.1);
view(3); view(-32,18); camlight; 
axis equal; grid on
xlabel('X'); ylabel('Y'); zlabel('Z')

function [curves Fc Vc] = getprojection(F,V,x,y)
% [CURVES Fc Vc] = GETPROJECTION(Fmodel,Vmodel,X,Y)
% Given a curve and a triangulated surface in 3D, project the curve 
% through the model along z, and find its intersection with the model surface.
% This is mostly repurposed from the stlslicer() demo (forum question #2175386)
%
% Fmodel, Vmodel describe a 3D surface
% X,Y are vectors describing a curve
%
% CURVES is a cell array of the projected curves in 3D
% Fc, Vc are a surface describing the projection

    % generate the slice surface
    nvc = numel(x);
    zl = [min(V(:,3)) max(V(:,3))] + [-1 1];
    Vc = [x(:) y(:) zl(2)*ones(nvc,1); % vertices
        x(:) y(:) zl(1)*ones(nvc,1)];
    va = (1:nvc).'; % top
    vb = va + nvc; % bottom
    vc = circshift(va,-1); % top
    vd = vc + nvc; % bottom
    Fc = [va vb vc; vd vc vb]; % faces
    Fc([nvc end],:) = []; % get rid of closure
    
    % prepare the surfaces as FV structs for SurfaceIntersection()
    thismodel.vertices = V;
    thismodel.faces = F;
    cutsurf.vertices = Vc;
    cutsurf.faces = Fc;
    
    % find the intersection
    curves = {}; % allocate
    [~,section] = SurfaceIntersection(thismodel,cutsurf);
    
    % SurfaceIntersection() will sometimes return near-duplicate vertices,
    % which breaks curve connectivity. so we need to find duplicate vertices 
    % and repair the edge list to match.
    section = fixdupverts(section);
    
    % the returned list of faces are all degenerate (basically the same as the edges)
    % the only thing this information represents is the boundary line,
    % not a triangulated surface.  we'll need to reorder it to form a curve(s).
    % so reorder the edge segments and identify unique curves
    E = unique(section.edges,'rows'); % only unique edges
    E = E(diff(E,1,2) ~= 0,:); % get rid of degenerate edges
    p = 1;
    while ~isempty(E)
        curves{p} = E(1,:); % pick the first edge
        E = E(2:end,:); % remove it from the pile
        
        % run forward along this curve until we reach an end
        % this is the same as the stlslicer() demo
        atend = false; % we're not at the end of this curve
        while ~atend % run until there are no more connecting segments
            nextv = curves{p}(end,2); % look at the last vertex on this curve
            idx = find(E(:,1) == nextv,1,'first'); % try to find an edge that starts at this vertex
            
            atend = isempty(idx); % we're at the end if there are no more connecting edges
            if ~atend
                curves{p}(end+1,:) = E(idx,:); %#ok<*SAGROW>
                E(idx,:) = [];
            end
        end
        
        % turn around and run backwards from the start of this curve
        % unlike the stlslicer() demo, we're not looking for strictly closed curves
        % so this is how i chose to address the lack of closure
        atend = false; % we're not at the end of this curve
        while ~atend % run until there are no more connecting segments
            nextv = curves{p}(1,1); % look at the first vertex on this curve
            idx = find(E(:,2) == nextv,1,'first'); % try to find an edge that ends at this vertex
            
            atend = isempty(idx); % we're at the end if there are no more connecting edges
            if ~atend
                curves{p} = [E(idx,:); curves{p}]; %#ok<*SAGROW>
                E(idx,:) = [];
            end
        end
        p = p+1;
    end
    
    % construct ordered vertex lists from ordered edge lists
    for p = 1:numel(curves)
        thisE = curves{p};
        thisE = [thisE(:,1); thisE(end,2)];
        curves{p} = section.vertices(thisE,:); %#ok<*AGROW>
    end
    curves = curves.';
end

function section = fixdupverts(section)
    tol = 1E-11; % this should probably be externalized
    [~,IA,IC] = uniquetol(section.vertices,tol,'byrows',true); % not stable sorting!
    
    % try to save some time if there aren't any obvious duplicates
    if numel(IA) == numel(IC)
        return;
    end
    
    % if there are duplicates, fix the edge list
    % could probably be simplified
    dupverts = setdiff(1:numel(IC),IA).';
    ndup = numel(dupverts);
    for k = 1:ndup
        refvert = find(IC == IC(dupverts(k)),1,'first');
        section.edges(section.edges == dupverts(k)) = refvert;
    end
end

This uses this file from the FEX (attached):

https://www.mathworks.com/matlabcentral/fileexchange/48613-surface-intersection

See also:

https://www.mathworks.com/matlabcentral/answers/2175386-surface-area-at-given-hight-of-an-stl-file