[Scilab-users] TR: A plane intersecting a surface

[Dang Ngoc Chan, Christophe]

Hello,



[sorry Paul, it was mistakedly sent only to you]



A far from optimal solution - not vectorised, the points are searched two times (almost each triangle side is common to another one -

using the Delaunay tessellation.



Requires the Atoms module cglab


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)
    epsilon = 1e-2;
    x = x1; y = y1; z = z1;
    d1 = plane(A, x1, y1, z1);
    d2 = plane(A, x2, y2, z2);
    if (abs(d1) <= epsilon) then
        x = x1;
        y = y1;
        z = z1;
    elseif (abs(d2) <= epsilon) then
        x = x2;
        y = y2;
        z = z2;
    else
        cont = %t;
        while cont
            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
            else
                if (sign(d1*d) == 1) then
                    x1 = x; y1 = y; z1 = z;
                else
                    x2 = x; y2 = y; z2 = z;
                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 = 10;
n = 50;
x = linspace(-2,2,n)';
y = linspace(-1,3,n)';

// 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 as the middle

for i = 1:nbtri
    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));
    d1 = plane(A, x1, y1, z1);
    d2 = plane(A, x2, y2, z2);
    d3 = plane(A, x3, y3, z3);
    s1 = (sign(d1*d2) == -1);
    s2 = (sign(d2*d3) == -1);
    s3 = (sign(d3*d1) == -1);
//    bool(i) = s1|s2|s3;
    if s1 then
        //points = [points ; 0.5*([x1, y1, z1] + [x2, y2, z2])];
        points = [points ; intercept(A, x1, y1, z1, x2, y2, z2)];
    end
    if s2 then
        //points = [points ; 0.5*([x2, y2, z2] + [x3, y3, z3])];
        points = [points ; intercept(A, x2, y2, z2, x3, y3, z3)];
    end
    if s3 then
        //points = [points ; 0.5*([x3, y3, z3] + [x1, y1, z1])];
        points = [points ; intercept(A, x3, y3, z3, x1, y1, z1)];
    end
end

// plot
scf(0);
clf();
plot3d(x, y, feval(x, y, saddle));
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;



--

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/ee715c50/attachment.htm>