Code covered by the BSD License  

Highlights from
geom3d

4.85714

4.9 | 39 ratings Rate this file 396 Downloads (last 30 days) File Size: 438 KB File ID: #24484
image thumbnail

geom3d

by

 

19 Jun 2009 (Updated )

Library to handle 3D geometric primitives: create, intersect, display, and make basic computations

Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

| Watch this File

File Information
Description

The aim of geom3d library is to handle and visualize 3D geometric primitives such as points, lines, planes, polyhedra... It provides low-level functions for manipulating 3D geometric primitives, making easier the development of more complex geometric algorithms.
Some features of the library are:
- creation of various shapes (3D points, 3D lines, planes, polyhedra...) through an intuitive syntax. Ex:
  createPlane(p1, p2, p3) to create a plane through 3 points.

- derivation of new shapes: intersection between 2 planes, intersection between a plane and a line, between a sphere and a line...

- functions for 3D polygons and polyhedra. Polyhedra use classical vertex-faces arrays (face array contain indices of vertices), and support faces with any number of vertices. Some basic models are provided (createOctaedron, createCubeoctaedron...), as well as some computation (like faceNormal or centroid)

- manipulation of planar transformation. Ex.:
   ROT = createRotationOx(THETA);
   P2 = transformPoint3d(P1, ROT);

- direct drawing of shapes with specialized functions. Clipping is performed automatically for unbounded shapes such as lines or rays. Ex:
    drawPoint3d([50 50 25; 20 70 10], 'ro'); % draw some points
    drawLine3d([X0 Y0 Z0 DX DY DZ]); % clip and draw straight line

The library is divided into two packages:
* geom3d for general computation in 3d
* meshes3d for representation and manipulation of polyhedral meshes

Several functions require the geom2d package, available on the FEx.
     
Additional help is provided in geom3d/Contents.m and meshes3d/Contents.m files, as well as summary files like 'points3d.m' or 'lines3d.m'.

Acknowledgements

Iscoplanar.M inspired this file.

This file inspired Image Ellipsoid 3 D, Perspective Projection, and Intersect Plane Surf Ii.

Required Products MATLAB
MATLAB
MATLAB release MATLAB 8.3 (R2014a)
MATLAB Search Path
/
/geom3d
/geom3d/geom3d
/geom3d/meshes3d
Other requirements Some functions require the geom2d library, also on the File Exchange (ID 7844)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (109)
17 Dec 2014 Dexter

Hej, David:

Thanks for your sharing of this wonderful library.
I am working on draw a polyhedron, the first Brillouin Zone in Physics. Which specific function should I look in this library for my purpose?

Best wishes!

22 Nov 2014 matteo

Hi David,

thank you so much for your work (the intersectLineSphere), I am working on my master thesis and you really saved me.

21 Nov 2014 David Legland

Hi Matthew,
you can check the 'intersectLineMesh3d' function. It uses line representation in the form [x0 y0 z0 dx dy dz], where (x0,y0,z0) is the origin of the line and (dx,dy,dz) represents a direction vector of the line. Meshes are represented as vertices + faces, with vertices being Nv-by-3 array of vertex coordinates, and faces being either Nf-by-3, Nf-by-4 array of vertex indices, or cell arrays of vertex indices allowing more generic faces.

There is a function for reading meshes in OFF format ('readMesh_off'). Otherwise, you may want to check other contributions to import (and eventually convert) you files. See graph toolbox for example.

19 Nov 2014 Matthew Brandsema

Hello, I just found this toolbox and I am wondering if you can tell me if what I want to do is possible with these functions.

I have a vector equation of a line in 3D space, and I want to determine the position that it intersects with a mesh.

The mesh will be a mesh created in an external modeling program. (I can obtain all the information about its vertex coordinates from opening up a .obj file and extracting the information. )

I also need to know the normal vector of the mesh at the point of intersection.

Can this be done???

01 Nov 2014 Fritz  
21 Oct 2014 David  
13 Oct 2014 David Legland

@Abdulrahman,
sorry, I missed your comment... There is no direct function for testing if a point is inside a 3D circle. What you can do is 1) extract the plane containing the circle, 2) test if the point belong to the plane, and 3) if yes, compute coordinates of point in 2D plane, and test if the 2D point is within the 2D circle (you may need geom2d also).
hope this helps ?
regards

13 Oct 2014 David Legland

Hi Jin,
yes, there aer numerical precision issus with this function... I have submitted a new version, that is more tolerant for finding intersection at edges of faces. However sometimes too many intersections are obtained. You can remove them by using the 'mergeClosePoints' function from the geom2d toolbox.
regards,
David

09 Oct 2014 Jin

Hi David,

Really great tool! It helped me with a lot of problems. Thank you for sharing it with us.

Recently I'm using intersectLineMesh3d quite often. Here is an example I feel not all the intersections are found.

V =

0.5036 0.3333 0.3063
0.4102 0.1667 0.2098
0.1742 0.3820 0.1983
0.3505 0.4296 0.1688
0.3183 0.3623 0.3165
0.4246 0.2660 0.1180

P =

2 3 6
2 3 5
3 4 6
3 4 5
1 2 6
1 2 5
1 4 5
1 4 6

LINE =

0.3538 0.3302 0.2503 -0.8359 0.1600 -0.5250

This command
[inters pos inds] = intersectLineMesh3d(LINE,V,P) only gives me one intersection, but it should have two instead. Could you help me with that? Did I misunderstand the meaning of input and output? Thank you so much.

16 Jun 2014 David Legland

@Johannes,
The new set of faces refer to indices of vertices actually used by the faces. So it is necessary to use the following syntax:
[V2 F2] = mergeCoplanarFaces(vertices, K);
drawMesh(V2, F2);

The demo file is somewhat outdated, I will update it, and fix the doc of the mergeCoplanarFaces function.

Concerning the tolerance value, it is used for comparing normalised normal vectors of planes containing faces. There can be some numerical issues, but using 1e-4 should be fine.

regards,
David

16 Jun 2014 Abdulrahman

David

Is there a function to find 3d points within a 2d circle !?

Great submission THANK YOU

13 Jun 2014 Johannes

Dear David,
I was very happy to find geom3d as it provides a multitude of tools i need. However, I am concerned to have discovered some oddities in mergeCoplanarFaces(). I have Voronoi polyhedra that are clipped to fit inside e.g. a cube. Here's an mwe:

%%% code %%%

close all;
clear all;
clc;

v1 = [0.7945 0.2065 0.3770
0.7851 0.2019 0.3754
0.6767 0.1483 0.3603
0.6747 0.1566 0.3578
0.3994 0.0620 0.8142
0.3305 0.4372 0.8683
0.5055 0.3719 0.8659
0.3032 0.2532 0.8922
0.3026 0.1366 0.8615
0.0841 0.2280 0.8028
0.1441 0.3496 0.8355
0.1113 0.2820 0.8271
0.8097 0.2518 0.7474
0.8414 0.2398 0.4822
0.7089 0.3694 0.8072
0.7858 0.3168 0.3958
0.3059 0.4221 0.8716
0.1471 0.3505 0.8370
0.2376 0.3283 0.8869
0.2138 0.3265 0.8775
0.3857 0.0514 0.7753
0.5503 0.0944 0.4249
0.2406 0.1644 0.5060
0.3526 0.4066 0.3743
0.3545 0.4235 0.3775
0.2192 0.2369 0.4857
0.1307 0.1782 0.5772
0.0560 0.1827 0.6723
0.0610 0.1760 0.6639
0.4537 0.5529 0.4090
0.2929 0.4950 0.4783
0.5109 0.6260 0.6706
0.1959 0.4789 0.6226
0.1916 0.4725 0.5978
0.5015 0.6228 0.6758
0.5646 0.5852 0.4243
0.4909 0.5704 0.4153
0.5468 0.6032 0.4606];

v2 = [ 0 1.0000 1.0000
0 0.8719 0.8367
0 0.7113 0.9914
0.0336 1.0000 0.7575
0.1827 0.9701 1.0000
0.1884 0.9592 1.0000
0.1394 1.0000 0.8965
0.1859 1.0000 0.9718
0.1910 1.0000 0.9642
0 1.0000 0.7134
0.1239 1.0000 0.8761
0.1670 1.0000 1.0000
0 0.7024 1.0000
0.1259 0.8741 1.0000];

vertices = v1;

% triangular faces can be obtained from convhulln
K = convhulln(vertices);

% display polyhedron
figure(1); clf; hold on;
set(gcf, 'renderer', 'opengl');
drawPolyhedron(vertices, K);

K2 = mergeCoplanarFaces(vertices, K, 1e-4);

% draw the polyhedron
figure(2); clf; hold on;
set(gcf, 'renderer', 'opengl');
drawPolyhedron(vertices, K2);

%%% code %%%

For v1, mergePlanarFaces only merges faces if the PRECISION is set to values > 1e-5. This startles me, since this particular polyhedron is not modified by me, but comes straight from convhull. For v2, the result seems to be completely garbled.

Am I doing something wrong here?

Cheers,
Johannes

10 Jun 2014 David Legland

Hi Chien Ting,
thanks, this is fixed! Yes, if you have any more updates, I an include them. Do not hesitate to contact me by email (on my author page)

regards,
David

03 Jun 2014 Chien Ting CHIN

I've been using your toolboxes for years, thanks for some great tools.

In drawEdge3d.m, line 29, I believe is a bug, it would never get executed. The correct line might be:
elseif nargin >= 6;

Also I've complete the documentation of drawEdge3d.m for my group's use. If you want see/use my version, let me know!

Thanks for some great tools!

03 Jun 2014 Chien Ting CHIN

Sorry I pasted the original line. The correct line might be:
elseif nargin >= 6;

24 Apr 2014 Jasper  
26 Feb 2014 David Legland

@John,

thanks for reporting! After checking, it appears that in your example, the plane is passing exactly through one of the vertices of the mesh. This can sometimes happen, and I do not know yet how to handle these cases. A simple workaround is to slightly change the position of the plane (or the mesh), by adding a random value at the origin:
plane = createPlane( randn(1,3)*1e3, [ 0 1 1 ] );

Hope this helps ?

24 Feb 2014 John Anderson

Hi David, great submission,

I think I have noticed a slight problem with calls to intersectPlaneMesh.m and wasn't sure if this had previously been flagged.

When I exectue the following to calculate the intersection of a plane with a sphere I find that when N = 10 I get the error

Error using horzcat
Dimensions of matrices being concatenated are not consistent.

Error in intersectPlaneMesh (line 85)
polyEdgeInds = [polyEdgeInds currentEdgeIndex]; %#ok<AGROW>

Error in testIntersectingMesh1 (line 27)
polys = intersectPlaneMesh( plane, vertices, faces );

However, when N=15 there are no exectuion problems

Code:

% create a plane
plane = createPlane( [ 0 0 0 ], [ 0 1 1 ] );

% create sphere
N = 10;
[ x, y, z ] = sphere( N );

% triangulation
tri = DelaunayTri( x(:), y(:), z(:) );

%vertices
vertices = tri.X;

%faces forming convex hull
faces = convexHull( tri );

cla;
% draw solid
patch( 'Vertices', vertices, 'Faces', faces, 'facecolor', 'y', 'edgecolor', 'y', 'facealpha', 0.1 );
hold on;

% intersecting polygon
polys = intersectPlaneMesh( plane, vertices, faces );

% draw polygon
drawPolygon3d(polys, 'color', 'b', 'LineWidth', 2);

% fill polygon
fillPolygon3d( polys, 'r', 'faceAlpha', 0.2 );

Best wishes,

27 Jan 2014 David Legland

Hi Germano,
thanks, I will check it! I remember I used normalisation before cross product for improving stability, but I have to check again. Anyway, should be fixed in next release!

21 Jan 2014 Germano Gomes

Great library, i use it all the time, thanks.
A small correction for:

function normals = faceNormal(nodes, faces)

if isnumeric(faces)
% compute vector of first edge
v1 = nodes(faces(:,2),1:3) - nodes(faces(:,1),1:3);
v2 = nodes(faces(:,3),1:3) - nodes(faces(:,1),1:3);

% normalize vectors
% v1 = normalizeVector3d(v1); % remove
% v2 = normalizeVector3d(v2); % remove

% compute normals using cross product
normals = cross(v1, v2, 2);
normals= normalizeVector3d(normals); % correction GTG 2014

else
normals = zeros(length(faces), 3);

for i=1:length(faces)
face = faces{i};
% compute vector of first edges
v1 = nodes(face(2),1:3) - nodes(face(1),1:3);
v2 = nodes(face(3),1:3) - nodes(face(1),1:3);

% normalize vectors
% v1 = normalizeVector3d(v1); % remove
% v2 = normalizeVector3d(v2); % remove

% compute normals using cross product
normals(i, :) = cross(v1, v2, 2);
end
normals= normalizeVector3d(normals);% correction GTG 2014
end

---------------
Always better to normalize at the end otherwise the final normals are not normalized because of roundoff etc errors in the cross computation, this also has the advantage of eliminating 2 calls to normalizeVector3d, making faceNormal a tad faster.

08 Jan 2014 samart  
13 Dec 2013 Khaled Khairy

After some minor initial tweaking, I find this a very helpful mission. Despite its size all parts seem to work seamlessly together. Thanks David.

08 Nov 2013 Shaoshuai

Very helpful! Thanks for sharing!

29 Oct 2013 J.R.! Menzinger

Very useful.

15 Oct 2013 David Legland

Hi Cedric,
thank for your interest!
Actually the centroid function is provided in the "geom2d" library, available on the Fex as well. I try to minimise dependencies between both submissions, but apparently some are still remaining...
You can also replace the centroid call in the demo script by the following:

faceCenter = mean(vertices(faceVertices, :), 1);

This should work as well.

Regards,
David

13 Oct 2013 Cedric

Hi David,
Thank you for your code!
I have a small problem.
I am not able to compile the DrawSoccerBall.m demo function.

there's a parameter missing when calling centroid function since the definition of centroid function is function centroid = polyhedronCentroid(vertices, faces)
Error message when I compile:
Undefined function 'centroid' for input arguments of type 'double'.

Error in drawSoccerBall (line 56)
faceCenter = centroid(vertices(faceVertices, :));

20 Aug 2013 David Legland

@Yong Min Yoon:
I just realized that I already implemented a "polygonArea3d" version, that is included in the geom3d archive... So I think this should do the job!
Regards,
David

20 Aug 2013 David Legland

@Yong Min Yoon:
It's great if ou could make the fucntion work better! Could you sent me the code so that I can include it in a future release ?
For polygon area, there is a function "polygonArea" in geom2d, that works for planar polygons (beware that it returns a signed area). You can transform 3D polygon to planar polygon by using "projPointOnPlane" function. I think this should be ok for computing area ?
regards

16 Aug 2013 Yong Min Yoon

@David
I modified your intersectPlaneMesh function code and make it work to find the multiple contours.

As a next step to make up my algorithm, I want to know that if there's a area computation of the polygon(the closed contour) in your works or not.

If you let me know, this will be great helpful for me. Thank you

14 Aug 2013 Yong Min Yoon

@David
Your code is totally what I want. So, I'm just going to try to fix your code. Can you explain what the problem is in your function?

13 Aug 2013 David Legland

@Yong Min Yoon:
first of all sorry for the delay to reply;
The function you are working with at the moment is not fully implemented. It expected to returns a set of polygons, but it currently returns only one. I have plan to extent to multiple polygons, but it will require some work (one idea would be to recursively link adjacent vertices until the first vertex is reached, then to start a new polygon). So, no easy solution for the moment, but stay tuned...

13 Aug 2013 Yong Min Yoon

@ David Legland

Hi, first of all, I'm very thankful for your codes. I have some questions in your "intersectPlaneMesh function".

I try to make the some cross section information of the vertebrae. Your code will be used to make up that.

when I used that codes with a vertebrae mesh data, 2 contours are expected to show the cross sections. However, only one contour was shown in the figure.

I think that your polys variable which is return in intersectPlaneMesh, is used to represent the multiple contours to show intersected region by the plane.

If it is right, can you explain how to set the multiple contours by the plane?

Images are in the below URL
URL: http://gall.dcinside.com/board/view/?id=medicalscience&no=162489

P.S. intersectionsPoints in your intersectPlaneMesh function seems to exactly show intersected points on the mesh data.

09 Aug 2013 Evan

@David:
Thanks for the response. What you've provided here has already been extremely helpful anyway.

I also had a small suggestion for computing vertex normals. If I'm reading it right, the function normalizes the facet normals before averaging them. But it turns out to be more accurate to weigh by area and angle (http://www.bytehazard.com/code/vertnorm.html). For area, I believe that just means simply not using normalizeVector3d when calling faceNormal. It shouldn't be too difficult to code the angle part on my own, but I'm just throwing it out there in case you were interested.

@Sven:
Awesome, thanks for the heads up.

09 Aug 2013 Sven

@Evan,
Please check out http://www.mathworks.com/matlabcentral/fileexchange/43013-unifymeshnormals
It should do what you're looking for.

09 Aug 2013 Sven

@Evan,
I'm working on something to do exactly this. Will keep you posted here.

08 Aug 2013 David Legland

@Evan:
depending on how the meshes are constructed, the orientation of normals may be consistent (ie, all pointing outwards of the mesh), or not.
The function "faceNormal" does not make any assumption on this, and simply computes the normal of each face.

However, some other functions require correct orientation of the normals. This is the case for vertexNormals, for example.

I do not have perfect solution... I have written "checkMeshAdjacentFaces" to check normal orientation consistency of a given mesh, and it can be safe to use it also. But I do not have direct solution to re-orient face normals in given direction.

hope this helps ?

07 Aug 2013 Evan

Hello David, I had a question about calculating normals. In the function faceNormal, you mention that no orientation is defined. But in vertexNormals, faceNormals is used to calculate the normals at the vertices. Won't this lead to some problems, or am I missing something?

I'm specifically interested in calculating normals of a surface (facet or vertex) that are oriented in the same direction, which I understand is difficult and non-trivial problem. Do you have any suggestions? I noticed the checkMeshAdjacentFaces function, but I'm not sure if that's useful for "re-orienting" specific normals.

06 Aug 2013 Yong Min Yoon

Hi, first of all, I'm very thankful for your codes. I have some questions in your "intersectPlaneMesh function".

I try to make the some cross section information of the vertebrae. Your code was used to make that.

when I used that codes with a vertebrae mesh data, 2 contours are expected to show the cross sections. However, only one contour was shown if the figure.

I think that your return(polys) is used to represent the multiple contours to show intersected region by the plane.

Is my thought right? If it is right, can you explain how to set the multiple contours by the plane?

Images are in the below URL
URL: http://gall.dcinside.com/board/view/?id=medicalscience&no=162489

P.S. intersectionsPoints in your intersectPlaneMesh function seems to exactly show intersected points on the mesh data.

01 Aug 2013 David Legland

Hi Benjamin,
you're not wrong, there was a mistake in the function surfToMesh...
I have submitted a fixed version, that should appear soon.
Regards,
David

31 Jul 2013 Benjamin

Dear David,

I'm a bit of a novice user, so it's likely this is a flub on my part, but I figured it worthy of mention.

I've been attempting to use the cylinderMesh function, but with little success. I even tried the "Draw a rotated cylinder" example, but all I get is
Index exceeds matrix dimensions.

Error in surfToMesh (line 91)
z = z(1:n1, 1:n2);

Error in cylinderMesh (line 63)
[vertices faces] = surfToMesh(x, y, z);

I'm not entirely sure what is causing this, any ideas?

Benjamin

10 Jun 2013 David Legland

Hi Khaled,
thanks for reporting problem. I will check behaviour for points given in spherical coords. The function works fine with points given as cartesian coordinates (assuming sphere is centered), so I hope this could help using the function.

Regards,
David

04 Jun 2013 Khaled Khairy

With all respect for this valuable submission. I do have some slight critique. I downloaded geom3d to "quickly" check against my code for the determination of the spherical angle, as the function "sphericalAngle" promised. I only called sphericalAngle exactly as indicated in its help. I ended up wasting two hours of quality time ploughing through various fixes. [1] calls to functions whose names have been changed. Sorry , but I just stopped keeping track(some have already been mentioned in these comments) and [2] createPlane, was called from within sphericalAngle apparently with the wrong dimensions. [3] intersectLinePlane complained about matrix dimension mismatch on line 71. Honestly, I consider perseverance one of my strenghts, but I gave up at this point, as I can't say how much more time this really simple operation would/should cost me should I continue with geom3d. Sorry, but a smooth program flow is really important, especially in extensive code submissions such as this one.

28 May 2013 David Legland

Hi Orestis,

You can check the function 'intersectPlaneMesh', that computes the 3D polygon resulting from the intersection of a polygonal mesh with a plane. The algorithm used is simple, so there can be numerical issues. Then you can compute the area with the function 'polygonArea3d'. By considering several planes, you should be able to do what you want.

There is no function in geom3d for reading STL files. So you should either write yout own, or find one.

regards,
David

27 May 2013 Orestis

Dear David

Could you let me know if your library can be used i order to slice a 3d stl file with arbitrary planes?

The aim is to calculate the area bounded by the stl file in each slice.

With best regards

11 Mar 2013 David Legland

Dear Dmitri,

thank you for your feedback! Actually, many functions require the geom2d contribution to work properly. I try to make geom3d as independent and self-contained as possible, but downloading geom2d as well is strongly recommended. I will update the different functions you mentioned to remove/reduce dependencies.

Regards,
David

08 Mar 2013 Dmitri Kamenetsky

Great package! Some minor bugs:

1. circleToPolygon function is missing. It is called in demoRevolutionSurface,
drawTorus and torusMesh.
2. vectorNorm is now vectorNorm3d. It is used in intersectLineTriangle3d, triangleArea3d.
3. drawEdge is now drawEdge3d. It is used in demoGeom3d.

08 Mar 2013 Dmitri Kamenetsky

One more bug: circleToPolygon function is missing. It is called in demoRevolutionSurface, drawTorus and torusMesh.

13 Feb 2013 panagiotis

@Sven
I am sorry , but I didn't know that point of view for plane creation. My background bibliography doesn't include such information.Thank you for informing me.

12 Feb 2013 Sven

@panagiotis: I highly recommend reading *again* the documentation for createPlane - it is very short and very clear, but your comment shows you still have a misunderstanding.

You seem to think that "3 points are a plane". This is not true. A plane is like a coordinate system - it has:

- An origin POINT
- An X-direction VECTOR
- A Y-direction VECTOR

You can *create* a plane from 3 points, but a plane is *not* just 3 points. When you read the help for createPlane, you will see the first example which shows you *exactly how* to create a plane from 3 points:

PLANE = createPlane(Pt1,Pt2,Pt3)

You will see that PLANE is NOT EQUAL to [Pt1,Pt2,Pt3]. There are *very* good reasons for this which you will understand once you read the help for createPlane.

12 Feb 2013 panagiotis

@David
You are right David , but didn't thought at all that there would be a misunderstaning for defining the plane by three points!I think it is commonly accepted. Haven't heard before for direction of plane from the third point.
Anyway, came up with something new.

In createPlane , you take into account only the 2nd & 3rd point for the cross product.You don't create vectors like p1-p2 & p1-p3.It looks to me as if you define the [0 0 0] , point for all planes. The normal in some cases looks wrong to me and didn't find any help.. If I use the vectors , results are fine.

plane=[1 1 1 ,2 2 1, 3 3 1]
planeNormal(plane)
ans = 0 0 0

I think that this should be a vector like [0 0 a].

12 Feb 2013 David Legland

@panagiotis:
Apparently you did not really even try to get som help by yourself. Either "help planes3d", "help projPointOnPlane", or even "help geom3d" gives you the definition used in the library for representing planes.

Giving a low rating just because you did not want to read the doc is somewhat unfair...

By the way, many thanks to Sven for the support!

Regards,
David

11 Feb 2013 panagiotis

Sven,thank you for informing me. The truth is that I didn't got that deep in help..
As soon as I use the library again , I will reconsider the 1-star rating!
Thanks!

10 Feb 2013 Sven

Hi Panagiotis, the help for projPointOnPlane says "... PLANE is a [N*9] array (see createPlane for details)". The help for createPlane is very clear about the definition.

Anyway, I'm glad that your issue has now been identified. One further note is that your "other" direction vector was specified as [0 1 1], which isn't "normalised" (ie, it has a magnitude of more than 1). This won't affect projPointOnPlane but it can affect other functions such as planePosition if you intend on using them.

Hope this makes your code work again, and if so, perhaps you might reconsider the 1-star rating :)

09 Feb 2013 panagiotis

@ Sven
There was no such notice for this definition. I didn't know that the third point defines the plane's direction .

08 Feb 2013 Sven

@panagiotis:

Your definition of the plane does not make sense:

plane=[0 1 0 0 1 1 0 0 0];

Notice that elements 7:9 are all zero. Elements 7:9 give the plane's Y-axis direction vector. A vector of all zeros has no direction at all, therefore it is impossible to project a point onto this plane, because the plane doesn't exist.

08 Feb 2013 panagiotis

@Sven(strong defense for this library!)

point=[0 0 0];
plane=[0 1 0 0 1 1 0 0 0];
projPointOnPlane(point,plane)
ans = NaN NaN NaN

07 Feb 2013 Sven

@panagiotis:

>> projPointOnPlane([100 100 0],[0 0 0 1 0 0 0 1 0])

ans =

100 100 0

... seems fine to me. Please provide a reproducible error.

07 Feb 2013 panagiotis

The main problem I found is in projectionPointOnPlane in geom3d library. We tried to project a point on a plane on which the point already belongs to. The result that came up was NaN instead of the point's coordinates. That led to some problems with other functions like intersection..

Panagiotis

06 Feb 2013 Sven

@panagiotis: please describe the actual issue(s) you found. I would bet that if you found a reproducible error, David would very quickly address it. Just saying "there are errors with intersections" isn't actually very helpful, and doesn't let anybody improve the package above the 1-star you feel it deserves at the moment.

06 Feb 2013 panagiotis

Dear David,
I am using your library for a couple of months and I have found some errors especially on projecting points and intersection..
Thought you might find it helpful.

Panagiotis

28 Jan 2013 David Legland

Hi Charles,
I have created a project page on sourceforge call "matgeom", that gather geom2d and geom3d. There are some forum facilities, so this should be easier to post data there I think.
http://matgeom.sourceforge.net/

Regards,
David

20 Jan 2013 Charles

Dear David,
I'm trying to use surfToMesh and then meshSurfaceArea function to calculate the approximate area covered by the dots with X and Y coordinates specified on X-Y plane. However, the area returns the value that does not make sense to me. I wonder if it's possible for me to attach the matrix here and get some help why it does not work in my code. The dots I have represents nearly a semi-circle and the result gives me the area close to a full circle of that diameter. Is there some tolerance setting on meshing? Thank you very much!

Best,
Charles

14 Dec 2012 Student Uni

HI Sir
Basically i have to draw that 3D rectangle and divinde its faces in
different infinitesimal areas, and divide that rectangle in small volums.
then i need to take the coordinates of Vertices of areas and volumes and
also coordinates of their central point
clc
clear all
close all
for x=0:1:2;
for y=0:1:2;
plot3([x,x],[0,2],[y,y]);
hold on
end
end
for x=0:1:2;
for y=0:1:2;
plot3([0,2],[x,x],[y,y]);
hold on
end
end
for x=0:1:2;
for y=0:1:2;
plot3([y,y],[x,x],[0,2]);
hold on
end
end
i need something like this but this is square and i dont know how to get
the required co ordinates
I appreciate in advance your kind attention
Best Regards

13 Dec 2012 Hoai-Nam LE

Many thanks David!!!
Nam

27 Nov 2012 David Legland

@Hoai-Nam,
I have submitted a new version with some new functions you should find useful:
* functions 'reversePlane' and 'reverseLine3d', that avoid explicitely permutting the directions vectors
* function intersectPlaneMesh, that returns the intersection between a plane and a 3D mesh. Not a full clipping function yet, but this should help!

regards,
David

21 Nov 2012 Yuri

David, thanks a lot! Yuri

13 Nov 2012 David Legland

@Yuri,
The function sphericalAngle assumes the input points are located on the surface of the unit sphere, that is not the case in your problem.
I think you should use the "angle3Points3d" function, that computes the angle between 3 points in space, using the second point as center of the angle.

anglePoints3d(p1, j1, l1)
ans =
1.5708

08 Nov 2012 Hoai-Nam LE

Dear David,

Thanks for your very useful toolbox!

I've used some of your functions and I've found some things :

+ My problem :
I have a 3D point set. I create a minimum convex hull of this point set, using convhulln, so I have a list of facets of the convex hull.

I'd like to clip this convex hull by a plane. And this leads me to your function.

+ There 2 ways to do this with your toolbox :

1/ Use the function clipConvexPolyhedronHP. It works for my problem!
However if we want to have a side of convex hull (below or above to the plane) to clip, we must have a right plane input!!!

For example:
chX = convhulln(Pts); %Pts : 3D point set

% Select 3 points on the plane z=vClip
P=[0 0 vClip;
0 1 vClip;
1 0 vClip];
% Create plane z=vClip
plane=createPlane(P);

% This function createPlane gives me : plane=[0,0,0.7,0,1,0,1,0,0]
% Then the function clipConvexPolyhedronHP returns me a wrong polyhedron. So I permute the 2 direction vectors of the plane : plane=[0,0,0.7,1,0,0,0,1,0] and it returns me a polyhedron with the right side that I desire.

So do you have a tip to use the function createPlane to have a good result (so I dont have to permute the 2 direction vectors of the plane) ?

2/ Use the function polyhedronSlice. It doesn't work because it doesnt return only the intersection points of the plane and the polyhedron, but also some points that are wrong (see the image https://dl.dropbox.com/u/8702546/WEB/image_1.png).

However, if the function works, I can control the side of polyhedron that I want to clip by a simple command of Matlab!

So I feedback you this to correct the function polyhedronSlice!

I dont know if someone has already asked about these things...

Regards,
Nam

06 Nov 2012 Yuri

David, very nice and useful package. Thanks.

I have a question about sphericalAngle function. I defined the following plane:
>> p1 = [1, 1, 1]; p2 = [100, 100, 1]; p3 = [1, 100, 1];
>> pl = createPlane(p1, p2, p3);

Then a point outside the plane and the projection point onto the plane:

>> l1 = [40,40,40];
>> j1 = projPointOnPlane(l1, pl);

I'd expect the angle between l1, j1 and any point on the plane to be pi/2. However, this is what I get:

>> alpha_rad = sphericalAngle(l1, j1, p3)

alpha_rad =

1.5745 (~ pi/2)

>> alpha_rad = sphericalAngle(l1, j1, p1)

alpha_rad =

0 + 2.1073e-08i (= 0)

>> alpha_rad = sphericalAngle(l1, j1, p2)

alpha_rad =

3.1416 - 0.0000i (= pi)

Can you, please, clarify the way the function works?

Thanks

25 Oct 2012 marc

Very nice. thanks

16 Oct 2012 AP

Great tools.

05 Jul 2012 Henri SHI  
29 Jun 2012 Myles  
18 Jun 2012 Haochen Tang

Hi David,

Thank you. I miss understand the definition. I thought it was the origin and the end point. Thanks a lot.

Best,
Haochen

18 Jun 2012 Eduardo

Hi David,
I´ll try to use some of that functions. First of all i want to do only one slice in order to get the intersection points, only after that i´ll do it with the other parallel planes.
Thanks.

Regards,
Eduardo

18 Jun 2012 David Legland

Hi Eduardo,
What you want to do seems possible, but to my knowledge there is no direct method.
There are some functions in the "meshes3d" sub-package that could help you, like intersectLineMesh3d or polyhedronSlice. You can also use intersectEdgePlane -> you can obtain the set of intersection points between a plane and all edges. If you want to represent the slice as 3D polygon, you will still need to link points in appropriate order.

@Haochen:
Beware of the definition of 3D line. In the toolbox 3D line is represented by [x0 y0 z0 dx dy dz]. In order to define a vertical line you should use [0.6 0 0 0 0 10]. In this case you will get the correct result.

Regards,
David

18 Jun 2012 David Legland

Hi Eduardo,
What youwant to do seems possible, but to my knowledge there is no direct method.
There are some functions in the "meshes3d" sub-package that could help you, like intersectLineMesh3d or polyhedronSlice. You can also use intersectEdgePlane -> you can obtain the set of intersection points between a plane and all edges. If you want to represent the slice as 3D polygon, you will still need to link points in appropriate order.

Regards,
David

14 Jun 2012 Haochen Tang

Hi David,

Thanks a lot for this great tool. I am using the intersectLinePolygon3d function. Some of the results seem not make sense. For example:

pts3d = [-1 -1 5;1 -1 5;1 1 5];
line=[0.6 0 0 0.6 0 10];
[inter inside] = intersectLinePolygon3d(line, pts3d);

And I got:
nter =

0.9000 0 5.0000

inside =

1

Should we expect to get inter = 0.6 0 5?

Thanks so much,
Haochen

11 Jun 2012 Eduardo

Hi David
I´m trying to slice a triangular mesh. I´m not an expert in Matlab and i´m having some difficulties in implementing this. I have the edges, and i´m trying to get the intersection points with several parallel planes...
Can you give me a hand?
Thank´s.

Regards,
Eduardo

15 May 2012 David Legland

Hi Amanda,
I do not know similar library for scilab. I suspect you need to re-code everything...
If you are looking for a free alternative, I can mention a port to octave of libraries geom3d+geom2d:
http://octave.sourceforge.net/geometry/index.html

regards,
David

11 May 2012 Amanda

Hi, David
Having been using this package for years. Really love it. Thanks for your gr8 job. Now I'm researching using Scilab to replace Matlab for some application. Searched scilab but didn't find geom3d available there. Any suggestion?

11 May 2012 Amanda  
10 May 2012 Tianjin Univ  
03 Mar 2012 David Legland

Hi Nir,
did you mean 3D polygon or 3D mesh ? It seems you speak about mesh, right ? In this case you case use function "intersectLineMesh" (contained in the updated submission). This should do the job.

@Benjamin: I have also added a function polygonArea3d in the update.

29 Feb 2012 Nir

Hey David,
It seems that the functions:
intersectRayPolygon3d
intersectLinepolygon3d
Don't work properly.
I have a rectangular 3D polygon (fairly simple, like an elongated cube) defined by both the vertices and the planes of it's sides.
A line or a ray, with a location in the middle of it in ANY direction should generate an intersection point. This does not work. Not for a line nor a ray.

Vertices:
0 0 0
0 0 2.7000
5.7000 0 0
5.7000 0 2.7000
5.7000 7.0000 0
5.7000 7.0000 2.7000
0 7.0000 0
0 7.0000 2.7000

Ray
Origin: 2.8508 3.5013 1.3500
Direction: -0.1 0.1 1

Help?

27 Feb 2012 Hannes

Thanks mate. Beer flies out to you :)

24 Feb 2012 David Legland

Hi Ben,

you can check the function "polyhedronSlice". It computes the 3D polygon that results from the intersection of a polyhedron with a plane. Example of use:
[nodes faces] = createCube;
plane = createPlane([.5 .5 .5], [2 3 4]);
poly = polyhedronSlice(nodes, faces, plane);
figure; hold on;
drawMesh(nodes, faces, 'facealpha', .5);
fillPolygon3d(poly, 'm');

There is no (not yet...) function for computing area of a 3D polygon. However it can be obtained by projecting on the intersecting plane:
poly2d = planePosition(poly, plane);
area = abs(polygonArea(poly2d))

Regards, David

16 Feb 2012 Benjamin Sanderse

Hi David,

Thanks for your very nice work. I am searching for a routine that can find the (area of a) polygon that results from the intersection between a plane and a cube. Would this be possible with your library?

Regards, Ben

10 Feb 2012 David Legland

Hi Patrik,

no, you haven't missed it, and it could be useful. I will add it in a future release.

regards

08 Feb 2012 Patrik Eschle

distanceLine3d works on infinite lines a+t*b. Sometimes (collision detection) the distance between line segments AB and CD is required. A distanceSegment3d would be useful.

The same is true for points (distancePointSegment ).

Or did I just miss it?

07 Feb 2012 Patrik Eschle

Nicely written, good comments - thank you!

13 Dec 2011 Shuhao Cao  
07 Aug 2011 Sakshi  
20 Jul 2011 David Legland

@Qais,
you should install both geom2d and geom3d, typically in the same directory. Then you have to add directories to the path, using te 'addpath' function, or the 'File->Set Path' menu). You need to add at least the directories 'geom2d', 'geom3d', 'polygons2d', 'meshes3d', and eventually 'polynomialCurves2d'.

Some demos are provided with each archive, you can run them directly, and have an idea on the syntax.

I will also try to add some info on the matGeom web page (that gathers geom2d and geom3d), at http://matgeom.sourceforge.net .

Hope this helps

19 Jul 2011 Qais AlKhazraji

Dear David,
I am just start using Matlab in my research for geometrical computation. I appreciate if you would like to demonstarate how to use this library(3d or 2d) since each function I used give me syntax error about other functions as undefined. I am a biginner in Matlab and may be my question is dump question. I am very thankful.

22 May 2011 David Legland

@Khaled and Qais,
What I had in mind was to use the following procedure :
* compute intersection of line and polygon plane
* project polyingon vertices onto the plane (ie, compute their 2D coords)
* project intersection point onto the plane
* check if 2D intersection point is within the 2D polygon.
The following piece of code should do the job (assuming poly3d is N-by-3 and line is 1-by-6):
plane = createPlane(poly3d(1:3, :));
inter = intersectLinePlane(line, plane);
poly2d = projPointOnPlane(poly3d, plane);
pInt2d = projPointOnPlane(inter, plane);
inside = xor(isPointInPolygon(pInt2d, poly2d), polygonArea(poly2d) < 0);
inter = inter(inside, :);

Both geom2d and geom3d are required. I will add functions 'intersectLinePolygon3d' and 'intersecRayPolygon3d' in the next release.

regards

21 May 2011 Qais AlKhazraji

Dear David:
I read your comment for Khaled question. It is not clear to me about what you said"The best is to use 3D plane, and eventually to project all points on the plane, then to use function for 2D polygons." Which points we should project on 3D plane and from my understanding line intersect a polygon in one point.
Thank you in advance

16 Feb 2011 David Legland

Hi Khaled,
to create a 3D polygon, you simply need to concatenate their coordinates into a N-by-3 array, e.g. poly=cat(1,p1,p2,p3);
There are no fucntions for intersections of 3D line with 3D polygon. The best is to use 3D plane, and eventually to project all points on the plane, then to use function for 2D polygons.

15 Feb 2011 Khaled

Dear eng. David,

Thanks, this powerful toolbox is very helpful.

I tried to create a 3D polygon from four or three 3D points, and couldn't. So instead I tried to create a plane, but when finding an intersection with a 3D line the answer will not be accurate as the plane was polygon.

Is there a possible way to create a 3D polygon from four or three points? and
Is there a possible way to find a intersection between a 3D line and 3D polygon?

your efforts are appreciated.

19 Jan 2011 Brandon  
01 Dec 2010 Brady

Incredibly helpful toolbox. Mostly intuitive function calls, and thorough documentation to make integration easy.

06 May 2010 Siyi Deng

very helpful.

11 Dec 2009 Sven

Excellent package for spatial manipulation. Extensive toolbox, that handles input intuitively. The only suggestion I have is, for completeness, to document some of the methods of calling functions for the more obscure functions. Eg., localToGlobal3d has one documented method of calling, but under the hood it parses input in many different ways. I find myself typing "edit" instead of "help" just to make sure I'm using things efficiently.

07 Dec 2009 David Legland

The format for cylinder is as follow: [xv yv zy R xc yc zc], with:
[xv yv zv] the direction vector of the cylinder,
R the radius of the cylinder
[xc yc zc] the origin of the cylinder.
the L parameter, which is given in help, is not used.
The cylinder is considered infinite, please use "linePosition3d" if you need finite cylinder.
Please note that this format is not consistent with the function "drawCylinder", for example. Therefore the function "intersectLineCylinder" is likely to change in a future release.
DL

03 Dec 2009 Paklah Abd Rahman

Hi,

I tried the "intersectLineCylinder" what are the definitions of: xa, ya, za, and L for cylinder??..Please anybody

01 Sep 2009 David Legland

Hi,
I have uploaded a new version, that fixes reported bugs, and have added some demos. Demos are available as published m-files, and can be accessed through the file presentation above. They demonstrates some usage possibilities.

Most drawing functions follow the 'plot' syntax. Data are given first, then optional couples of parameter name and value. If not implemented, try using "h=drawMyShape(...)", then set up options using "set(h, 'color', 'k')" for example.

For bug reports and suggestions, the best is to sent me directly and e-mail (can be found in my author page, replace "DOT" by ".").

DL

03 Aug 2009 Jon Sporring

This package aims at solving visualization problems, I often have. However, I have not been able to I'm trying to use this package, but experiencing a number of problems:

1. How are the colors controlled of the objects drawn?
2. I want to draw cylindrical representation of arrows, since quiver3 together with drawPlane does not work well, possibly due to lack of control of drawing order. Instead I want to draw thick arrows, and I'm using drawCylinder.m. However, in line 89 is used the non-existing function local2Global.m. It appears to be replaceable by localToGlobal3d, if localToGlobal3d line 29 is replaced with "center = varargin{1};".

For suggestions to corrections, etc. what is the best procedure?

Thanks, Jon

17 Jul 2009 Nathan Thern

This library and its companion, geom2d, are shining jewels - indispensable for doing non-trivial calculations on objects in 2d and 3d space.

By the way, leave all the files in the geom3d directory and add the directory to your path. Then you can type "help geom3d" or "doc geom3d" and get properly linked help text in the command window or the help window. The Contents.m file is the proper "MATLAB way" of providing documentation.

19 Jun 2009 us

it is too bad that this (apparently) full-fledged software package does not come with an appropriate and modern help/intro/example module called a published m-file...

as of now, it is somewhat tedious for a naive user to have to wade through the help sections of the individual help files (or the densly written contents.m) - and it will (most likely) not give your package as much attention as it probably deserves...

us

Updates
26 Aug 2009

some bug fixes, demo files as published m-files.

31 Aug 2009

add missing demos

29 Jul 2010

fixed bugs in rotations and in drawing of some shapes, added rotation by Euler angles, various update in doc and in code

14 Jan 2011

Add new functions for meshes and polyhedra, and for 3D transforms (rotations, basis transform). See file changelog.txt

30 Jun 2011

use degrees for drawing shapes, split lib into packages 'geom3d' and meshes3d

13 Oct 2011

various code cleanup and speed improvements, mainly contributed by Sven Holcombe (many thanks!)

02 Mar 2012

bug fixes, new functions for spherical polygons, for reading meshes (off foramt), for point-edge distance, and for computing surface area of 3D polygons, meshes, or ellipsoids.

05 Apr 2012

fix bugs in demos

27 Nov 2012

several bug fixes, new functions (fitLine3d, parallelPlane, reversePlane, intersectPlaneMesh, meshCentroid...), and fonctions for creating meshes from patches or geometric primitives.

08 Jul 2013

minor bug fixes, new functions for manipulating polygons, speed improvements (thanks to Sven Holcombe!)

02 Aug 2013

fixed bug in surfToMesh

11 Jun 2014

fix various small bugs

13 Oct 2014

added trimMesh function (that removes non connected vertices from a mesh), update mergeCoplanarFaces, and update tolerance management in interesectLineMesh3d

Contact us