[Scilab-users] A plane intersecting a surface

classic Classic list List threaded Threaded
10 messages Options
Carrico, Paul Carrico, Paul
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] A plane intersecting a surface

Dear All

I’ve not been using scilab for a while, but I’ve a good opportunity to dive into it once again ;-)

Is there a tool implemented into Scilab to determine the cross section of a 3D (experimental) surface and a plane?

Note that :

-          the curve has not a Cartesian equation and it is composed of a point cloud coming for experimental measurement,

-          ideally the tool looks for the closest out of plane points in order to perform interpolations

Before reinventing the wheel, I’m wondering if something exists.

Nb: I built a saddle surface, but of course only points (not necessary equally spaced) exist in the real life.

Thanks for any advice and suggestion

Paul

function [z]=saddle(x, y)

    z = x^2 - y^2

endfunction

 

// surface making ... of course in the real life the surface comes from exprimental data (no cartesian equation is attached on))

n = 50;

x = linspace(-2,2,n)';

y = linspace(-1,3,n)';

z = feval(x,y,saddle);

plot3d(x,y,z);

 

// plane equation: ax + by + cz + d = 0

 

 

 

EXPORT CONTROL :
Cet email ne contient pas de données techniques
This email does not contain technical data

 


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: A plane intersecting a surface

Hi Paul,

 

I am not aware of such tool.

 

To extract the points in the experimental cloud that are within a given distance from the plane, use equation in:

https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_plane

 

If your cloud of points behaves well enough, you can interpolate it first into a dense grid (for instance using cshep2d) and then extract only the points that are very close to the plane.

 

What would you like to do next with those points on the 2D plane?

If you just need to highlight them in 3D view, then scatter3 is your buddy.

 

Regards,

Rafael

 


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: A plane intersecting a surface

In reply to this post by Carrico, Paul
Hello Paul

Do you have an example of typical 3D data ? Can it be safely projected on a principal plane ? Is it noisy ? 

S.

Le 6 sept. 2018 à 21:50, Carrico, Paul <[hidden email]> a écrit :

Dear All

I’ve not been using scilab for a while, but I’ve a good opportunity to dive into it once again ;-)

Is there a tool implemented into Scilab to determine the cross section of a 3D (experimental) surface and a plane?

Note that :

-          the curve has not a Cartesian equation and it is composed of a point cloud coming for experimental measurement,

-          ideally the tool looks for the closest out of plane points in order to perform interpolations

Before reinventing the wheel, I’m wondering if something exists.

Nb: I built a saddle surface, but of course only points (not necessary equally spaced) exist in the real life.

Thanks for any advice and suggestion

Paul

function [z]=saddle(x, y)

    z = x^2 - y^2

endfunction

 

// surface making ... of course in the real life the surface comes from exprimental data (no cartesian equation is attached on))

n = 50;

x = linspace(-2,2,n)';

y = linspace(-1,3,n)';

z = feval(x,y,saddle);

plot3d(x,y,z);

 

// plane equation: ax + by + cz + d = 0

 

 

 

EXPORT CONTROL :
Cet email ne contient pas de données techniques
This email does not contain technical data

 


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Christophe Dang Ngoc Chan Christophe Dang Ngoc Chan
Reply | Threaded
Open this post in threaded view
|

Re: A plane intersecting a surface

In reply to this post by Carrico, Paul
Hello,

> De : users [mailto:[hidden email]] De la part de Rafael Guerra
> Envoyé : samedi 8 septembre 2018 14:52
>
> If your cloud of points behaves well enough, you can interpolate it first into a dense

If nobody is expert in this field, then I could invoke a memory when I was a student.
I've heard about an algorithm using intercept with tetrahedrons,
it was used for surface rendering.

So you might perform a Delaunay tessellation of your cloud,
determine which tetrahedrons are cut
and determine the coordinates of the intercepts.

Or ask some CGI  specialists.

HTH

Regards


--
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.
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Carrico, Paul Carrico, Paul
Reply | Threaded
Open this post in threaded view
|

Re: A plane intersecting a surface

Dear all

 

Thanks Christophe, Rafael and Stéphane for the first feedback;

 

Only obvious things in the code hereafter, but it highlights I guess what I would like to do (to cross section the surface); the results are not really noisy and their number is of about few hundred.

 

Concerning the Delaunay approach, I thought about it but I've been thinking a simplest solution may exist if I can plot the surface (interpolation from the grid) ?

 

Paul

 

###########################################

mode(0)

 

function [z]=saddle(x, y)

    z = x^2 - y^2

endfunction

 

function [z]=x_square(x, d)

    z = x^2 - d^2

endfunction

 

function [z]=y_square(y, d)

    z = y^2 - d^2

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)';

z = feval(x,y,saddle);

scf(0);

plot3d(x,y,z);

 

// obvious cases

// n = (0 1 0) then z = x^2 - d^2

d = 0;

z1 = x_square(x,d);

scf(1);

plot(x,z1);

 

// n = (1 0 0) then z = d^2 - y^2

d = 0;

z2 = y_square(y,d);

scf(2);

plot(x,z2);

 

 

 

 

-----Message d'origine-----
De : users [mailto:[hidden email]] De la part de Dang Ngoc Chan, Christophe
Envoyé : lundi 10 septembre 2018 09:15
À : Users mailing list for Scilab
Objet : [EXTERNAL] Re: [Scilab-users] A plane intersecting a surface

 

Hello,

 

> De : users [mailto:[hidden email]] De la part de Rafael Guerra

> Envoyé : samedi 8 septembre 2018 14:52

> 

> If your cloud of points behaves well enough, you can interpolate it first into a dense

 

If nobody is expert in this field, then I could invoke a memory when I was a student.

I've heard about an algorithm using intercept with tetrahedrons,

it was used for surface rendering.

 

So you might perform a Delaunay tessellation of your cloud,

determine which tetrahedrons are cut

and determine the coordinates of the intercepts.

 

Or ask some CGI  specialists.

 

HTH

 

Regards

 

 

--

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.

_______________________________________________

users mailing list

[hidden email]

https://urldefense.proofpoint.com/v2/url?u=http-3A__lists.scilab.org_mailman_listinfo_users&d=DwIFAw&c=0hKVUfnuoBozYN8UvxPA-w&r=4TCz--8bXfJhZZvIxJAemAJyz7Vfx78XvgYu3LN7eLo&m=P1df3Jy1yla9yY6mABJopK4mYcW7v-fqNLKDRWKGljw&s=Rx09ENmvHowtyMmodOUeRVEG2RZPKzK3Sy4zc7Mw0VM&e=

EXPORT CONTROL :
Cet email ne contient pas de données techniques
This email does not contain technical data

 


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: A plane intersecting a surface

In reply to this post by Christophe Dang Ngoc Chan
Paul,

The attached home-made example illustrates the very simple procedure that I
have suggested earlier (i.e., using distance to plane + cshep2d function).

<http://mailinglists.scilab.org/file/t495698/Scilab_plane_intersection_with_cloud_of_dots_RG.png>

PS: if you have multi-valued input data points then cshep2d will not work.

Regards,
Rafael




--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: A plane intersecting a surface

In reply to this post by Carrico, Paul
Hello Paul,

Le 10/09/2018 à 09:29, Carrico, Paul a écrit :

Dear all

 

Thanks Christophe, Rafael and Stéphane for the first feedback;

 

Only obvious things in the code hereafter, but it highlights I guess what I would like to do (to cross section the surface); the results are not really noisy and their number is of about few hundred.

 

Concerning the Delaunay approach, I thought about it but I've been thinking a simplest solution may exist if I can plot the surface (interpolation from the grid) ?

in your set of (x_i,y_i,z_i) 3D points, are the (x_i,y_i) organized on a grid, i.e. cartesian product [discrete values of x] times [discrete values of y] ?

S.

 

Paul

 

###########################################

mode(0)

 

function [z]=saddle(x, y)

    z = x^2 - y^2

endfunction

 

function [z]=x_square(x, d)

    z = x^2 - d^2

endfunction

 

function [z]=y_square(y, d)

    z = y^2 - d^2

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)';

z = feval(x,y,saddle);

scf(0);

plot3d(x,y,z);

 

// obvious cases

// n = (0 1 0) then z = x^2 - d^2

d = 0;

z1 = x_square(x,d);

scf(1);

plot(x,z1);

 

// n = (1 0 0) then z = d^2 - y^2

d = 0;

z2 = y_square(y,d);

scf(2);

plot(x,z2);

 

 

 

 

-----Message d'origine-----
De : users [[hidden email]] De la part de Dang Ngoc Chan, Christophe
Envoyé : lundi 10 septembre 2018 09:15
À : Users mailing list for Scilab
Objet : [EXTERNAL] Re: [Scilab-users] A plane intersecting a surface

 

Hello,

 

> De : users [[hidden email]] De la part de Rafael Guerra

> Envoyé : samedi 8 septembre 2018 14:52

> 

> If your cloud of points behaves well enough, you can interpolate it first into a dense

 

If nobody is expert in this field, then I could invoke a memory when I was a student.

I've heard about an algorithm using intercept with tetrahedrons,

it was used for surface rendering.

 

So you might perform a Delaunay tessellation of your cloud,

determine which tetrahedrons are cut

and determine the coordinates of the intercepts.

 

Or ask some CGI  specialists.

 

HTH

 

Regards

 

 

--

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.

_______________________________________________

users mailing list

[hidden email]

https://urldefense.proofpoint.com/v2/url?u=http-3A__lists.scilab.org_mailman_listinfo_users&d=DwIFAw&c=0hKVUfnuoBozYN8UvxPA-w&r=4TCz--8bXfJhZZvIxJAemAJyz7Vfx78XvgYu3LN7eLo&m=P1df3Jy1yla9yY6mABJopK4mYcW7v-fqNLKDRWKGljw&s=Rx09ENmvHowtyMmodOUeRVEG2RZPKzK3Sy4zc7Mw0VM&e=

EXPORT CONTROL :
Cet email ne contient pas de données techniques
This email does not contain technical data

 


-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: A plane intersecting a surface

In reply to this post by Rafael Guerra
The new attached snapshot is taken from a better angle to illustrate the plane cutting the 3D data cloud.

PS: The Scilab 3D plots all use orthographic views. An option to include perspective is missing.  

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

Scilab_plane_intersection_with_point_cloud_RG.png (94K) Download Attachment
Christophe Dang Ngoc Chan Christophe Dang Ngoc Chan
Reply | Threaded
Open this post in threaded view
|

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

In reply to this post by Carrico, Paul

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.
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Christophe Dang Ngoc Chan Christophe Dang Ngoc Chan
Reply | Threaded
Open this post in threaded view
|

Re: A plane intersecting a surface

In reply to this post by Christophe Dang Ngoc Chan

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.
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users