Update upstream source from tag 'upstream/4.1.0'

Update to upstream version '4.1.0'
with Debian dir f973f3f0739564c39cd0469e3c0bdbb441165c2d

13 files changed, 409 insertions, 21 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 2073c07..ed59347 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Name: geometry
-Version: 4.0.0
-Date: 03-02-2020
+Version: 4.1.0
+Date: 16-03-2024
 Author: Juan Pablo Carbajal <ajuanpi+dev@gmail.com>,
  Philip Nienhuis <prnienhuis@users.sf.net>,
  Simeon Simeonov <sss41@cam.ac.uk>,
diff --git a/INDEX b/INDEX
index 0f6d790..c375214 100644
--- a/INDEX
+++ b/INDEX
@@ -32,6 +32,10 @@ geometry >> Computational Geometry
  transformShape
 2D Others
  beltProblem
+3D Polygons
+ inpolyeder
+ inpolyeder2
+ isIntersectionInPolygon3D 
 Input [broken]
  @svg/svg
  @svg/plot
diff --git a/NEWS b/NEWS
index 54da527..9133f45 100644
--- a/NEWS
+++ b/NEWS
@@ -1,5 +1,13 @@
 Summary of important user-visible changes for releases of the geometry package
 ===============================================================================
+geometry-4.1.0   Release Date: 16-03-2024  Release Manager: Juan Pablo Carbajal
+===============================================================================
+* Added new functions to detect intersetions and points inside 3d polygons.
+* ispolycw that allows for two cell array inputs (px and py), see https://savannah.gnu.org/patch/?10297
+* Fixed compilation error with g++ v11, see https://savannah.gnu.org/bugs/?60884
+* Fixed failing test in joinPolygons
+
+===============================================================================
 geometry-4.0.0   Release Date: 03-02-2020  Release Manager: Juan Pablo Carbajal
 ===============================================================================
 As of this version, geometry does not contain any matgeom functionality.
diff --git a/inst/@svg/svg.m b/inst/@svg/svg.m
index dcf366a..0848136 100644
--- a/inst/@svg/svg.m
+++ b/inst/@svg/svg.m
@@ -1,4 +1,4 @@
-## Copyright (C) 2012-2019 Juan Pablo Carbajal
+## Copyright (C) 2012-2020 Juan Pablo Carbajal
 ##
 ## This program is free software; you can redistribute it and/or modify it under
 ## the terms of the GNU General Public License as published by the Free Software
@@ -61,22 +61,23 @@ function svg = svg(name='')
 
 endfunction
 
-%!test
+%!xtest
 %!  dc = svg('drawing5.svg');
 %!  dc.getpath();
 %!  dc.pathid();
 %!  dc.getpath('path3756');
-%!
+
+%!xtest
 %!  dc = svg('drawing.svg');
 %!  ids = dc.pathid();
 %!  dc.getpath({ids{[1 3]}});
 
-%!test
+%!xtest
 %!  dc = svg('drawing6.svg');
 %!  ids = dc.pathid();
 %!  P = dc.path2polygon(ids{1});
 
-%!test
+%!xtest
 %! dc = svg('drawing6.svg');
 %! dc.plot();
 %! dc.plot('color','r','linewidth',2);
diff --git a/inst/inpolyeder.m b/inst/inpolyeder.m
new file mode 100644
index 0000000..3a49fa5
--- /dev/null
+++ b/inst/inpolyeder.m
@@ -0,0 +1,136 @@
+## Copyright (C) 2014 Andreas Emch and Eduardo Hahn Paredes
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+% -*- texinfo -*-
+% @deftypefn {Function File} {[@var{in}, @var{on}] =} inpolyeder (@var{x}, @var{y}, @var{z}, @var{xv}, @var{yv}, @var{zv}, @var{doPlot}, @var{doPlotHelpers})
+% For a polyeder defined by vertex points @code{(@var{xv}, @var{yv}, @var{zv})},
+% determine if the points @code{(@var{x}, @var{y}), @var{z})} are inside or
+% outside the polyeder.
+% The variables @var{x}, @var{y}, @var{z} must have the same dimension. The optional
+% output @var{on} gives the points that are on the polyeder.
+% @seealso{inpolygon}
+% @end deftypefn
+
+% Author: Andreas Emch, Eduardo Hahn Paredes
+% Created: December 2012
+
+function [in, on] = inpolyeder (x, y, z, xv, yv, zv, doPlot, doPlotHelpers)
+
+    if (nargin < 6)
+        print_usage ();
+    endif
+
+    if (!(isreal(x) && isreal(y) && isreal(z)
+            && ismatrix(y) && ismatrix(y) && ismatrix(z)
+            && size_equal(x, y, z)))
+        error("inpolyeder: The first 3 vectors need to have the same size (test points)");
+    elseif (! (isreal(xv) && isreal(yv) && isreal(zv)
+            && isvector(xv) && isvector(yv) && isvector(zv)
+            && size_equal(xv, yv, zv)))
+        error("inpolyeder: The last 3 vectors need to have the same size (polyeder corners)");
+    endif
+
+    if (!isbool(doPlot))
+        error("inpolyeder: doPlot has to be a boolean value");
+    endif
+
+    if (!isbool(doPlotHelpers))
+        error("inpolyeder: doPlotHelpers has to be a boolean value");
+    endif
+
+    starttTime = cputime;
+
+    X = [xv, yv, zv];
+    K = convhulln(X);
+
+    P = [x y z];
+
+    on = zeros(size(x), "logical")';
+    in = zeros(size(x), "logical")';
+
+    if (doPlot == true)
+        clf
+        hold on
+        t = trisurf(K, X(:,1), X(:,2), X(:,3));
+        set(t,'facealpha',0.5)
+        m=unique(K);
+        plot3(X(m,1), X(m,2), X(m,3), 'ko', 'markerfacecolor', 'b');
+    endif
+
+    counter = 0;
+
+    for p1 = P'
+        counter++;
+
+        left = 0;
+        right = 0;
+
+        do
+            p2 = p1 + rand(3,1);
+        until (p1 != p2)
+
+        if (doPlotHelpers)
+            point1 = p1 + (-20 .* (p2-p1));
+            point2 = p1 + (20 .* (p2-p1));
+
+            plot3([point1(1,1); point2(1,1)], [point1(2,1); point2(2,1)], [point1(3,1); point2(3,1)], 'k-');
+        endif
+
+        for poly = K'
+            a = X(poly(1,1),:)';
+            b = X(poly(2,1),:)';
+            c = X(poly(3,1),:)';
+            [pointIn pointDistance point] = isIntersectionInPolygon3D(p1, p2, a, b, c);
+            if (pointIn)
+                if (pointDistance > 0)
+                    left += 1;
+                elseif (pointDistance < 0)
+                    right += 1;
+                elseif (pointDistance == 0)
+                    on(1,counter) = true;
+                endif
+            endif
+            if (doPlotHelpers && !isnan(point))
+                if (pointIn)
+                    plot3(point(1,1), point(2,1), point(3,1), 'ko', 'markerfacecolor', 'y');
+                else
+                    plot3(point(1,1), point(2,1), point(3,1), 'k*', 'markerfacecolor', 'y');
+                endif
+            endif
+        endfor
+        if (mod(left,2) == 1 && mod(right,2) == 1)
+            in(1,counter) = true;
+        endif
+    endfor
+
+    if (doPlot == true)
+        plot3(x(in'), y(in'), z(in'),'ko', 'markerfacecolor', 'g',
+            x(~in'), y(~in'), z(~in'),'ko', 'markerfacecolor', 'r',
+            x(on'), y(on'), z(on'),'ko', 'markerfacecolor', 'k')
+
+        x_max = max(max(X(:,1)), max(P(:,1))) + 0.5;
+        x_min = min(min(X(:,1)), min(P(:,1))) - 0.5;
+        y_max = max(max(X(:,2)), max(P(:,2))) + 0.5;
+        y_min = min(min(X(:,2)), min(P(:,2))) - 0.5;
+        z_max = max(max(X(:,3)), max(P(:,3))) + 0.5;
+        z_min = min(min(X(:,3)), min(P(:,3))) - 0.5;
+        axis ([x_min x_max y_min y_max z_min z_max]);
+        hold off
+    endif
+    printf('Total cpu time for checking %d points in a polyeder of %d polygons: %f seconds\n', size(x)(1,1), size(K)(1,1), cputime-starttTime);
+endfunction
diff --git a/inst/inpolyeder2.m b/inst/inpolyeder2.m
new file mode 100644
index 0000000..0744507
--- /dev/null
+++ b/inst/inpolyeder2.m
@@ -0,0 +1,145 @@
+## Copyright (C) 2014 Andreas Emch and Eduardo Hahn Paredes

+##

+## This file is part of Octave.

+##

+## Octave is free software; you can redistribute it and/or modify it

+## under the terms of the GNU General Public License as published by

+## the Free Software Foundation; either version 3 of the License, or (at

+## your option) any later version.

+##

+## Octave is distributed in the hope that it will be useful, but

+## WITHOUT ANY WARRANTY; without even the implied warranty of

+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

+## General Public License for more details.

+##

+## You should have received a copy of the GNU General Public License

+## along with Octave; see the file COPYING.  If not, see

+## <http://www.gnu.org/licenses/>.

+

+% -*- texinfo -*-

+% @deftypefn {Function File} {[@var{in}, @var{on}] =} inpolyeder (@var{x}, @var{y}, @var{z}, @var{xv}, @var{yv}, @var{zv}, @var{xk}, @var{yk}, @var{zk}, @var{doPlot}, @var{doPlotHelpers})

+% For a polyeder defined by vertex points @code{(@var{xv}, @var{yv}, @var{zv})},

+% determine if the points @code{(@var{x}, @var{y}), @var{z})} are inside or 

+% outside the polyeder.

+% The variables @var{x}, @var{y}, @var{z} must have the same dimension. The optional

+% output @var{on} gives the points that are on the polyeder.

+% @seealso{inpolygon}

+% @end deftypefn

+

+% Author: Andreas Emch, Eduardo Hahn Paredes

+% Created: December 2012

+

+function [in, on] = inpolyeder2 (x, y, z, xv, yv, zv, xk, yk, zk, doPlot, doPlotHelpers)

+

+    % Check number of arguments

+    if (nargin < 9)

+        print_usage ();

+    endif

+

+    if (!(isreal(x) && isreal(y) && isreal(z)

+            && ismatrix(y) && ismatrix(y) && ismatrix(z)

+            && size_equal(x, y, z)))

+        error("inpolyeder2: The first 3 vectors need to have the same size (test points)");

+    elseif (! (isreal(xv) && isreal(yv) && isreal(zv)

+            && isvector(xv) && isvector(yv) && isvector(zv)

+            && size_equal(xv, yv, zv)))

+        error("inpolyeder2: Those 3 vectors need to have the same size (polyeder corners)");

+    elseif (! (isreal(xk) && isreal(yk) && isreal(zk)

+            && isvector(xk) && isvector(yk) && isvector(zk)

+            && size_equal(xk, yk, zk)))

+        error("inpolyeder2: The last 3 vectors need to have the same size (surface definers)");

+    endif

+

+    if (!isbool(doPlot))

+        error("inpolyeder: doPlot has to be a boolean value");

+    endif

+

+    if (!isbool(doPlotHelpers))

+        error("inpolyeder: doPlotHelpers has to be a boolean value");

+    endif

+

+    starttTime = cputime;

+

+    X = [xv, yv, zv];

+    K = [xk, yk, zk];

+    P = [x y z];

+

+    on = zeros(size(x), "logical")';

+    in = zeros(size(x), "logical")';

+

+    if (doPlot == true)

+        clf 

+        hold on

+        t = trisurf(K, X(:,1), X(:,2), X(:,3));

+        set(t,'facealpha',0.5, 'markerfacecolor', 'g')

+        m=unique(K);

+        plot3(X(m,1), X(m,2), X(m,3), 'ko', 'markerfacecolor', 'b');

+    endif

+

+    counter = 0;

+

+    for p1 = P'

+        counter++;

+

+        left = 0;

+        right = 0;

+

+        do

+            p2 = p1 + rand(3,1);

+        until (p1 != p2)

+    

+        if (doPlotHelpers)

+            point1 = p1 + (-20 .* (p2-p1));

+            point2 = p1 + (20 .* (p2-p1));

+        

+            plot3([point1(1,1); point2(1,1)], [point1(2,1); point2(2,1)], [point1(3,1); point2(3,1)], 'k-');

+        endif

+

+        for poly = K'

+            a = X(poly(1,1),:)';

+            b = X(poly(2,1),:)';

+            c = X(poly(3,1),:)';

+        

+            [pointIn pointDistance point] = isIntersectionInPolygon3D(p1, p2, a, b, c);

+        

+            if (pointIn)

+                if (pointDistance > 0)

+                    left += 1;

+                elseif (pointDistance < 0)

+                    right += 1;

+                elseif (pointDistance == 0)

+                    on(1,counter) = true;

+                endif

+            endif

+        

+            if (doPlotHelpers && !isnan(point))

+                if (pointIn)

+                    plot3(point(1,1), point(2,1), point(3,1), 'ko', 'markerfacecolor', 'y');

+                else

+                    plot3(point(1,1), point(2,1), point(3,1), 'k*', 'markerfacecolor', 'y');

+                endif

+            endif

+        endfor

+        if (mod(left,2) == 1 && mod(right,2) == 1)

+            in(1,counter) = true;

+        endif

+    endfor

+

+    if (doPlot == true)

+        plot3(

+            x(in'), y(in'), z(in'),'ko', 'markerfacecolor', 'g',

+            x(~in'), y(~in'), z(~in'),'ko', 'markerfacecolor', 'r',

+            x(on'), y(on'), z(on'),'ko', 'markerfacecolor', 'k'

+        )

+

+        x_max = max(max(X(:,1)), max(P(:,1))) + 0.5;

+        x_min = min(min(X(:,1)), min(P(:,1))) - 0.5;

+        y_max = max(max(X(:,2)), max(P(:,2))) + 0.5;

+        y_min = min(min(X(:,2)), min(P(:,2))) - 0.5;

+        z_max = max(max(X(:,3)), max(P(:,3))) + 0.5;

+        z_min = min(min(X(:,3)), min(P(:,3))) - 0.5;

+        axis ([x_min x_max y_min y_max z_min z_max]);

+        hold off

+    endif

+    printf('Total cpu time for checking %d points in a polyeder of %d polygons: %f seconds\n', size(x)(1,1), size(K)(1,1), cputime-starttTime);

+endfunction

diff --git a/inst/isIntersectionInPolygon3D.m b/inst/isIntersectionInPolygon3D.m
new file mode 100644
index 0000000..e2b61ca
--- /dev/null
+++ b/inst/isIntersectionInPolygon3D.m
@@ -0,0 +1,86 @@
+## Copyright (C) 2014 Andreas Emch and Eduardo Hahn Paredes

+##

+## This file is part of Octave.

+##

+## Octave is free software; you can redistribute it and/or modify it

+## under the terms of the GNU General Public License as published by

+## the Free Software Foundation; either version 3 of the License, or (at

+## your option) any later version.

+##

+## Octave is distributed in the hope that it will be useful, but

+## WITHOUT ANY WARRANTY; without even the implied warranty of

+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

+## General Public License for more details.

+##

+## You should have received a copy of the GNU General Public License

+## along with Octave; see the file COPYING.  If not, see

+## <http://www.gnu.org/licenses/>.

+

+% -*- texinfo -*-

+% @deftypefn {Function File} {[@var{in}, @var{distance}] =} isIntersectionInPolygon3D (@var{p1}, @var{p2}, @var{a}, @var{b}, @var{c})

+% For a line, defined by two vector-points @code{(@var{p1}, @var{p2})},

+% calculate the intersection-point with the plane, defined by three vector-points @code{(@var{a}, @var{b}), @var{c})},

+% and determine if this intersection-point is inside or on the polygon defined by three vector-points @code{(@var{a}, @var{b}), @var{c})}

+% The variables @var{p1}, @var{p2}, @var{a}, @var{b}, @var{c} are all points in the 3D-dimension.

+% @seealso{inpolygon}

+% @end deftypefn

+

+% Author: Andreas Emch, Eduardo Hahn Paredes

+% Created: January 2013

+

+function [in distance point] = isIntersectionInPolygon3D (p1, p2, a, b, c)

+    u = b-a;

+    v = c-a;

+

+    E = [a u v];

+    n = cross(u',v');

+    

+    if (u == v)

+        in = false;

+        distance = NaN;

+        point = NaN;

+        return;

+    endif

+    

+    x = n(1,1);

+    y = n(1,2);

+    z = n(1,3);

+    d = -(a(1,1) * x + a(2,1) * y + a(3,1) * z);

+    

+    plane = 0;

+    

+    if (E(1,2) == 0 && E(1,3) == 0)

+        plane = 1;

+    elseif (E(2,2) == 0 && E(2,3) == 0)

+        plane = 2;

+    elseif (E(3,2) == 0 && E(3,3) == 0)

+        plane = 3;

+    else

+        plane = 0;

+    endif

+    

+    P = [p1 p2-p1];

+    t = -(x * P(1,1) + y * P(2,1) + z * P(3,1) + d) / (x * P(1,2) + y * P(2,2) + z * P(3,2));

+

+    dP = p1 + (t .* (p2-p1));

+    xPoly = [a(1,1) b(1,1) c(1,1) a(1,1)];

+    yPoly = [a(2,1) b(2,1) c(2,1) a(2,1)];

+    zPoly = [a(3,1) b(3,1) c(3,1) a(3,1)];    

+

+    [in1 on1] = inpolygon(dP(1,1), dP(2,1), xPoly, yPoly);

+    [in2 on2] = inpolygon(dP(2,1), dP(3,1), yPoly, zPoly);

+    [in3 on3] = inpolygon(dP(3,1), dP(1,1), zPoly, xPoly);

+    

+    if ((plane == 0 && (in1 == 1 || on1 == 1) && (in2 == 1 || on2 == 1) && (in3 == 1 || on3 == 1))

+        ||(plane == 1 && (in2 == 1 || on2 == 1))

+        ||(plane == 2 && (in3 == 1 || on3 == 1))

+        ||(plane == 3 && (in1 == 1 || on1 == 1)))

+        in = true;

+        distance = t;

+        point = dP;

+    else

+        in = false;

+        distance = NaN;

+        point = dP;

+    endif

+endfunction

diff --git a/inst/ispolycw.m b/inst/ispolycw.m
index 45dbffc..253fb4e 100644
--- a/inst/ispolycw.m
+++ b/inst/ispolycw.m
@@ -1,5 +1,6 @@
 ## Copyright (C) 2016 - Juan Pablo Carbajal
 ## Copyright (C) 2017 - Piyush Jain
+## Copyright (C) 2022 - Nir Krakauer 
 ##
 ## This program is free software; you can redistribute it and/or modify
 ## it under the terms of the GNU General Public License as published by
@@ -17,19 +18,19 @@
 ## -*- texinfo -*-
 ## @deftypefn {} {@var{ccw} =} ispolycw (@var{p})
 ## @deftypefnx {} {@var{ccw} =} ispolycw (@var{px}, @var{py})
-## Returns true if the polygon @var{p} are oriented Clockwise.
+## Returns true if the polygon @var{p} are oriented clockwise.
 ##
 ## @var{p} is a N-by-2 array containing coordinates of vertices. The coordinates
-## of the vertices of the polygon can also be given as two N-by-1 arrways
+## of the vertices of the polygon can also be given as two N-by-1 arrays
 ## @var{px}, @var{py}.
 ##
 ## If polygon is self-crossing, the result is undefined.
 ##
-## If x and y contain multiple contours, either in NaN-separated vector form or in cell array form, ispolycw returns a logical array containing one true or false value per contour.
+## If @var{px}, @var{py} contain multiple contours, either in NaN-separated vector form or in cell array form, ispolycw returns a logical array containing one true or false value per contour.
 ##
 ## If a contour contains two or fewer vertices, ispolycw returns true.
 ##
-## If @var{points} is a cell, each element is considered a polygon, the
+## If @var{p} (or @var{px}, @var{py}) is a cell, each element is considered a polygon, the
 ## resulting @var{cw} array has the same shape as the cell.
 ##
 ## @seealso{polygonArea}
@@ -38,7 +39,11 @@
 function cw = ispolycw (px, py)
 
   if iscell (px)
-    cw = cellfun (@ispolycw, px);
+    if nargin == 2
+      cw = cellfun (@ispolycw, px, py);
+    else
+      cw = cellfun (@ispolycw, px);
+    endif
   else
     
     if nargin == 2;
@@ -70,6 +75,8 @@ end
 %!assert (ispolycw (pcw));
 %!assert (ispolycw ({pccw;pcw}), [false;true]);
 %!assert (ispolycw ({pccw,pcw}), [false,true]);
+%!assert (ispolycw ({pccw(:, 1);pcw(:, 1)}, {pccw(:, 2);pcw(:, 2)}), [false;true]);
+%!assert (ispolycw ({pccw(:, 1),pcw(:, 1)}, {pccw(:, 2),pcw(:, 2)}), [false,true]);
 %!assert (ispolycw(ph),[false;true]);
 
 %!test
diff --git a/inst/joinPolygons.m b/inst/joinPolygons.m
index 59c2127..3b9d711 100644
--- a/inst/joinPolygons.m
+++ b/inst/joinPolygons.m
@@ -57,5 +57,5 @@ endfunction
 %! XY = joinPolygons ({[1 6; 2 5; 3 4]; [4 3; 5 2; 6 1]});
 %! assert (XY, [1 6; 2 5; 3 4; NaN NaN; 4 3; 5 2; 6 1]);
 
-%!error <joinPolygons: cell array expected> joinPolygons ([1 2 NaN 3 4], [56 NaN 78])
+%!error <joinPolygons: function called with too many inputs> joinPolygons ([1 2 NaN 3 4], [56 NaN 78])
 %!error <joinPolygons: column vectors expected>  joinPolygons ({[1,0], [0,2]});
diff --git a/src/clipper.cpp b/src/clipper.cpp
index 1b10300..7c9ef05 100644
--- a/src/clipper.cpp
+++ b/src/clipper.cpp
@@ -747,7 +747,8 @@ void DisposeOutPts(OutPt*& pp)
 
 inline void InitEdge(TEdge* e, TEdge* eNext, TEdge* ePrev, const IntPoint& Pt)
 {
-  std::memset(e, 0, sizeof(TEdge));
+  //std::memset(e, 0, sizeof(TEdge));
+  *e = {};
   e->Next = eNext;
   e->Prev = ePrev;
   e->Curr = Pt;
diff --git a/src/martinez.cpp b/src/martinez.cpp
index 52558b3..d981597 100644
--- a/src/martinez.cpp
+++ b/src/martinez.cpp
@@ -26,7 +26,7 @@ void Martinez::print (SweepEvent& e)
 
 // Compare two sweep events
 // Return true means that e1 is placed at the event queue after e2, i.e,, e1 is processed by the algorithm after e2
-bool Martinez::SweepEventComp::operator() (SweepEvent* e1, SweepEvent* e2) {
+bool Martinez::SweepEventComp::operator() (SweepEvent* e1, SweepEvent* e2) const {
 	if (e1->p.x > e2->p.x) // Different x-coordinate
 		return true;
 	if (e2->p.x > e1->p.x) // Different x-coordinate
@@ -40,7 +40,7 @@ bool Martinez::SweepEventComp::operator() (SweepEvent* e1, SweepEvent* e2) {
 }
 
 // e1 and a2 are the left events of line segments (e1->p, e1->other->p) and (e2->p, e2->other->p)
-bool Martinez::SegmentComp::operator() (SweepEvent* e1, SweepEvent* e2) {
+bool Martinez::SegmentComp::operator() (SweepEvent* e1, SweepEvent* e2) const {
 	if (e1 == e2)
 		return false;
 	if (signedArea (e1->p, e1->other->p, e2->p) != 0 || signedArea (e1->p, e1->other->p, e2->other->p) != 0) {
diff --git a/src/martinez.h b/src/martinez.h
index 3275171..307f5bd 100644
--- a/src/martinez.h
+++ b/src/martinez.h
@@ -38,7 +38,7 @@ private:
 
 	struct SweepEvent;
 	struct SegmentComp : public binary_function<SweepEvent*, SweepEvent*, bool> { // for sorting edges in the sweep line
-		bool operator() (SweepEvent* e1, SweepEvent* e2);
+		bool operator() (SweepEvent* e1, SweepEvent* e2) const;
 	};
 	
 	struct SweepEvent {
@@ -65,7 +65,7 @@ private:
 	static void print (SweepEvent& e); // This function is intended for debugging purposes
 
 	struct SweepEventComp : public binary_function<SweepEvent*, SweepEvent*, bool> { // for sortening events
-		bool operator() (SweepEvent* e1, SweepEvent* e2);
+		bool operator() (SweepEvent* e1, SweepEvent* e2) const;
 	};
 
 	/** @brief Event Queue */
diff --git a/src/polygon.cpp b/src/polygon.cpp
index 63628c5..17b4cff 100644
--- a/src/polygon.cpp
+++ b/src/polygon.cpp
@@ -94,7 +94,7 @@ void Polygon::move (double x, double y)
 namespace { // start of anonymous namespace
 	struct SweepEvent;
 	struct SegmentComp : public binary_function<SweepEvent*, SweepEvent*, bool> {
-		bool operator() (SweepEvent* e1, SweepEvent* e2);
+		bool operator() (SweepEvent* e1, SweepEvent* e2) const;
 	};
 
 	struct SweepEvent {
@@ -117,7 +117,7 @@ namespace { // start of anonymous namespace
 	};
 
 	struct SweepEventComp : public binary_function<SweepEvent*, SweepEvent*, bool> {
-		bool operator() (SweepEvent* e1, SweepEvent* e2) {
+		bool operator() (SweepEvent* e1, SweepEvent* e2) const {
 			if (e1->p.x < e2->p.x) // Different x coordinate
 				return true;
 			if (e2->p.x < e1->p.x) // Different x coordinate
@@ -133,7 +133,7 @@ namespace { // start of anonymous namespace
 
 } // end of anonymous namespace
 
-bool SegmentComp::operator() (SweepEvent* e1, SweepEvent* e2) {
+bool SegmentComp::operator() (SweepEvent* e1, SweepEvent* e2) const {
 	if (e1 == e2)
 		return false;
 	if (signedArea (e1->p, e1->other->p, e2->p) != 0 || signedArea (e1->p, e1->other->p, e2->other->p) != 0) {