File Exchange

- You will see updates in your activity feed
- You may receive emails, depending on your notification preferences

version 1.24 (922 KB) by
David Legland

Geometric computing library for 3D shapes: meshes, points, lines, planes...

4.9

153 Downloads

Updated 26 Sep 2019

View Version History1.24 | new functions for processing 3D meshes (repairing, simplification), renaming of inertia ellipsoid into equivalent ellipsoid. |
||

1.23.1.1 | package as mltbx file |
||

1.23.1.0 | package as toolbox. |
||

1.23.0.0 | add meshes3d/meshFaceAreas.m, add geom3d/isTransform3d.m (thanks to Oqilipo), updated return type of several drawing functions, fixed bug in orientedBox3d |
||

1.22.0.0 | added several functions fore reading/writing 3D meshes in PLY and OFF formats, and new functions for working with 3d meshes (distancePointMesh, isPointInMesh). |
||

1.21.0.0 | integrate several contrinutions (oqilipo, Roozbeh). add clip/split/concatenate meshes, new transforms from vectors. |
||

1.20.0.0 | remove deprecated functions, update some display functions |
||

1.19.2.0 | fix bug in clipConvexPolyhedronHP, fix bug in nomalizePlane.m for multiple planes (thanks to Zubiao Xiong), add orientedBox3d |
||

1.19.1.0 | fix bug in polyhedronCentroid |
||

1.19.0.0 | Several updates for processing of meshes (intersectPlaneMesh, intersectLineSphere), for display of meshes, for computation and display of inertia ellipsoids |
||

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

1.17.0.0 | fix various small bugs |
||

1.16.0.0 | fixed bug in surfToMesh |
||

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

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

1.13.0.0 | fix bugs in demos |
||

1.12.0.0 | 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. |
||

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

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

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

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

1.4.0.0 | add missing demos |
||

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

1.0.0.0 |

**Editor's Note:** This file was selected as MATLAB Central Pick of the Week

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.

David Legland (2020). geom3d (https://www.mathworks.com/matlabcentral/fileexchange/24484-geom3d), MATLAB Central File Exchange. Retrieved .

Created with
R2019a

Compatible with any release

**Inspired by:**
iscoplanar.m

**Inspired:**
BrillouinzoneBCCLattice(), Rotate point cloud on orizontal plane - script (ITA, ENG), STL Lattice Generator, intersectPlaneSurf II, Perspective projection, Optimal Step Nonrigid ICP, Geometric measures in 2D/3D images, GJK algorithm distance of closest points in 3D, image ellipsoid 3D, 3D_voronoi_cuboid_bounded

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

- demoDrawLine3d
- demoDrawPlane3d
- demoDrawPoint3d
- demoDrawTubularMesh
- demoGeom3d
- demoInertiaEllipsoid
- demoRevolutionSurface
- demoTransformPoint3d
- drawSoccerBall

- createTrefoilKnot.m
- demoClipMeshVertices.m
- demoConcatenateMeshes.m
- demoCutMeshByPlane.m
- demoPolyhedra
- demoRemoveMeshFaces.m
- demoRemoveMeshVertices.m
- demoSplitMesh.m
- demoTriangulateFaces
- demoVoronoiCell

- anglePoints3d
- angles3d
- angleSort3d
- boundingBox3d
- box3dVolume
- boxes3d
- cart2cyl
- cart2sph2
- cart2sph2d
- circle3dOrigin
- circle3dPoint
- circle3dPosition
- circles3d
- clipConvexPolygon3dHP
- clipEdge3d
- clipLine3d
- clipPoints3d
- clipPolygon3dHP
- composeTransforms3d
- Contents.m
- createBasisTransform3d
- createEdge3d
- createLine3d
- createPlane
- createRotation3dLineAngle
- createRotationOx
- createRotationOy
- createRotationOz
- createRotationVector3d
- createRotationVectorPoint3d
- createScaling3d
- createSphere
- createTranslation3d
- crossProduct3d
- cyl2cart
- cylinderSurfaceArea
- dihedralAngle
- distanceLines3d
- distancePointEdge3d
- distancePointLine3d
- distancePointPlane
- distancePoints3d
- distancePointTriangle3d
- drawAxis3d
- drawAxisCube
- drawBox3d
- drawCircle3d
- drawCircleArc3d
- drawCube
- drawCuboid
- drawCylinder
- drawEdge3d
- drawEllipse3d
- drawEllipseCylinder
- drawEllipsoid
- drawGrid3d
- drawLabels3d
- drawLine3d
- drawPartialPatch
- drawPlane3d
- drawPlatform
- drawPoint3d
- drawPolygon3d
- drawPolyline3d
- drawSphere
- drawSphericalEdge
- drawSphericalPolygon
- drawSphericalTriangle
- drawSurfPatch
- drawTorus
- drawVector3d
- edgeLength3d
- edges3d
- edgeToLine3d
- ellipsoidSurfaceArea
- equivalentEllipsoid
- eulerAnglesToRotation3d
- fillPolygon3d
- fillSphericalPolygon
- fillSphericalTriangle
- fitAffineTransform3d
- fitCircle3d
- fitEllipse3d
- fitLine3d
- fitPlane
- hypot3
- inertiaEllipsoid
- intersectBoxes3d
- intersectEdgePlane
- intersectLineCylinder
- intersectLinePlane
- intersectLinePolygon3d
- intersectLineSphere
- intersectLineTriangle3d
- intersectPlanes
- intersectPlaneSphere
- intersectRayPolygon3d
- intersectThreePlanes
- isBelowPlane
- isCoplanar
- isParallel3d
- isPerpendicular3d
- isPlane
- isPointOnLine3d
- isTransform3d
- linePosition3d
- lines3d
- lineToEdge3d
- medianPlane
- mergeBoxes3d
- midPoint3d
- normalizePlane
- normalizeVector3d
- oblateSurfaceArea
- orientedBox3d
- parallelLine3d
- parallelPlane
- planeNormal
- planePoint
- planePosition
- planes3d
- planesBisector
- points3d
- polygon3dNormalAngle
- polygonArea3d
- polygonCentroid3d
- polygons3d
- projLineOnPlane
- projPointOnLine3d
- projPointOnPlane
- prolateSurfaceArea
- randomAngle3d
- randomPointInBox3d
- recenterTransform3d
- registerPoints3dAffine
- reverseLine3d
- reversePlane
- revolutionSurface
- rotation3dAxisAndAngle
- rotation3dToEulerAngles
- sph2cart2
- sph2cart2d
- spheres
- sphericalAngle
- sphericalVoronoiDomain
- surfaceCurvature
- transformLine3d
- transformPlane3d
- transformPoint3d
- transforms3d
- transformVector3d
- triangleArea3d
- vectorAngle3d
- vectorCross3d
- vectorNorm3d
- vectors3d

- angle3Points
- circleToPolygon
- createLine
- createRotation
- createTranslation
- deg2rad
- distancePointLine
- distancePoints
- findClosestPoint
- isAxisHandle
- isPointInPolygon
- lineAngle
- localToGlobal3d
- minDistancePoints
- normalizeAngle
- orientedBox
- parseDrawInput
- polygonArea
- polygonCentroid
- rad2deg
- rotateVector
- splitPolygons
- splitPolygons3d
- transformPoint
- vectorAngle

- boxToMesh
- checkMeshAdjacentFaces
- clipConvexPolyhedronHP
- clipMeshVertices
- collapseEdgesWithManyFaces
- concatenateMeshes
- Contents.m
- createCube
- createCubeOctahedron
- createDodecahedron
- createDurerPolyhedron
- createIcosahedron
- createMengerSponge
- createOctahedron
- createRhombododecahedron
- createSoccerBall
- createStellatedMesh
- createTetrahedron
- createTetrakaidecahedron
- curveToMesh
- cutMeshByPlane
- cylinderMesh
- distancePointMesh
- drawFaceNormals
- drawMesh
- drawPolyhedron
- ellipsoidMesh
- ensureManifoldMesh
- intersectLineMesh3d
- intersectPlaneMesh
- isManifoldMesh
- isPointInMesh
- mergeCoplanarFaces
- mergeMeshVertices
- meshAdjacencyMatrix
- meshBoundary
- meshBoundaryEdgeIndices
- meshBoundaryVertexIndices
- meshDihedralAngles
- meshEdgeFaces
- meshEdgeLength
- meshEdges
- meshFace
- meshFaceAdjacency
- meshFaceAreas
- meshFaceCentroids
- meshFaceEdges
- meshFaceNormals
- meshFaceNumber
- meshFacePolygons
- meshSurfaceArea
- meshVertexClustering
- meshVertexNormals
- meshVolume
- minConvexHull
- polyhedra
- polyhedronCentroid
- polyhedronMeanBreadth
- polyhedronNormalAngle
- polyhedronSlice
- readMesh_off
- readMesh_ply
- removeDuplicateFaces
- removeInvalidBorderFaces
- removeMeshEars
- removeMeshFaces
- removeMeshVertices
- smoothMesh
- sphereMesh
- splitMesh
- steinerPolytope
- subdivideMesh
- surfToMesh
- tetrahedronVolume
- torusMesh
- transformMesh
- triangulateCurvePair
- triangulateFaces
- triangulatePolygonPair
- trimeshEdgeFaces
- trimeshMeanBreadth
- trimeshSurfaceArea
- trimMesh
- vertexNormal
- writeMesh_off
- writeMesh_ply

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

Select web siteYou can also select a web site from the following list:

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

- América Latina (Español)
- Canada (English)
- United States (English)

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

Nguyen NhungDavid Legland@LC: not for the moment, but this could be easily included I think. I'll add this to the task list!

LCThanks for the great package. I see that there are codes to find the intersection between a 3D line with a 3D polygon. Is there a code to find the intersection between a 3D line SEGMENT and a 3D polygon? Also, a 3D line segment with a 3D polyhedron?

David Legland@Hamish Marsden: unfortunatley not... This is a feature I would like to have as well! The computational complexity may be however very high. Maybe one day...

Hamish MarsdenHi David, great toolbox!! I am trying to merge two intersecting meshes into a single mesh. Are there any functions that could do this, or do you have any advice on how this can be done

David Legland@Xuguang Wang: no, there is no direct function to do this. One possible way to do is to concatenate the two vertex arrays, then relabel the face index array of the second mesh, and concatenate the two face index arrays. This would result in a single mesh, but this will not resolve intersecting faces nor multiple vertices...

Xuguang WangHi David, I would like to merge two meshes into one single mesh. Have you a function to do it? Thanks.

Mohammad AskariAdilAdilThis library is extremely useful and very well documented. 10/10 would recommend.

David Legland@Mert Aslantürk: The function is intended to compute the coordinates of a geometric data structure in a new basis. When applying the transform, the input should be the coordinates of the object in the "source" basis, and the output is the coordinates in the "target" basis. Another example:

plane1 = [0 0 0 1 0 0 0 1 0]; % global basis

plane2 = [2 3 4 0 1 0 1 0 1]; % another basis

transfo = createBasisTransform3d(plane1, plane2);

plane2T = transformPlane3d(plane2, transfo) % results in identity

I agree this is not very intuitive, but I could not find better for the moment...

Ana Guerra LanganMert AslantürkHi David,

Great work! I enjoy using this library. I have one question which is puzzling me about createBasisTransform3d.m

%%%%%%%%%%%%%%%%%%%%%%%%%

source = [0 0 0 1 0 0 0 0 1];

target = [0 0 0 0 1 0 0 0 1];

tForm = createBasisTransform3d(source,target);

out = transformPlane3d(source,tForm)

out =

0 0 0 0 0 1 -1 0 0

I have expected to have out = target. I would really appreciate if you can explain what I am missing here. Thank you for your time.

Best

David LeglandHi John Zhang,

for creating the "faces" elements of the mesh from the vertices, there are different ways, but this requires to add assumptions about which vertice are neighbors. There is no such implementation in geom3d (unless creating the convex hull is enough for you?).

I'll check your proposition for clipConvexPolyhedronHP (and one day, I will implement a cleaner clipMesh function...);

David LeglandHi Ramaswamy Kasyap Pantangi,

there is a function to draw a cylinder (called drawCylinder). There is no "direct" way to draw a semi-cylinder (I am not sure about what you mean by that), but Iguess that checking the code for drawCylinder should help you implement you own function.

Best,

John ZhangHey David,

I believe there may be a fix possible in the function clipConvexPolyhedronHP, where in the last loop, we can add a

face(face == 0) = [];

to make sure that the face we are using does not include any of the 0 index vertices that we used the plane to cut off.

John ZhangHey David, I have a question about how I could create a 3D mesh from an Nx3 array of vertices (sorry if this is an easy question). I tried the steinerPolytope function and the surfToMesh function, but those seem to not fit what I had in mind. The vertices in the array are in no particular order, and I want to return the Faces-Vertices representation so that I can use the other functions in the library.

Thank you for the library! I am relatively new to Matlab, and I find these 3D functions fast and easy to use.

Ramaswamy Kasyap Pantangihow Can I draw a semi cylindrical foil with a nbym grid

cao xiangzhouzhen yejdiva6t9David LeglandHi Tan Phan,

I suspect this is due to some changes that occured to the way Matlab managed Graphical widgets since 2008. What is the exact error message? Anyway, I recommend using the most recent version.

Tan PhanWhy when i using drawCube function on 2008b it's throw error with Patch function but when i run on 2016a it's run perfectly? Thank you and sorry for my bad English!

Xiaofei Hu@David Legland Hello David, hopefully to see your contribution for 3D solid meshing. So far, do you know any other matlab library can do this ?

David Legland@Xiaofei Hu: 3D solid meshes are not managed by the geom3d library. I have made some attempts by considering vertices+edges+faces+solids structures, but did not obtain stable enough result.

Xiaofei HuHello, David, excellent job. I was wondering if it possible to generate a mesh for a 3D object with 3D solid element (tetrahedron for example), just like in a finite element mesh. Could you please show me a way to do that, thanks.

Jordan GrahamThat has worked perfectly David thank you.

David Legland@Jordan: I see; one solution could be to use the second output argument of the "intersectLineMesh3d" function. It corresponds to the position of the intersection point on the straight line. Postion 0 corresponds to line origin. if position is <0, it corresponds to the 180° case. If position is >1, it corresponds to a translation by more than one direction vector. By adjusting the norm of the direction vector of the "line" input argument (params #4 to #6), you should get what you want.

Jordan GrahamThanks for the reply David, the edges function worked well in plotting a finite line, however i have since found that it has not solved my problem. I am attempting to use this finite line to determine intersection points in a mesh. My issue is that the lines are originating from a central point in a circular pattern (Imagine a bicycle wheel laid flat). Say for example the mesh is a box lying beside the point of origin, I want to be able to calculate how many spokes of the wheel intersect the mesh, I can do this, but I also get intersections from the lines originating at 180° as they are not finite in one direction from the starting point. I don't think there is an intersectEdgeMesh3d function.

Jingjing MengDavid LeglandHi Jordan,

I think you can use the functions related to "edges" (the terminology I used in this library for "finite" lines between two points). You can create an edge by concatenating the coordinates of the two points (-> E = [P1 P2]). You ca ndraw the edge with the drawEdge function. Have a look at the file "edges.m" for a more detailed doc.

Jordan GrahamDavid, thank you for the great library. When creating a 3D line is there a way to have this line exist only between two certain points rather than having a seemingly infinite line which goes through the two specified points? for example if i have points P1 and P2 i want the line to exist only between these two points and not extend beyond either point in any direction. I hope this is understandable, happy to provide more details if needed. Thanks again.

Florian ReinboldLionel Juillen AFDavid Legland@Sadeep, no, there is no such function in the library.

Sadeep ThilakarathnaDoes this have any function to check whether two random ellipsoids intersect ? Thank you

David Legland@Sergio ABP, @the neuromechanists: Thanks for your feedbacks, I have updated the main repository accordingly.

@Golak Mahanta: yes, the centroid function is part of the geom2d toolbox. You just need to download it as well. I will include the requested file in the next release to facilitate the use of geom3d.

Christopher Riehsneuromechanist@Golak Mahanta

"centroid" is a Matlab core function, it is introduced in R2017b though. It is also a function in Geom2d. Given that minConvexHull was there long before 2017, I think you just need to import Geom2d toolbox as well.

Matlab's centroid function: https://www.mathworks.com/help/matlab/ref/polyshape.centroid.html

Geom2d's centroid: https://www.mathworks.com/matlabcentral/fileexchange/7844-geom2d

Golak MahantaUndefined function or variable 'centroid'.

Error in minConvexHull (line 61)

pointsCentroid = centroid(points);

Error in createSoccerBall (line 37)

f = minConvexHull(n)';

neuromechanistHi David

First thanks for your excellent code. It already saved me days and maybe weeks of coding.

I was working on your "distancePointMesh", and I think I noticed a subtle inconsitency in the variable names in the linear algorithm. Once I ran the code, it threw this error:

Undefined function or variable 'minDist'.

Error in distancePointMesh (line 124)

if distf < minDist

which makes sense, because minDist is not defined before at all (there is a place on line 113, where "dist" is getting preallocated with inf).

So the fixes I made that worked for me include the following:

1- change "dist" to "dists" on the lines 1, 91, 93 and 104.

2- change "dist" to "minDist" on the line 113.

Thanks and Merry Christmas.

Sergio ABPI added this lines to the file drawCylinder.m right after "%% Display cylinder mesh" to allow transparency:

alpha = 1;

for k=1:length(varargin)

if (strcmp(varargin(k),'FaceAlpha'))

alpha = varargin{k+1};

end

end

And changed this lines:

patch(x2(1,:)', y2(1,:)', z2(1,:)', color, 'edgeColor', 'none','FaceAlpha',alpha);

patch(x2(2,:)', y2(2,:)', z2(2,:)', color, 'edgeColor', 'none','FaceAlpha',alpha);

Nice library!

David Legland@Tao Zhang: thanks for the fixes! I have included them to the master branch of the MatGeom repository (https://github.com/mattools/matGeom.git). They will be included in next release.

Tao Zhangother typos:

"createRotationOx" line51: elseif length(varargin)==3 should be ==4, and line56 should be theta = varargin{4}

"createRotationOy" line55: dy = 0 should be dy=varargin{2}, and line56 should be dz = 0

"createRotationOz" line53: elseif length(varargin)==3 should be ==4, and line58 should be theta = varargin{4};

Tao ZhangHi David,

Just noted some minor typos:

"polygonArea3d" line65: vn = vectorNorm(cp) can be replaced by vn = vectorNorm3d(cp) so that the toolbox is self-contained.

"polygonCentroid3d" line33: elseif nargin == 2 should be elseif nargin == 3

David Legland@Chitan Gandhi: the intersectPlaneMesh function curruntly only works with closed surface meshes. When the mesh is open, some edges are attached to only one face, and one obtains the error you got. Enhancement is on the TODO list, but not that straightforward to manage...

@Laslo Kormoczi: I have submitted a new version packaged as toolbox. This was not the case for the previous version. So this should better work now.

Chintan GandhiHi David,

I was able to use triangulate function and worked fine. Now , I am using intersectplanemesh function on mesh(arranged in the required format from .STL data).Its giving me error : Error using intersectPlaneMesh (line 85)

crossing edge 1 (7,8) is associated to 1 faces. Could you please help me regarding the error.

Laszlo KormocziHi David,

In the current release the toolbox installer does not set the matlabPath to include the toolbox folders. Please review configuration.xml in the mltbx file.

David Legland@Chitan Gandhi: sure! This function was created to manage polygonal meshes whoses faces have many vertices. You can try the following:

[V, E, F] = createDodecahedron;

F2 = triangulateFaces(F);

drawMesh(V, F2);

axis equal; view(3);

Chintan GandhiHello David,

I want to use triangulate face function , is it possible to triangulate face having more than 4 vertices ?

David Legland@Resmi Johnson: there seems to be some limitations on code generation when functions or scripts contains cell functions. See for example https://uk.mathworks.com/help/simulink/ug/cell-array-restrictions-for-code-generation.html. I will try to have a look, but I am not sure if I can solve this easily...

RESMI JOHNSONHello David,

In my project , I am using Embedded Coder for generating C code from matlab .So while generating C code from matlab, I am getting error with respect to 'cell' function for the function 'intersectPlaneMesh(plane, v, f)' . So could you help me to solve this.

David Legland@resmi Johnson: the cell function is part of Matlab for many years, so I do not really understand where does the problem come from... Its aim is to allocate an array of given size (here, 1-by-nFaces), containing data with various types or dimensions. Also, what do you mean by "to generate the code"?

RESMI JOHNSONHello David,

I am using function 'intersectPlaneMesh(plane, v, f)' in one of my simulation. And when I tried to generate the code for this, it is giving me an error because of " cell " function. Could you suggest an alternate solution for this to make it code generatable.

Thanks in advance.

pcidreirNawar AlweshGreat David, thanks for your help and quick reply

David LeglandHi Nawar,

sorry, there was a bug in the orientedBox3d function... As a quick fix, you can replace line 55 by the following:

markers2d = [0 0; L1/2 0; 0 L2/2];

I make a new release soon.

David

Nawar AlweshHi David

Thanks so much for your code.

I was trying to use the orientedBox3d function but I always end up with the pts outside the box, what am I doing wrong? Thanks

eg

X = rand(300,3);

box3d = orientedBox3d(X);

Barbara WorthamDavid Legland@claudio: the function 'isPointInPolygon' is included in the "geom2d" contribution, also on FileExchange. If you downloadit and install it, this will solve the pbm.

I will try to include the missing files in a future release to avoid this kind of bad surprises...

David

David Legland@Sunita: sorry for the long delay... The function "clipConvexPolygon3dHP" does not work for a cube directly, but it is possible to create six planes corresponding to each face of the cube, and apply the function six times. I will try to include a more packaged version of this approach (this was the initial plan, in fact...). It does *not* work for triangulated surfaces (while it should be possible to call it for each triangle, but would be rather inefficient).

You can check the "clipMeshVertices" function -> it works for triangulated surfaces, by specifying a bounding box (you just need to convert the cube to bounding box). However, it does not clip the triangles with the box -> it keeps only vertices located within the box, and the faces attached to them.

claudio carraraHi David, when I use the functions 'intersectRayPolygon3d' and 'intersectLinePolygon3d', I have this error message:

Undefined function or variable 'isPointInPolygon'.

How to solve this?? where is the function ''isPointInPolygon''?

Sunita SahaSunita Saha@David thanks for your suggestion and your code is really helpful.

@David will 'clipConvexPolygon3dHP' function work to clip the area of a triangulated surface inside a cube? Does it give the intersection points of the triangulated surface with the cube i.e. inside the cube?

David Legland@Patrick Hsieh: A first possibility is to combine several rigid transformation matrices: one for translation to plane center, one for rotating first direction vector, and another one for rotation second direction vector. The library provids functions for creating rotation matrices around X, Y, or Z axes. One difficulty would be to compute the rotation angle(s) from the direction vectors of the planes. The atan2 function should help for that.

Also, you can check the "createBasisTransform3d" function. The principle is to compute the (rigid) transform matric that will map points from one plane to another plane. I suspect you could use it in a rather straightforward way, if I udnerstand well your problem.

Patrick HsiehI should rephrase my question. Which sequences of functions should I use in order to transform a triangular plane onto another such that they have one vertex overlapped and the direction of any one of the two edges containing the common vertex aligned. For example, both of them have a vertex on (0,0,0), lie on the x-y plane, and have one edge sitting on X axis?

Patrick HsiehI am doing study on the deformation of vesicles in 3-dimension. Your codes are extremely helpful. Is transformPlane3d the correct choice if I want to transform a trangular surface to the xy plane?

David LeglandHi Sunita,

the plane argument must be 1-by-9 array, and the sphere argument must be 1-by-4. It is possible to process multiple input at a time, but in this case the two inputs must have the same number of rows, or one of the two inputs must have one row.

To compute the area of the circle, you need to extract he fourth column of the CIRC output, and use the area=pi*r^2 formula.

Hope this helps?

David

Sunita SahaAfter getting the circles how to calculate the area of those circles

Sunita SahaError using intersectPlaneSphere (line 46)

data should have same length, or one data should have length 1

Error in Finalsurface (line 299)

circ=intersectPlaneSphere(Plane,S);

How to solve this??

David LeglandHi Jihee,

well, this is again a function that I want to include for a long time... One quick workaround is to use the "intersectLineMesh3d" function. When choosing an orientation, creating the line, computing the intersections, and removing the points with "position" below than 0, one obtain the intersections with a 3D ray emanating from the point. If the number of intersections is odd, then the point is inside the mesh. I have not tested the stability however.

Jihee HanThanks for the quick response.

I have one more question.

What code should I use to check if a point is inside a polyhedron?

Thanks in advance.

David LeglandHi Jihee,

thanks for your feedback!

At the moment, there is no totally generic function for computing the intersection of two polyhedra. As you guessed, this is a complicated problem... However, you can consider the "clipMeshVertices" and "cutMeshByPlane" functions. If one of the mesh is convex, clipping the other mesh with the planes corresponding to each face could be a viable solution.

And... yes, I would like to include such a function in the future!

Jihee HanJihee HanThank you for the sharing.

This file is very useful for me to find any intersection points. Examples in your codes are easy to understand and use the codes.

However, I couldn't find a code that gives the intersection between two polyhedra.

I guess it should be a complicated problem since a polyhedron is defined by multiple faces and vertices, and then there are two of them.

Is there any idea of finding the intersection between two polyhedra by combining the current codes in geom3d/mesh3d?

Or, if possible, could you consider a code that finds the intersection between two polyhedra in a future release?

Thank you again for your help!

David Legland@Faez: for the moment, it is not possible to use multiple polygons with the function "intersectRayPolygon3d". However, I think it is possible to implement a combination of the two functions (intersectRayPolygon3d and TriangleRayIntersection) that avoid the use of repmat and allows multiple output. I will try to investigate this for a future release.

Faez AlkadiHello David,

Is it possible to find intersection of number of rays with different origins with polyhedron (faces and vertices) using the function (intersectRayPolygon3d). Because for what I understand, this function could work for multiple rays but single polygon. But what i'm looking for is multiple rays, multiple polygons (polyhedron).

Jaroslaw Tuszynski has done that in a function called (TriangleRayIntersection). But this function is kind of slow when using repmat.

Thank you so much.

David LeglandHi Katherine,

I do not work with STL files, so it will be complicated for me to answer precisely. The 'intersectLinePolygon3d' takes as second input the list of vertex coordinates of the polygon. If STL is stored as a face-vertex array, maybe there was a problem in generating this list? You can sent me a minimal working example by email and I'll try to have a look.

Katherine BeaulieuHello David,

I am trying to use your function intersectLinePolygon3d to find the point of intersection between a line and an stl object. I am just using the vertices for the variable 'poly' but I don't get an intersection when there should be 2. I was wondering if my input for 'poly' works in this case?

Thank you!

David Legland@David: unfortunately I do not use the PDE toolbox, so I do not know if it is possible or not to create such models using the functions in geom3d...

David NaviauxThanks David for a well thought out set of functions. I have a question regarding the creation of a stack cuboids of various widths and depths that can be used to model a power semiconductor for thermal analysis through the PDE toolbox. Each cuboid needs to have its own thermal characteristics. I know that I can use the matlab "multicuboid" function to create a multi-cell geometry if all of the cuboids have the same width and depths.

Can I use your library to create my semiconductor model?

Thank you!

Abdelrhman HanyHannahThanks David, that is brilliant. I had sorted out a solution for when the cylinders were upright but this has solved the issues I had just come across when I rotated my cylinders.

David Legland@sarvenaz: the 4*PI factors comes from two things. First, mean breadth is proportional to the "integral of mean curvature" with a ratio 2*PI. Then, the mean curvature is concentrated on edges (no curvature within faces, and no contibution of vertices). For each edge, one of the main curvatures is zero, the other one is the dihedral angle between adjacent faces. Hence an additional factor of 2, resulting in 2*2*PI=4*PI. You can also check the refs in "polyhedronMeanBreadth".

sarvenaz babakhaniHi David

Thank you for your practical code

I used "trimeshMeanBreadth" that's awesome

But I have a question about final formula

I wanna know about dividing to "4pi" in eight hand side, what is it for?

David Legland@Miten: I have submitted a new version, that should fix the strange behaviour you mentioned.

@Hannah: I have added detection of intersection between line and cylinder end-cap. So this should work better for you now!

Also, many functions have been updated and several new ones introduced. Check out GitHub repo for details (https://github.com/mattools/matGeom)

David Legland@Hannah: indeed, the cylinder is considered to be "open". If you need the intersection with the extremity, one workaround is to compute intersection of the line with the 3D plane containing the end of cylinder (in practice, you would need to consider 2 planes, one for each extremity). Then you can test if the distance between the intersection point and the cylinder main axis is less than the cylinder radius. I will try adding this behaviour in a next release.

HannahHi David,

I'm trying to calculate path lengths through a cylinder using the intersectLineCylinder. It is working very well until the line goes in the side of the cylinder but is angled enough to come out at the top. The output is only one set of coordinates. Is this because the top of the cylinder is not recognised as a surface? If this is the case, do you have any suggestions how I could work around this?

Many thanks, Hannah

Patrick BolanExcellent! I've always gone back to VTK for 3d computational geometry (comparing MR spectroscopy with oblique MR imaging planes) but this toolkit makes it easy. Very clear, concise.

Miten PatelHi David,

I'm trying to understand how to use geom3d library to translate a series of co-planar points onto a plane defined by a different series of co-planar points.

For example, if I have a series of points that define a unit circle on the YZ plane centred on 0, 0, 0 and I had a plane at a random 3D angle and random origin then I would like the 3D global coordinates that define the same unit circle on that plane.

As a test of my understanding, I set up a matrix 'z' of following 4 points z = [ 0 1 0; 0 0 -1; 0 -1 0; 0 0 1]. These are on the global plane. I want to translate them to the plane defined as [ 0.0544 0.0000 0.0000 0.0000 0.9959 -0.0901 0.0000 -0.0901 -0.9959 ]. This is simplified test where destination plane is in the YZ plane and just translated along the X axis.

I used the following statements:

z = [ 0 1 0; 0 0 -1; 0 -1 0; 0 0 1];

lp = [ 0.0544 0.0000 0.0000 0.0000 0.9959 -0.0901 0.0000 -0.0901 -0.9959 ];

t = createBasisTransform3d('global', lp)

transformPoint3d(z, t) gives

0.9959 -0.0901 0.0544

0.0901 0.9959 0.0544

-0.9959 0.0901 0.0544

-0.0901 -0.9959 0.0544

This is not what I expected to get. Firstly the destination plane is translated by 0.0544 on X axis whereas the new points have this translation reflected on the Z' coordinate.

Also, the resulting X', Y' coordinates are what I thought should be in the Y', Z' coordinates. Also, the rotation around the four points is counter-clockwise (if X' > Y and Y' > Z) whereas the original points were in a clockwise rotation.

Can you advise where I am going wrong?

Many Thanks, Miten

Halil Ibrahim KayimHi David,

Thanks for this library, it's helping me a lot in my master's study. The only trouble I had with this awesome toolbox was with the mergeCoplanarFaces method. It seems you added the threshold variable in a later version, however I don't think it works as it should. After investigating how this threshold was used, I saw that it is used for two purposes. First one is to test the value of a cross product of two normal vectors. The second one is later used inside mergeCoplanarFaces when the candidate nodes are tested for coplanarity in the isCoplanar method. The trouble with this is that while the first threshold is only meaningful for values between 0 and 1, the second threshold can be any value greater than 0, and does not behave the same for geometries of different scale.

Now, I don't know how to go about solving this behavior, but it doesn't seem to break anything right now to just comment out the isCoplanar check inside mergeCoplanarFaces. This may be a horrible workaround and may ruin stuff in the future, so can you please have a look at this issue?

Thanks for this toolbox and for your troubles.

Best,

Halil

David Legland@ThT: I do not know why this behaviour is happening. Please send me a sample data set by email, and I will try to have a look.

@Leonardo: If I understand well you would like to compute intersection between a 3D triangle and a 3D polyline? Unfortunately there is no such function for this at the moment. I think compmuting intersection with line, and post-processing the result could give some results, but this would require some coding efforts.

Regards,

David

Leonardo ColavittiHi David,

thanks for the library and your efforts. I'm using the function TriangleRayIntersection of the package geom3D and I was wondering if in the input, instead of giving the direction of a straight line, it exists an implementation if my line is a slightly curve line and I know it in every point (basically, I have a vector in the 3 dimensions). I need this for a ray-tracing problem...

Thanks for your kindly help,

Best regards,

Leonardo

ThTHi David,

Thanks for the library and your efforts. I have a question though regarding an issue that I am facing. I am loading a CAD model in matlab from an .stl file. This gives two matrices with the vertices and the faces of the model, already in triangles. However, the existing trianglulation is a kind of weird and for that reason I am trying to reformat it in a better way. Therefore, I am using the "mergeCoplanarFaces()" in order to create faces with more that 3 vertices and then re-triangulate the new mesh in a more proper way. However, the output of the "mergeCoplanarFaces()" function gives me some multi-vertex faces with "NaN" points included and if I try to plot the new faces and vertices there are some gaps within the mesh, which were not existing in first place. Do you know what it might be happening, and if there is any workaround in order to avoid this behaviour.

Thanks.

Theo

David LeglandHi Lina,

The (dx,dy,dz) components of the line representation simply correspond to the differences in coordinates of the two 3D points.You can check the utility function "createLine3d", that encapsulates this process.

You can try the following example:

sphere = [40 50 60 30];

figure; drawSphere(sphere);

hold on; axis equal;axis([0 100 0 100 0 100]);

p1 = [50 50 50]; p2 = [51 52 53];

line = createLine3d(p1, p2); drawLine3d(line, 'b');

pts = intersectLineSphere(line, sphere);

drawPoint3d(pts, 'k*')

Otherwise, you can send me your data by email, I,can try to have a look;

Regards,

David

Tsvetelina GermanovaHi David,

I'm using the Intersect line sphere function but I seem to not obtain the proper angle of the intersection. I think I may be giving a wrong input for the dx dy dz component of the line (or at least I cannot find other explanation). Can you please give me an example of what would be the the dx dy and dz if I have a line defined by two points?

Thanks,

Lina

David LeglandHI Annie,

I would suggest trying some functions from the geom2d package (also available on FEx):

* the "removeMultipleVertices" function merge adjacent vertices whose distance is small enough

* the "mergeClosePoints" merge points with distance small enough, without taking into consideration point order (ie, it is possible to merge point with not adjacent indices)

These two functions were designed for planar geometry, but should for for 3D as well.

Regards,

David

Annie TranHi David,

I used the intersectPlaneMesh and when calculating the intersection between a plane that coplanar with faces of the mesh, the output polygon contained repeated vertices. Could you suggest idea to eliminate the repeated vertices since it is helpful for me in further process.

HaT

Olabanji ShonibareSuzie OmanThe lines, sphere intersect isn't very useful. The example is not easy to follow.

Li ZhengDavid LeglandHI Alec,

you can check the sphericalAngle function, that computes angle between three points on a sphere. Maybe you will need to convert spherical coordinates.

If have used the principle to compute spherical angle of Voronoi domains on the sphere for another contribution (estimation of surface area in 3D images). I think you can grab some pieces of code from the source:

https://github.com/dlegland/matImage/blob/master/matImage/imMinkowski/private/sphericalVoronoiDomainArea.m

Regards!

Alec DayHello, I'm trying to use this software (its great by the way) to calculate the area of a polygon stretched over a sphere. If I have a range of azimuth and elevation values that connect to form a polygon over a sphere, which function could i use to calculate its surface area (on say a unit sphere, or better yet, a sphere of radius r)?

Thanks!

Alec.

David Legland@Charles,

sorry for delay... Thanks for your suggestion. I do not use Doxygen, so maybe this will be difficult for me to apply this solution. But this should open ideas for scripting or automatising. Thanks!

David Legland@Chirag,

I do not know exactly where your problem can come from. Could you identify some specificities of the meshes causing troubles? Can you sent me by email a problematic mesh?

Some ways of investigation:

* maybe use the trimMesh function before computing may help.

* appy the triangulateFaces before computing smoothing

Regards,

David

David LeglandHi S SS, and sorry for the delay...

The intersectPlaneMesh function works for the moment only with infinite planes. One solution would be to project the intersection on the plane, to apply 2D intersection algorithms on the projection, and to project back to 3D. I am not sure this is quite simple however...

David LeglandHi Pauline,

the intersectPlaneMesh function requires for the moment a "closed" mesh, that is a mesh without free boudary. Torus and cube mesh fullfil the requirement, but if you apply on open meshes the error you obtained is expected. I will try to work on an enhanced version that manages also the case of open meshes.

Chirag BhuvaHi David ,

I am using smoothMesh function,it is working good but sometimes it is not worked,the error is Inner matrix dimensions must agree.

Error in smoothMesh (line 55)

v2 = adj * v2;

Error in DemoMain (line 31)

[FV2.vertices, FV2.faces] = smoothMesh(FV.vertices, double(FV.faces));

because of adjacency matrix and vertices's size is not same .so could you please guide me.

Pauline HuetThank you David for the code !

I am quite new to matlab/geom3D and I have some trouble with the intersectPlaneMesh function. I am trying to get the intersection of a mesh read from an stl file as faces and vertices (successfully displayed with drawMesh) with a plane. Unfortunately I always get the same error 'crossing edge 1 (4,5) is associated to 1 faces' whereas with standard privitives like cubes or torus all goes right. Any idea of what could have caused this error ? Thanks in advance

S SSDear David, thank you very much for your submission. It is very useful for me.

I am using the intersect plane mesh function to have the intersections between planes and a mesh. Actually I have several intersections because the plane is infinite. I want to have only one intersection. For this I want to have a limited plane. Have you a solution for me? Thank you very much

Charles BrillonDear David, I sent you a personal e-mail (in french). Have you considered using this : https://www.mathworks.com/matlabcentral/fileexchange/25925-using-doxygen-with-matlab to generate documentation?

That might prove helpful.

Regards,

-Charles

David LeglandHi Gary,

yes, keeping consistent numbering of version and doc is complicated... But I agree, this kinf of information is important. I'll try to enhance this, maybe this can be automated in some way.

regards,

David

Gary KenwardDavid:

NP. Glad I could help.

One additional note: you might want to update the version and date in the Contents file. I understand what a pain maintaining documentation can be, however, I just spent a half hour trying to figure out why I had a 2011 version of GEOM3D loaded into Matlab. It can be troublesome to suspect that you are using a very old version of a toolbox.

Also, you might want to update the version to 1.21 or such.

Anyway, great toolbox. Thanks for sharing it with us.

Gary

David LeglandHi Gary,

thank you for reporting the problem, and quickly providing the solution! I have updated a new version that uses your fix (also in other similar functions).

regards,

David

Gary KenwardI fixed the problem I encountered with fillPolygon3d with two changes, marked with the characters " %<======= " below. I may have missed something, but the changes appear to work.

function varargout = fillPolygon3d(varargin)

%FILLPOLYGON3D Fill a 3D polygon specified by a list of points

%

% fillPolygon3d(COORD, COLOR)

% packs coordinates in a single [N*3] array.

% COORD can also be a cell array of polygon, in this case each polygon is

% drawn using the same color.

%

% fillPolygon3d(PX, PY, PZ, COLOR)

% specify coordinates in separate arrays.

%

% fillPolygon3d(..., PARAM, VALUE)

% allow to specify some drawing parameter/value pairs as for the plot

% function.

%

% H = fillPolygon3d(...) also return a handle to the list of line objects.

%

% See Also:

% polygons3d, drawPolygon, drawPolyline3d

%

% ------

% Author: David Legland

% e-mail: david.legland@nantes.inra.fr

% Created: 2007-01-05

% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

% check case we want to draw several curves, stored in a cell array

var = varargin{1};

if iscell(var)

hold on;

h = [];

for i = 1:length(var(:))

h = [h; fillPolygon3d(var{i}, varargin{2:end})]; %#ok<AGROW>

end

if nargout>0

varargout{1}=h;

end

return;

end

% extract curve coordinate

if size(var, 1)==1 %<======= check if var is 1 x N

% first argument contains x coord, second argument contains y coord

% and third one the z coord

px = var;

if length(varargin)<3

error('Wrong number of arguments in fillPolygon3d');

end

py = varargin{2};

pz = varargin{3};

varargin = varargin(4:end);

else

% first argument contains both coordinate

px = var(:, 1);

py = var(:, 2);

pz = var(:, 3);

varargin = varargin(2:end);

end

% extract color information

if isempty(varargin)

color = 'c';

else

color = varargin{1};

varargin = varargin(2:end);

end

% ensure polygon is closed

%<======= fill3 closes polygons if needed; code below generates concatenation dimension inconsistency error

% px = [px; px(1)];

% py = [py; py(1)];

% pz = [pz; pz(1)];

% fill the polygon

h = fill3(px, py, pz, color, varargin{:});

if nargout>0

varargout{1}=h;

end

Gary KenwardI'm having a problem with fillPolygon3d(PX, PY, PZ, COLOR), where PZ, PY, and PZ specify coordinates in separate arrays.

fillPolygon3d([0.4 0.6 0.5], [0.5 0.5 0.5], [0.2 0.2 0.4], 'blue');

Error using fill3

String argument is an unknown option.

Error in fillPolygon3d (line 76)

h = fill3(px, py, pz, color, varargin{:});

Any suggestions?

BinuDavid LeglandHi Dimitris,

unfortunatley there is for the moment no function for computing offset of 3D polygons or meshes. Depending on what you want to do, you can (1) project the 3D polygon on a plane, compute offset in 2D, then project the result back in 3D, or (2) compute an extrusion of the polygon, by duplicating the polygon, translating one of the two polygons by a given distance in the direction perpendicular to the supporting plane of the polygon, then associate the vertices to form a 3D mesh.

Dimitris MylonasHello how to offset polygons in 3d?

Manar Al AsadHello,

I'm experiencing problems with the centroid-related functions (calling non-existing functions). Has anyone else encountered problems with them. If so, does anyone have a fix?

anusha gorrilaCaleb BrownMANIK BANSALCollin StillmanWei ZhaoDavid LeglandHi Anusha,

it is difficult to figure out what could be wrong. Usually, for 'closed' meshes (such as boundary representations), intersesctions with lines produce even number of intersections. The "pos" output corresponds to the relative position on the line. Negative value means the position is before the line origin. You can also investigate the "ind" output, that corresponds to index of intersecting face.

It is still possible there are some numerical computation errors, that cause apparition of unexpected points, but I did not encounter yet...

anusha gorrilaHello David, Thanks for such a wonderful code. I've been using the intersectLineMesh3d function for finding the intersection point on the mesh from a line and i get 2 intersection points and the position of the 1st one is a negative value and the position of 2nd one is a positive value. Could you explain the output as my line intersects the mesh only at one point. The output is something as below:

points =

536.1685 361.9645 40.1591

573.2433 399.9937 77.9874

pos =

-28.0442

10.0410

ind =

19

201

Thanks,

Anusha.

David LeglandHi Hannah,

hmm, strange behaviour I agree... In the current version, the algorithm starts by triangulating the mesh. By consequent, each square face of the cube is transformed into two triangular faces. It the origin of the line is in the middle of the cube, it will intersect the two triangles corresponding to a single square.

By the way, if you use additional ouptput arguments, they will refer to the triangulated mesh! (I will add warning or checking in future version).

regards,

David

HannahThanks for this code, I'm finding it very useful. I just have one question, when I'm using intersectLineMesh3d, why do the intersection points repeat themselves, for example for a cube with a straight line through the middle you get 4 intersection points?

Many thanks, Hannah

David LeglandHi Sandy,

I'm glad you could find the way of using isPointInPolygon!

For other kind of intersections, you can consider 'intersectLineMesh3d', as well as 'intersectPlaneMesh' functions.

regards!

Grant SandyPlease disregard comment on August 5. I see now that the geom2D package is required. So far everything is working perfectly with the addition of geom2D functions.

Grant SandyThere seems to be a missing function called isPointInPolygon which is called by intersectLinePolygon3d. Is this something that was overlooked?

Grant SandyQuite a treasure! I am using the line/plane intersection function now. I will also use the other line/surface intersection functions.

David LeglandHi Tariq,

yes, the inertiaEllipsoid function computes an equivalent ellipsoid from data, based on inertia of the input points.

You can get the centroids either from the first three columns of the result of inertiaEllipsoid, or by using the 'centroid' function from the geom2d package.

For computing surface area of a bounding box, I would suggest summing the area of each face. The area of each face can be obtained by multiplying the corresponding side lengths.

regards

tariq bdairHi David,

Is the inertiaEllipsoid function calculte the best fit ellipse like in 2d data?

Also, how can i get the centroids for double data?

How can i calculate the surface area for bounding box?

Thank you

David LeglandHi Md Shadir,

You will need a function to parse the PLY file, and extract the vertex coordinates and vertex indices of each face. There are several tools in other libs, for example the graph toolbox of Gabriel Peyre (function 'read_ply'), or the computer vision toolbox (function 'pcread').

Note however that to my knowledge, PLY files provide only vertex coordinates. So you will need to compute the faces of your mesh. In this case, I would recommend using delaunay triangulation, and then try to remove non necessary tetraedra (this can be a complicated task however...), and extract their outer face.

regards,

Md Shakir KhanHi David,

great submission. I have been looking into the "geom3d" and trying to implement the "intersectLineMesh3d" for my project. But I need to work with 3D mesh data (.ply file) and check for the intersection. The .ply (3D mesh) file contains element vertex as property floats x,y,z and element faces. Is there any way to feed the .ply file into your code and run the "intersectLineMesh3d" algorithm? In fact, I tried by directly copying the vertices and faces from the .ply file (as below) but got the error: Subscript indices must either be real positive integers or logicals.

vertices = [-0.289329 -0.230469 -0.550781

-0.285156 -0.230469 -0.543497

-0.285156 -0.229169 -0.550781

-0.277344 -0.230469 -0.549533

-0.285156 -0.229169 -0.550781

-0.285156 -0.230469 -0.543497

-0.277344 -0.230214 -0.550781

-0.285156 -0.229169 -0.550781

-0.277344 -0.230469 -0.549533

-0.276716 -0.230469 -0.550781]

However, it works well with the positive data.

Can you please help me out? I would really be very grateful to you. Thanks.

An example (part) of my data:

ply

format ascii 1.0

comment file created by Microsoft Kinect Fusion

element vertex 2959173

property float x

property float y

property float z

element face 986391

property list uchar int vertex_index

end_header

0.300781 0.625023 -1.113281

0.296898 0.628906 -1.113281

0.300781 0.628906 -1.117165

-0.695368 -0.660156 -1.496094

-0.691406 -0.659174 -1.488281

-0.691406 -0.660156 -1.491550

-0.695368 -0.660156 -1.496094

-0.699219 -0.658886 -1.496094

-0.691406 -0.659174 -1.488281

..........................................................

..........................................................

..........................................................

3 2959146 2959147 2959148

3 2959149 2959150 2959151

3 2959152 2959153 2959154

3 2959155 2959156 2959157

3 2959158 2959159 2959160

3 2959161 2959162 2959163

3 2959164 2959165 2959166

3 2959167 2959168 2959169

3 2959170 2959171 2959172

Md Shakir KhanBut I need to work with 3D mesh data (.ply file) and check for the intersection. The .ply (3D mesh) file contains element vertex as property floats x,y,z and element faces. Is there any way to feed the .ply file into your code and run the "intersectLineMesh3d" algorithm? In fact, I tried by directly copying the vertices and faces from the .ply file (as below) but got the error: Subscript indices must either be real positive integers or logicals.

vertices = [-0.289329 -0.230469 -0.550781

-0.285156 -0.230469 -0.543497

-0.285156 -0.229169 -0.550781

-0.277344 -0.230469 -0.549533

-0.285156 -0.229169 -0.550781

-0.285156 -0.230469 -0.543497

-0.277344 -0.230214 -0.550781

-0.285156 -0.229169 -0.550781

-0.277344 -0.230469 -0.549533

-0.276716 -0.230469 -0.550781]

However, it works well with the positive data.

Can you please help me out? I would really be very grateful to you. Thanks.

An example (part) of my data:

ply

format ascii 1.0

comment file created by Microsoft Kinect Fusion

element vertex 2959173

property float x

property float y

property float z

element face 986391

property list uchar int vertex_index

end_header

0.300781 0.625023 -1.113281

0.296898 0.628906 -1.113281

0.300781 0.628906 -1.117165

-0.695368 -0.660156 -1.496094

-0.691406 -0.659174 -1.488281

-0.691406 -0.660156 -1.491550

-0.695368 -0.660156 -1.496094

-0.699219 -0.658886 -1.496094

-0.691406 -0.659174 -1.488281

..........................................................

..........................................................

..........................................................

3 2959146 2959147 2959148

3 2959149 2959150 2959151

3 2959152 2959153 2959154

3 2959155 2959156 2959157

3 2959158 2959159 2959160

3 2959161 2959162 2959163

3 2959164 2959165 2959166

3 2959167 2959168 2959169

3 2959170 2959171 2959172

Sleh Eddine BrikaNiels SchroeterDavid LeglandHi Andy,

I am not sur what you want to do exactly; one possiblity is to identify faces that are located within a polygonal mesh. In this case, the you will need to detect if all vertices are within the mesh. There is for the moment no such function in geom3d, but I think there is on FEx. Another possibility is to display faces with different transparencies. In that case I think the simpler is to manipulate direclty the graphical handles of the meshes.

Hope this helps ?

Andreas HornDear David, this is a really awesome library!

However, I'm specifically looking for a function that removes internal faces in a complex patch (exactly the functionality as described here: http://meshlabstuff.blogspot.com/2009/04/how-to-remove-internal-faces-with.html). Do you have any idea how I could do this within Matlab?

Thanks so much for any hint!

Andy

dingkun liandingkun liandingkun lianDavid LeglandHi Weihua,

Extracting the position of 3D pixels in a 3D image is not complicated (usign for example the find function), for managing self occlusion there is no straightforward solution in geom3d... One possible solution would be to reconstruct the set of 3D position (using for example marching cubes), and then consider intersections of reconstructed structuers with rays parallel to the direction of the view. By sorting intersections with respect to relative position on ray, a virtual view of the structrues could be obtained. This could be rather computer-intensive however.

UTAHi David

Thank you for your code. I also interesting to extract the pixels positions of the 3D object which in a specific view direction. That is to say, I want to get the 3D pixels locations after consider the self-Occlusion. Would you please give me some help on this?

Thanks ahead!

weihua

David LeglandHi Alex,

unfortunately there is for the moment no such functionality if the geom3D toolbox. But you may check the following (very recent !) submission: http://www.mathworks.com/matlabcentral/fileexchange/55803-surfaces-intersect

If you transform a 3D polygon into a set of triangular faces, I suppose you can use "surface interesect" for computing intersection points.

Regards,

Alex LifshitzHi David,

Thank you for the explanation and correcting the code.

I have another question. Is it possible somehow to clip a 3d polygon by another shape, for example a cylinder rather than just a plane?

Thanks

Alex

David LeglandHi Alex,

The computation of the volume of a mesh is correct: you need to consider signed volume to take into account non convex meshes. You can try with a torus for example... As a quick workaround, you can consider using abs(det), but the function will return correct results only for convex (or star-shaped with respect to origin, I guess) meshes.

Actually, the bug was in the "clipConvexPolyhedronHP" function. A new face is created during the process, and its orientation was not kept consistent with the rest of the mesh. I have updated the package, that should solve the pbm.

regards,

Alex LifshitzFollowing my previous comment below it seems that there is a bug in the implementation of meshVolume function.

The function divides the volume into tetrahedra. Then the total volume equals to the sum of the volumes of all tetrahedra. The code computes the volume by calculating the determinant of the vertices.

vols(i) = det(tetra) / 6;

It turns out that sometimes this determinant can be negative due to a different order of the rows (vertices coordinates).

The correct formula should include absolute value of the determinant.

vols(i) = abs(det(tetra)) / 6;

That solves the problem.

Alex LifshitzHi

I am interested to calculate a volume of the ellipsoid after slicing it with an arbitrary plane. I am using the following code (here for simplicity I will use a sphere and cut it in the center).

[x,y,z] = sphere(30);

[f,v,c]=surf2patch(x,y,z,z,'triangles');

meshVolume(v, f)

This calculates a volume of a sphere as expected (4/3*pi)

Now, I slice it in half.

P = createPlane([0,0,0], [0, 0, 1]);

[v1 f1] = clipConvexPolyhedronHP(v, f, P);

Now if I calculate the volume, I get incorrect value. If I draw it, I see that the reason for this is that the surface along which the sphere was cut, remains open.

drawMesh(v1,f1,'facealpha',0.5)

My question is how can I close the shape such that the volume can be calculated correctly.

Thanks!

David LeglandHi Michael,

sorry, there was a bug in the polyhedronCentroid function... I have just submitted a revised version that should solve this.

David LeglandHi Umair,

to compute the volume, you need to specify how the vertices are connected. There is no "simple" method, but you can search for alpha shapes for example.

Umair Bin AsimIs there any way here, with which we can find the volume of an arbitrary 3D shape defined by a set of 3D points?

Michael StrohmeierHi David,

thank you for your answer.

I have to bother you with another issue.

Wenn i compile polyhedronCentroid.m i get the error 'Undefined function or variable "vt".' Where do i get the data for "vt"?

Thank you in advance,

Michael

David LeglandHi Michael,

unfortunately I think this will not be possible... There is already a vectorization procedure for managing all faces of the mesh in a global way. Then it seems difficult to process also lines in a vectorized way.

Tristan ChaplinHi David,

Thank you very much, it works just as I'd hoped.

Thanks,

Tristan

Michael StrohmeierHi David,

your library is very useful! Thank you for sharing.

I am using intersectLineMesh3d and have to loop the function for several lines in order to get all intersection points.

I was wondering if it is possible to modify the code in order to use multiple lines...like you did for intersectLineSphere recently (last comment of Tristan Chaplin).

This would be very helpful since a loop around the function is too slow for my purpose.

Thank you very much!

Michael

David LeglandHi Tristan,

Yes this is possible! I just updated a new version, with a modified version of the intersectLineSphere function. There is an example in the header of the function that demonstrates how to use it with multiple lines.

regards!

David

Tristan ChaplinHi David,

Thanks for this library, it has been very useful for me.

I was wondering, is there a away to use intersectLineSphere with multiple lines? Similar to the way intersectLinePlane works? It's too slow for my purposes to use a for loop and do each line one by one.

Thanks,

Tristan

David LeglandHi Dale,

in fact, both 'line' and 'ray' are infinite. If P0 is the origin point and V is the direction vector, both shapes represent a set of points satisfying P = t * V + P0. In the case of a line, t ranges from -inf to +inf, whereas in the case of a ray, t ranges from 0 to +inf. Note also that the LINE or RAY variables represent only origin point and direction vector, and not the 'type' of the shape. The distinction between line or ray is made only by calling the appropriate function.

If you want to find intersection with a finite line (a "line segment"), you can get the position 't' of the intersection point on the line from the 'linePosition3d' function. If value is between 0 and 1, the points lies on the line segment comprised between P0 and P0+V.

hope this helps,

David

Dale FriedHi David,

Thanks for sharing this great library. I have used it quite a bit already, and am now trying a new application.

I am missing what is the functional difference between intersectLinePolygon3d and intersectRayPolygon3d.

I would expect that the line function would understand a finite length, but the documentation just mentions an origin and a direction vector. Should that read, instead, "direction vector giving the length of the line segment"?

What am I misunderstanding?

Thanks.

-Dale

Krisztian Szuchermuna majeedHi David,

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

please, I ask you can I crop or segment a mesh of type obj into part by this tools, please I need that

David LeglandHi Yiyu,

actually I modified this function yesterday also... I have used another approach, that makes it much faster, and it now returns multiple polygons as well. You can try the new version here: https://github.com/dlegland/matGeom/blob/master/matGeom/meshes3d/intersectPlaneMesh.m . Maybe you can compare with yours? If you have any advice/hint i'm still interested!

regards!

Yiyu HongHi David,

I modified your intersectPlaneMesh function code and make it work to find multiple contours.Do you need the code?

David Legland@Homa Rasouli:

for computing surface area of a polygonal mesh, you need an array with vertex coordinates, and an array that stores vertex indices of each face. Such arrays can be obtained from the 'isosurface' function, for example. Otherwise, you can read them from mesh file formats such as OFF or PLY.

Homa RasouliI have problem in computing the surface of triangular mesh. The function works only with triangle mesh as an input. How should I compute the vertices?

David Legland@SM:

you're right, I use XYZ convention: the three angles correspond to extrinsic rotations performed around X, Y and Z axes, in that order. Note that the result row vector is ordered as [thetaZ thetaY thetaX].

See also the "eulerAnglesToRotation3d" function for details on the reverse operation.

regards, David

SMHi David,

What convention are you following in the Rotation3dToEulerAngles file? I am not completely familiar with the equations, but it looks to me like this is an XYZ convention.

Thanks, Sem

David LeglandHi Johannes,

the volume of polyhedron is defined with respect to the orientation of faces. If faces point outwards, volume is positive, otherwise it is negative. You can simply consider the absolute value of the result. Be careful to ensure that orientation of all faces is consistent: you can use "checkMeshAdacentFaces" for that purpose.

David LeglandHi Dexter,

for drawing a 3D polyhedron, you can use the "drawMesh" function. It uses a classical Face-Vertices representation. For creating the polyhedron, you may have a look at the 'steinerPolytope' function, it seems it may fit your needs.

regards,

David LeglandHi Dexter,

for drawing a 3D polyhedron, you can use the "drawMesh" function. It uses

Johannes NeumannHi David,

thanks for your quick reply. I am using the following code now:

K = minConvexHull(tmpVerts,myEPS)';

tmpInd = unique(cell2list(K));

[K, vol] = minConvexHull(tmpVerts(tmpInd,:),myEPS);

I modified your minConvexHull to pass the volume from qhull and cell2list is cell2mat essentially.

Concerning the volume: I guess it should be

abs(det(tetra) / 6)

in your meshVolume function? Otherwise, I do not get positive volumes in general.

DexterHej, 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!

matteoHi David,

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

David LeglandHi 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.

Matthew BrandsemaHello, 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???

FritzDavidDavid 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

David LeglandHi 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

JinHi 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.

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

AbdulrahmanDavid

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

Great submission THANK YOU

Johannes NeumannDear 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

David LeglandHi 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

Chien Ting CHINI'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!

Chien Ting CHINSorry I pasted the original line. The correct line might be:

elseif nargin >= 6;

JdeenDavid 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 ?

John AndersonHi 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,

David LeglandHi 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!

Germano GomesGreat 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.

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

ShaoshuaiVery helpful! Thanks for sharing!

J.R.! MenzingerVery useful.

David LeglandHi 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

CedricHi 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, :));

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

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

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

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?

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...

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.

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.

Sven@Evan,

Please check out http://www.mathworks.com/matlabcentral/fileexchange/43013-unifymeshnormals

It should do what you're looking for.

Sven@Evan,

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

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 ?

EvanHello 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.

Yong Min YoonHi, 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.

David LeglandHi 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

BenjaminDear 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

David LeglandHi 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

Khaled KhairyWith 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.

David LeglandHi 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

OrestisDear 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

David LeglandDear 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

Dmitri KamenetskyGreat 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.

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

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.

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.

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].

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

panagiotisSven,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!

SvenHi 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 :)

panagiotis@ Sven

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

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.

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

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.

panagiotisThe 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

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.

panagiotisDear 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

David LeglandHi 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

CharlesDear 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

Student UniHI 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

Hoai-Nam LEMany thanks David!!!

Nam

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

YuriDavid, thanks a lot! Yuri

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

Hoai-Nam LEDear 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

YuriDavid, 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

marcVery nice. thanks

APGreat tools.

Henri SHIMylesHaochen TangHi David,

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

Best,

Haochen

EduardoHi 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

David LeglandHi 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

David LeglandHi 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

Haochen TangHi 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

EduardoHi 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

David LeglandHi 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

AmandaHi, 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?

AmandaTianjin UnivDavid LeglandHi 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.

NirHey 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?

HannesThanks mate. Beer flies out to you :)

David LeglandHi 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

Benjamin SanderseHi 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

David LeglandHi Patrik,

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

regards

Patrik EschledistanceLine3d 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?

Patrik EschleNicely written, good comments - thank you!

Shuhao CaoSakshiDavid 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

Qais AlKhazrajiDear 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.

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

Qais AlKhazrajiDear 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

David LeglandHi 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.

KhaledDear 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.

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

Siyi Dengvery helpful.

SvenExcellent 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.

David LeglandThe 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

Paklah Abd RahmanHi,

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

David LeglandHi,

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

Jon SporringThis 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

Nathan ThernThis 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.

usit 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