diff --git a/matGeom/meshes3d/intersectPlaneMesh.m b/matGeom/meshes3d/intersectPlaneMesh.m
index 051eadff..c410d460 100644
--- a/matGeom/meshes3d/intersectPlaneMesh.m
+++ b/matGeom/meshes3d/intersectPlaneMesh.m
@@ -1,15 +1,20 @@
-function [rings, openPolys] = intersectPlaneMesh(plane, v, f)
+function [polys, closedFlag] = intersectPlaneMesh(plane, v, f)
 %INTERSECTPLANEMESH Compute the polylines resulting from plane-mesh intersection.
 %
-%   RINGS = intersectPlaneMesh(P, V, F)
-%   [RINGS, CURVES] = intersectPlaneMesh(P, V, F)
+%   PLOYS = intersectPlaneMesh(P, V, F)
+%   [POLYS, CLOSED] = intersectPlaneMesh(P, V, F)
 %   Computes the intersection between a plane and a mesh. 
 %   The plane P is given as:
 %   P = [X0 Y0 Z0  DX1 DY1 DZ1  DX2 DY2 DZ2]
 %   The mesh is given as numeric array V of vertex coordinates and an array
 %   of (triangular) face vertex indices.
-%   The first output is a cell array of closed polylines, the second output
-%   is a cell array of open polylines ("line strings"). 
+%   The output POLYS is a cell array of polylines, where ech cell contains
+%   a N-by-3 numeric array of coordinates. The (optional) output CLOSED is
+%   a logical array the same size as the POLYS, indicating whether the
+%   corresponding polylines are closed (true), or open (false). Use the
+%   functions 'drawPolygon3d' to display closed polylines, and
+%   'drawPolyline3d' to display open polylines.
+%
 %
 %   Example
 %     % Intersect a cube by a plane
@@ -44,7 +49,7 @@
 %     drawPolygon3d(polySet, 'lineWidth', 2, 'color', 'y')
 %
 %     % Demonstrate ability to draw open mesh intersections
-%     poly = circleArcToPolyline([10 0 5 0 180], 33);
+%     poly = circleArcToPolyline([10 0 5 90 180], 33);
 %     [x, y, z] = revolutionSurface(poly, linspace(-pi, pi, 65));
 %     [v, f] = surfToMesh(x, y, z);
 %     f = triangulateFaces(f);
@@ -53,8 +58,8 @@
 %     drawMesh(v, f, 'linestyle', 'none', 'facecolor', [0.0 0.8 0.0], 'faceAlpha', 0.7);
 %     drawPlane3d(plane, 'facecolor', 'm', 'faceAlpha', 0.5);
 %     % compute and display intersection
-%     [curves1, curves2] = intersectPlaneMesh(plane, v, f);
-%     drawPolyline3d(curves2, 'linewidth', 2, 'color', 'b')
+%     [curves, closed] = intersectPlaneMesh(plane, v, f);
+%     drawPolyline3d(curves(~closed), 'linewidth', 2, 'color', 'b')
 %
 %
 %   See also 
@@ -240,3 +245,8 @@
     poly = intersectionPoints(polyEdgeInds, :);
     rings = [rings, {poly}]; %#ok<AGROW>
 end
+
+
+%% Format output array
+polys = [rings, openPolys];
+closedFlag = [true(1, length(rings)), false(1, length(openPolys))];
diff --git a/tests/meshes3d/test_intersectPlaneMesh.m b/tests/meshes3d/test_intersectPlaneMesh.m
index b0c8a11a..1b992b81 100644
--- a/tests/meshes3d/test_intersectPlaneMesh.m
+++ b/tests/meshes3d/test_intersectPlaneMesh.m
@@ -1,5 +1,5 @@
 function test_suite = test_intersectPlaneMesh
-%TEST_INTERSECTPLANEMESH  Test case for the file intersectPlaneMesh
+%TEST_INTERSECTPLANEMESH  Test suite for the file intersectPlaneMesh.
 %
 %   Test case for the file intersectPlaneMesh
 %
@@ -11,12 +11,13 @@
 %
 % ------
 % Author: David Legland
-% e-mail: david.legland@grignon.inra.fr
+% e-mail: david.legland@inrae.fr
 % Created: 2012-07-31,    using Matlab 7.9.0.529 (R2009b)
 % Copyright 2012 INRA - Cepia Software Platform.
 
 test_suite = functiontests(localfunctions);
 
+
 function test_cube_horizPlane(testCase)
 plane = createPlane([5 5 5], [0 0 1]);
 [v, f] = createCube;
@@ -24,9 +25,9 @@ function test_cube_horizPlane(testCase)
 v = v * 10;
 polys = intersectPlaneMesh(plane, v, f);
 
-testCase.assertTrue(iscell(polys));
-testCase.assertEqual(1, length(polys));
-testCase.assertEqual(8, size(polys{1}, 1));
+assertTrue(testCase, iscell(polys));
+assertEqual(testCase, 1, length(polys));
+assertEqual(testCase, 8, size(polys{1}, 1));
 
 
 function test_cube_diagPlane(testCase)
@@ -36,9 +37,9 @@ function test_cube_diagPlane(testCase)
 v = v * 5;
 polys = intersectPlaneMesh(plane, v, f);
 
-testCase.assertTrue(iscell(polys));
-testCase.assertEqual(1, length(polys));
-testCase.assertEqual(6, size(polys{1}, 1));
+assertTrue(testCase, iscell(polys));
+assertEqual(testCase, 1, length(polys));
+assertEqual(testCase, 6, size(polys{1}, 1));
 
 
 function test_cube_planeOutside(testCase)
@@ -48,8 +49,8 @@ function test_cube_planeOutside(testCase)
 v = v * 5;
 polys = intersectPlaneMesh(plane, v, f);
 
-testCase.assertTrue(iscell(polys));
-testCase.assertEqual(0, length(polys));
+assertTrue(testCase, iscell(polys));
+assertEqual(testCase, 0, length(polys));
 
 
 function test_torus(testCase)
@@ -61,7 +62,7 @@ function test_torus(testCase)
 plane = [0 0 0  1 0 0  0 0 1];
 polys = intersectPlaneMesh(plane, v, f);
 
-testCase.assertEqual(2, length(polys));
+assertEqual(testCase, 2, length(polys));
 
 
 function test_verticalTriangle(testCase)
@@ -71,11 +72,12 @@ function test_verticalTriangle(testCase)
 f = [1 2 3];
 plane = [0 0 5  1 0 0  0 1 0];
 
-[rings, curves] = intersectPlaneMesh(plane, v, f);
+[polys, flags] = intersectPlaneMesh(plane, v, f);
 
-assertEqual(testCase, 0, length(rings));
-assertEqual(testCase, 1, length(curves));
-assertEqual(testCase, [2 3], size(curves{1}));
+assertEqual(testCase, 1, length(polys));
+assertEqual(testCase, 1, length(flags));
+assertEqual(testCase, false, flags);
+assertEqual(testCase, [2 3], size(polys{1}));
 
 
 function test_openCube(testCase)
@@ -86,9 +88,10 @@ function test_openCube(testCase)
 f = triangulateFaces(f(1:4, :));
 plane = [0 0 0.5  1 0 0  0 1 0];
 
-[rings, curves] = intersectPlaneMesh(plane, v, f);
+[polys, flags] = intersectPlaneMesh(plane, v, f);
 
-assertEqual(testCase, 0, length(rings));
-assertEqual(testCase, 2, length(curves));
-assertEqual(testCase, [3 3], size(curves{1}));
-assertEqual(testCase, [3 3], size(curves{2}));
+assertEqual(testCase, 2, length(polys));
+assertEqual(testCase, 2, length(flags));
+assertEqual(testCase, [false false], flags);
+assertEqual(testCase, [3 3], size(polys{1}));
+assertEqual(testCase, [3 3], size(polys{2}));