[Scilab-users] A plane intersecting a surface

[Dang Ngoc Chan, Christophe]

Same code but a bit more comments


mode(0)

function [z]=saddle(x, y)
    z = x^2 - y^2
endfunction

A = [1 2 3 4];

function [d]=plane(A, x, y, z)
    d = A(1)*x + A(2)*y + A(3)*z + A(4); // cartesian equation
endfunction

// Dichotomy search of the intercept
function [X]=intercept(A, x1, y1, z1, x2, y2, z2)
    // P1(x1, y1, z1) ; P2(x2, y2, z2) ; [P1P2] is a triangle side
    epsilon = 1e-2;
    d1 = plane(A, x1, y1, z1); // algebraic distance to the plane
    d2 = plane(A, x2, y2, z2);
    if (abs(d1) <= epsilon) then // P1 close enough to the plane
        x = x1;
        y = y1;
        z = z1;
    elseif (abs(d2) <= epsilon) then // P2 close enough to the plane
        x = x2;
        y = y2;
        z = z2;
    else // dichotomy search of a point that is close enough
        cont = %t;
        while cont
            // M(x, y, z) middle point of [P1P2]
            x = 0.5*(x1 + x2); y = 0.5*(y1 + y2); z = 0.5*(z1 + z2);
            d = plane(A, x, y, z);
            if abs(d) <= epsilon then cont = %f // M is close enough
            else
                if (sign(d1*d) == 1) then
                    // M is on the same side of the plane as P1
                    x1 = x; y1 = y; z1 = z; // replace P1 by M
                else
                    x2 = x; y2 = y; z2 = z; // replace P2 by M
                end
            end
        end
    end
    X = [x, y, z];
endfunction

// surface making ... of course in the real life the surface comes from experimental data (no Cartesian equation is attached on))
n = 50;
x = linspace(-2,2,n)';
y = linspace(-1,3,n)';

tic();
// reshape
[X, Y] = ndgrid(x, y);
X = matrix(X, size(X, "*"), 1)';
Y = matrix(Y, size(Y, "*"), 1)';
Z = saddle(X, Y);

// tesselation
tri = delaunay_2(X, Y);
nbtri = size(tri, "r");
points = [];

// detect when the plane crosses the side of a triangle
// then set the intercept point

for i = 1:nbtri
    // extract the points of the tirangle #i
    x1 = X(tri(i, 1)); y1 = Y(tri(i, 1)); z1 = Z(tri(i, 1));
    x2 = X(tri(i, 2)); y2 = Y(tri(i, 2)); z2 = Z(tri(i, 2));
    x3 = X(tri(i, 3)); y3 = Y(tri(i, 3)); z3 = Z(tri(i, 3));
    // compute the algebraic distances
    d1 = plane(A, x1, y1, z1);
    d2 = plane(A, x2, y2, z2);
    d3 = plane(A, x3, y3, z3);
    // check for each pair of points if they are in different half-spaces
    s1 = (sign(d1*d2) == -1);
    s2 = (sign(d2*d3) == -1);
    s3 = (sign(d3*d1) == -1);
    // determine the intercept point if need be
    if s1 then
        points = [points ; intercept(A, x1, y1, z1, x2, y2, z2)];
    end
    if s2 then
        points = [points ; intercept(A, x2, y2, z2, x3, y3, z3)];
    end
    if s3 then
        points = [points ; intercept(A, x3, y3, z3, x1, y1, z1)];
    end
end
duration=toc();

// plot
scf(0);
clf();
plot3d(x, y, feval(x, y, saddle));
// plot3d(X, Y, Z); // doesn't work, don't know why
surf = gce();
param3d(points(:, 1), points(:, 2), points(:, 3));
curve=gce();
curve.mark_mode = "on";
curve.mark_style = 0;
curve.line_mode = "off";
curve.mark_foreground = 5;
curve.mark_size = 1;

disp("t = "+string(duration)+" s.");



--

Christophe Dang Ngoc Chan

Mechanical calculation engineer



This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20180910/73f00401/attachment.htm>