diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orientation.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orientation.h
index 4c39e7684b0..75dae3dd8da 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orientation.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orientation.h
@@ -26,7 +26,6 @@
 
 
 #include <algorithm>
-#include <CGAL/Polygon_mesh_processing/compute_normal.h>
 #include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
 #include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
 #include <CGAL/boost/graph/helpers.h>
@@ -40,23 +39,28 @@ namespace CGAL {
 namespace Polygon_mesh_processing {
 
 namespace internal{
-  template <class Less_xyz, class VPmap>
-  struct Compare_vertex_points_xyz_3{
-    Less_xyz less;
+
+  template <class GT, class VPmap>
+  struct Compare_vertex_points_z_3
+  {
+    typename GT::Compare_z_3 compare_z;
     VPmap vpmap;
 
-   Compare_vertex_points_xyz_3(VPmap const& vpmap)
-	: vpmap(vpmap){}
+    Compare_vertex_points_z_3(VPmap const& vpmap, const GT& gt)
+      : vpmap(vpmap)
+    {
+      compare_z = gt.compare_z_3_object();
+    }
 
     typedef bool result_type;
     template <class vertex_descriptor>
     bool operator()(vertex_descriptor v1, vertex_descriptor v2) const
     {
-      return less(get(vpmap, v1), get(vpmap, v2));
+      return CGAL::SMALLER == compare_z(get(vpmap, v1), get(vpmap, v2));
     }
-
   };
 
+
   template<typename Kernel, typename PM, typename NamedParameters>
   bool is_outward_oriented(typename boost::graph_traits<PM>::vertex_descriptor vd,
                            const PM& pmesh,
@@ -117,17 +121,55 @@ bool is_outward_oriented(const PolygonMesh& pmesh,
   VPMap vpmap = choose_param(get_param(np, vertex_point),
                              get_const_property_map(vertex_point, pmesh));
   //Kernel
-  typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type Kernel;
-
-  internal::Compare_vertex_points_xyz_3<typename Kernel::Less_xyz_3, VPMap >
-    less_xyz(vpmap);
+  typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GT;
+  GT gt = choose_param(get_param(np, geom_traits), GT());
 
+  //find the vertex with maximal z coordinate
   typename boost::graph_traits<PolygonMesh>::vertex_iterator vbegin, vend;
   cpp11::tie(vbegin, vend) = vertices(pmesh);
-  typename boost::graph_traits<PolygonMesh>::vertex_iterator v_min
-    = std::min_element(vbegin, vend, less_xyz);
 
-  return internal::is_outward_oriented<Kernel>(*v_min, pmesh, np);
+  internal::Compare_vertex_points_z_3<GT, VPMap> less_z(vpmap, gt);
+  typename boost::graph_traits<PolygonMesh>::vertex_iterator v_max_it
+    = std::max_element(vbegin, vend, less_z);
+  typename boost::graph_traits<PolygonMesh>::vertex_descriptor v_max = *v_max_it;
+
+  //among the incident edges to v_max,
+  // find one of the edges e with the minimal slope
+  typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
+  halfedge_descriptor min_slope_he = halfedge(v_max, pmesh);
+  CGAL_assertion(v_max == target(min_slope_he, pmesh));
+
+  typename GT::Compare_slope_3 compare_slope = gt.compare_slope_3_object();
+  BOOST_FOREACH(halfedge_descriptor he, halfedges_around_target(v_max, pmesh))
+  {
+    CGAL_assertion(v_max == target(min_slope_he, pmesh));
+    CGAL_assertion(v_max == target(he, pmesh));
+    if(CGAL::SMALLER == compare_slope(get(vpmap, v_max),
+                          get(vpmap, source(min_slope_he, pmesh)),
+                          get(vpmap, v_max),
+                          get(vpmap, source(he, pmesh))))
+    {
+      min_slope_he = he;
+    }
+  }
+
+  //around that edge, there are two faces(because the mesh is without borders),
+  //take the face f for which the other edge incident to f and v has the minimal slope
+  //Note we only keep the corresponding halfedge
+  if(CGAL::SMALLER != compare_slope(get(vpmap, v_max),
+                          get(vpmap, source(min_slope_he, pmesh)),
+                          get(vpmap, v_max),
+                          get(vpmap, source(prev(opposite(min_slope_he, pmesh), pmesh), pmesh))))
+    min_slope_he = opposite(min_slope_he, pmesh);
+
+  //check that the normal to f has z > 0
+  typename GT::Angle_3 angle_3 = gt.angle_3_object();
+  typename GT::Vector_3 vertical(0, 0, 1);
+
+  return CGAL::ACUTE == angle_3(get(vpmap, source(min_slope_he, pmesh)),
+                                get(vpmap, target(min_slope_he, pmesh)),
+                                get(vpmap, target(next(min_slope_he, pmesh), pmesh)),
+                                vertical);
 }
 
 ///\cond SKIP_IN_MANUAL