integrated PolyFit

Alpha_shapes_2/include/CGAL/Alpha_shape_2.h

@@ -319,7 +319,7 @@ public:
       return n;
     }
 
-private :
+protected :
 
   //--------------------- INITIALIZATION OF PRIVATE MEMBERS -----------
 

Documentation/doc/Documentation/Installation.txt

@@ -622,6 +622,24 @@ and, by extension, of surface meshes (see Section \ref BGLPartitioning of the pa
 More information is available on the METIS library
 at <A HREF="http://glaros.dtc.umn.edu/gkhome/metis/metis/overview">`http://glaros.dtc.umn.edu/gkhome/metis/metis/overview`</A>.
 
+\subsection thirdpartyGLPK GLPK
+
+\sc{GLPK} (GNU Linear Programming Kit) is a library for solving linear programming (LP), mixed integer programming (MIP), and other related problems.
+
+In \cgal, \sc{GLPK} provides an optional linear integer program solver in the \ref PkgPolygonalSurfaceReconstruction package.
+
+The \sc{GLPK} web site is <A HREF="https://www.gnu.org/software/glpk/">`https://www.gnu.org/software/glpk/`</A>.
+
+
+\subsection thirdpartySCIP SCIP 
+
+\sc{SCIP} (Solving Constraint Integer Programs) is currently one of the fastest open source solvers for mixed integer programming (MIP) and mixed integer nonlinear programming (MINLP).
+
+In \cgal, \sc{SCIP} provides an optional linear integer program solver in the \ref PkgPolygonalSurfaceReconstruction package.
+
+The \sc{SCIP} web site is <A HREF="http://scip.zib.de/">`http://scip.zib.de/`</A>.
+
+
 \section secbuilding Building CGAL 
 
 The results of a successful configuration are build files that control the build step.

Documentation/doc/Documentation/dependencies

@@ -2,6 +2,7 @@ Algebraic_foundations
 Advancing_front_surface_reconstruction
 AABB_tree
 Polygon
+Polygonal_surface_reconstruction
 Number_types
 Alpha_shapes_2
 Alpha_shapes_3

Documentation/doc/Documentation/packages.txt

@@ -97,6 +97,7 @@
 \package_listing{Poisson_surface_reconstruction_3}
 \package_listing{Scale_space_reconstruction_3}
 \package_listing{Advancing_front_surface_reconstruction}
+\package_listing{Polygonal_surface_reconstruction}
 \package_listing{Optimal_transportation_reconstruction_2}
 
 \cgalPackageSection{PartGeometryProcessing,Geometry Processing}

Documentation/doc/biblio/geom.bib

@@ -109344,6 +109344,13 @@ problems in computational geometry is presented."
 , update =	"96.01 held+mitchell"
 }
 
+@article{nan2017polyfit,    
+     title = {PolyFit: Polygonal Surface Reconstruction from Point Clouds},    
+     author = {Nan, Liangliang and Wonka, Peter},    
+     journal = {ICCV},        
+     year = {2017} 
+}
+
 @article{nhl-merp-84
 , author =	"A. Naamad and W.-L. Hsu and D. T. Lee"
 , title =	"On the maximum empty rectangle problem"

Installation/AUTHORS

@@ -68,6 +68,7 @@ Kurt Mehlhorn
 Naceur Meskini
 Andreas Meyer
 Michal Meyerovitch
+Liangliang Nan
 Oren Nechushtan
 Gabriele Neyer
 Ralf Osbild

Installation/include/CGAL/license/Polygonal_surface_reconstruction.h

@@ -0,0 +1,67 @@
+// Copyright (c) 2016  GeometryFactory SARL (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public License as
+// published by the Free Software Foundation; either version 3 of the License,
+// or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0+
+//
+// Author(s) : Andreas Fabri
+//
+// Warning: this file is generated, see include/CGAL/licence/README.md
+
+
+#ifndef CGAL_LICENSE_POLYGONAL_SURFACE_RECONSTRUCTION_H
+#define CGAL_LICENSE_POLYGONAL_SURFACE_RECONSTRUCTION_H
+
+#include <CGAL/config.h>
+#include <CGAL/license.h>
+
+
+
+
+#ifdef CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_COMMERCIAL_LICENSE
+
+#  if CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_COMMERCIAL_LICENSE < CGAL_RELEASE_DATE
+
+#    if defined(CGAL_LICENSE_WARNING)
+
+       CGAL_pragma_warning("Your commercial license for CGAL does not cover "
+                           "this release of the Polygonal Surface Reconstruction package.")
+#    endif
+
+#    ifdef CGAL_LICENSE_ERROR
+#      error "Your commercial license for CGAL does not cover this release \
+of the Polygonal Surface Reconstruction package. \
+You get this error, as you defined CGAL_LICENSE_ERROR."
+#    endif // CGAL_LICENSE_ERROR
+
+#  endif // CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_COMMERCIAL_LICENSE < CGAL_RELEASE_DATE
+
+#else // no CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_COMMERCIAL_LICENSE
+
+#  if defined(CGAL_LICENSE_WARNING)
+     CGAL_pragma_warning("\nThe macro CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_COMMERCIAL_LICENSE is not defined."
+                          "\nYou use the CGAL Polygonal Surface Reconstruction package under "
+                          "the terms of the GPLv3+.")
+#  endif // CGAL_LICENSE_WARNING
+
+#  ifdef CGAL_LICENSE_ERROR
+#    error "The macro CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_COMMERCIAL_LICENSE is not defined.\
+ You use the CGAL Polygonal Surface Reconstruction package under the terms of \
+the GPLv3+. You get this error, as you defined CGAL_LICENSE_ERROR."
+#  endif // CGAL_LICENSE_ERROR
+
+#endif // no CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_COMMERCIAL_LICENSE
+
+#endif // CGAL_LICENSE_CHECK_POLYGONAL_SURFACE_RECONSTRUCTION_H

Polygonal_surface_reconstruction/doc/Polygonal_surface_reconstruction/Doxyfile.in

@@ -0,0 +1,3 @@
+@INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS}
+PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - Polygonal Surface Reconstruction"
+EXCLUDE = ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/internal

Polygonal_surface_reconstruction/doc/Polygonal_surface_reconstruction/PackageDescription.txt

@@ -0,0 +1,36 @@
+// Polygonal Surface Reconstruction should equal the project title in Doxyfile.in
+
+/// \defgroup PkgPolygonalSurfaceReconstruction Polygonal Surface Reconstruction Reference
+
+/*!
+\addtogroup PkgPolygonalSurfaceReconstruction
+
+\cgalPkgDescriptionBegin{Polygonal Surface Reconstruction, PkgPolygonalSurfaceReconstructionSummary}
+\cgalPkgPicture{polyfit.png}
+
+\cgalPkgSummaryBegin
+\cgalPkgAuthors{Liangliang Nan}
+\cgalPkgDesc{This package provides a method for piecewise planar object reconstruction from point clouds. 
+The method takes as input an unordered point set (and optionally planar segments) of a piecewise planar 
+object and outputs a lightweight and watertight surface mesh interpolating the input point set. The
+method can handle arbitrary piecewise planar objects and is capable of recovering sharp features and is
+robust to noise and outliers.}
+\cgalPkgManuals{Chapter_PolygonalSurfaceReconstruction, PkgPolygonalSurfaceReconstruction}
+\cgalPkgSummaryEnd
+
+\cgalPkgShortInfoBegin
+\cgalPkgSince{4.14}
+\cgalPkgDependsOn{\ref PkgAlphaShape2Summary and \ref PkgSolverSummary}
+\cgalPkgBib{cgal:x-x}
+\cgalPkgLicense{\ref licensesGPL "GPL"}
+
+\cgalPkgShortInfoEnd
+
+\cgalPkgDescriptionEnd
+
+\cgalClassifedRefPages
+
+## Classes ##
+- `CGAL::Polygonal_surface_reconstruction<Gt>`
+
+*/

Polygonal_surface_reconstruction/doc/Polygonal_surface_reconstruction/Polygonal_surface_reconstruction.txt

@@ -0,0 +1,228 @@
+namespace CGAL {
+/*!
+
+\mainpage User Manual
+\anchor Chapter_PolygonalSurfaceReconstruction
+
+\cgalAutoToc
+
+\author Liangliang Nan
+
+\section secIntro Introduction
+
+This package implements a hypothesis-and-selection based method for piecewise planar object
+reconstruction from point clouds \cgalCite{nan2017polyfit}. The method takes as input an
+unordered point set sampled from a piecewise planar object. The output is a compact and
+watertight surface mesh interpolating the input point set. The method assumes that
+all necessary major planes are provided (or can be extracted from the input point set using the
+shape detection method described in \ref PkgPointSetShapeDetection3Summary, or any other
+alternative methods).
+
+The existing surface reconstruction methods available in \cgal (i.e.,
+\ref PkgPoissonSurfaceReconstructionSummary, \ref PkgAdvancingFrontSurfaceReconstructionSummary, and
+\ref PkgScaleSpaceReconstruction3Summary) are
+suitable for point set representing objects described by smooth surfaces. For man-made objects
+such as buildings, the results are not satisfactory due to the imperfections
+and complexity of the reconstructed models (i.e., gigantic meshes, missing regions, noises, and
+undesired structures). This is mainly because these methods tend to closely follow the surface
+details. The algorithm introduced in this package is intended to produce simplified and watertight
+surface meshes. It can recover sharp features of the objects, and it can handle a large amount of
+noise and outliers, complementing other surface reconstruction methods.
+
+A mixed integer programming (MIP) solver is required to solve the constrained optimization problem
+formulated by the method (see \ref PkgSolverSummary).
+
+\section secAlgo Algorithm
+
+Unlike traditional piecewise planar object reconstruction methods that focus on either extracting good
+geometric primitives or obtaining proper arrangements of primitives, the emphasis of this method lies
+in intersecting the primitives (i.e. planes) and seeking
+for an appropriate combination of them to obtain a manifold and watertight polygonal surface model.
+
+The method casts surface reconstruction as a binary labeling problem based on a hypothesizing-and-selection
+strategy. The reconstruction consists in the following steps:
+-# extract planes from the input point set (can be skipped if planes are known or provided by other means);
+-# generate a set of candidate faces by intersecting the extracted planar primitives;
+-# select an optimal subset of the candidate faces through optimization under hard constraints that enforces the final polygonal surface to
+be topologically correct.
+
+\cgalFigureBegin{polyfit, polyfit_pipeline.png}
+Polygonal surface reconstruction
+
+(a) Input point set; (b) Extracted planar segments; (c) Candidate faces generated by pairwise intersection;
+(d) Faces selected through optimization; (e) Reconstructed model.
+\cgalFigureEnd
+
+
+
+
+
+\subsection subsecAlgoEnergyTerms Energy Terms
+
+Given \f$ N \f$  candidate faces \f$F = \{f_i | 1 \leq i \leq N\}\f$ generated by pairwise intersection, the
+optimal subset of the candidate faces that best describes the geometry of the object and form a manifold and
+watertight polygonal surface is selected through optimization.
+
+Let \f$X = \{x_i | 1 \leq i \leq N\}\f$ denote the selection of the faces (i.e., \f$f_i\f$ is chosen
+if \f$x_i = 1\f$ or not chosen if \f$x_i = 0\f$), the objective function consists of three energy terms:
+data-fitting, model complexity, and point coverage.
+
+- Data-fitting. This term is intended to evaluate the fitting quality of the faces to the point cloud. It is
+defined to measure a confidence-weighted percentage of points that do not contribute to the final reconstruction.
+
+    \f$
+  		\begin{equation}
+			E_f = 1 - \frac{1}{|P|} \sum_{i=1}^N x_i \cdot support(f_i),
+			\label{eq:datafitting}
+		\end{equation}
+    \f$
+
+  \f$|P|\f$ is the total number of points in \f$P\f$. \f$support(f_i)\f$ is the number of supporting points
+  accounting for an appropriate notion of confidence.
+
+
+- Model complexity. This term is to encourage simple structures (i.e., large planar regions).
+It is defined as the ratio of sharp edges in the model.
+
+    \f$
+  		\begin{equation}
+		E_m = \frac{1}{|E|}\sum_{i=1}^{|E|} corner(e_i),
+		\label{eq:complexity}
+		\end{equation}
+	\f$
+
+	\f$|E|\f$ denotes the total number of pairwise intersections in the candidate face set.
+	\f$corner(e_i)\f$ is an indicator function whose value is determined by the configuration of the two selected
+	faces connected by \f$e_i\f$. \f$corner(e_i)\f$ will have a value of 1 if the faces associated with \f$e_i\f$
+	introduce a sharp edge in the final model. Otherwise, \f$corner(e_i)\f$ has a zero value meaning the two faces
+	are coplanar.
+
+- Point coverage. This term is intended to handle missing data. The idea is to keep the unsupported regions
+  (i.e., regions not covered by points) of the final model as small as possible. This term is defined as the ratio of
+  uncovered regions in the model.
+
+    \f$
+  			\begin{equation}
+			E_c = \frac{1}{area(M)}\sum_{i=1}^N x_i \cdot (area(f_i) - area(M_i^\alpha)),
+			\label{eq:coverage}
+			\end{equation}
+	\f$
+
+  Here \f$area(M)\f$, \f$area(f_i)\f$, and \f$area(M_i^\alpha)\f$ denote the surface areas of the final model,
+  a candidate face \f$f_i\f$, and the \f$\alpha\f$-shape mesh \f$M_i^\alpha\f$ of \f$f_i\f$, respectively.
+  The \f$area(M_i^\alpha)\f$ is to approximate the area of the face covered by the 3D points.
+
+
+  For more details on the energy terms and the formulation, please refer to the original paper \cgalCite{nan2017polyfit}.
+
+
+\subsection subsecAlgoFaceSelection Face Selection
+
+With the above energy terms, the optimal set of faces can be obtained by minimizing a weighted sum of these terms
+under certain hard constraints.
+
+\f$
+    \begin{equation}
+    \begin{aligned}
+    \underset{\mathbf{X}}{\text{min}} \quad & \lambda_f \cdot E_f + \lambda_m \cdot E_m +  \lambda_c \cdot E_c   \\
+    \text{s.t.}                       \quad & \begin{cases}
+    \begin{aligned}
+    & \sum_{j \in \mathcal{N}(e_i)} {x_j} = 2 \quad \text{or} \quad 0,     && 1 \leq i \leq |E|   \\
+    & x_i \in \{0, 1\},                                                && 1 \leq i \leq N
+    \end{aligned}
+    \end{cases}
+    \end{aligned}
+    \label{eq:optimization}
+    \end{equation}
+\f$
+
+Here \f$ \sum_{j \in \mathcal{N}(e_i)} {x_j} \f$ counts the number of faces connected by an edge \f$ e_i \f$.
+This value is enforced to be either 0 or 2, meaning none or two of the faces can be selected. These hard constraints
+guarantee that each edge of the final model connects only two adjacent faces.
+
+\cgalFigureBegin{intersection, intersection.png}
+Intersection of two planes
+
+Two planes intersect with each other resulting in four parts (a). (b) to (g) show all possible combinations to
+ensure a 2-manifold surface. The edge in (b) and (c) connects two co-planar parts, while in (d) to (g) it introduces
+sharp edges in the final model.
+\cgalFigureEnd
+
+
+\section secPSRexamples Examples
+
+\subsection subsecExampleNoInputPlane Reconstruction from Point Clouds
+
+The method assumes that all necessary planes can be extracted from the input point set. The following example
+first extracts planes from the input point cloud and then reconstructs the surface model.
+
+\cgalExample{Polygonal_surface_reconstruction/polyfit_example_without_input_planes.cpp}
+
+
+
+
+\subsection subsecExampleUserProvidedPlanes Reconstruction with User-Provided Planes
+
+The following example shows the reconstruction using user-provided planar segments.
+
+\cgalExample{Polygonal_surface_reconstruction/polyfit_example_user_provided_planes.cpp}
+
+
+
+\subsection subsecExampleDetailControl Model Complexity Control
+In addition to favoring clean and compact reconstruction results by encouraging large planar regions,
+the model complexity term also provides control over the model details, i.e., increasing the influence of
+this term results in fewer details in the reconstructed 3D models.
+
+\cgalFigureBegin{complexity, complexity.png}
+The effect of the model complexity term
+
+Reconstruction by gradually increasing the influence of the model complexity term.
+The values under each subfigure are the weights used in the corresponding optimization.
+\cgalFigureEnd
+
+
+The following example shows how to control the model complexity by tuning the weight of the model complexity term.
+
+\remark
+This example also shows how to reuse the intermediate results from the candidate generation step.
+
+\cgalExample{Polygonal_surface_reconstruction/polyfit_example_model_complexty_control.cpp}
+
+
+\section secPerformances Performance
+
+\subsection subsecComplexity The Problem Complexity
+
+The method is intended for reconstructing single objects with reasonable geometric complexity.
+The current implementation computes pairwise intersections of the planar segments, which is sufficient
+but not necessary to ensure topologically accurate reconstructions. Running on large complex objects
+may result in an extremely large number of candidate faces, and thus a huge integer programming problem
+to be solved.
+
+The reconstruction time of a single object with moderate complexity is typically within a few seconds.
+Among the three steps, the face selection step deminates the reconstruction pipeline when the number
+of candidate faces is large (e.g., more than 5,000).
+
+
+
+
+
+\subsection subsecSolver The Numerical Solver
+The current implementation incorporates two open source solvers: \ref thirdpartyGLPK and \ref thirdpartySCIP
+(see \ref PkgSolverSummary). It should be noted that GLPK only manages to solve small problems, i.e., objects
+with reasonably simple structure.
+In case you are reconstructing more complex objects, you may need to consider more efficient open source
+solvers (e.g., <a href = "https://projects.coin-or.org/Cbc">CBC</a>) or even commercial solvers (e.g.,
+<a href = "http://www.gurobi.com/">Gurobi</a>, <a href = "https://www.ibm.com/analytics/cplex-optimizer">CPLEX</a>).
+The following table gives a rough idea of the performance of some solvers.
+
+| Model            | Problem Size <br>variables/constraints | Gurobi | CBC | SCIP | GLPK | LP_SOLVE |
+| :----:           | :----:       | :----  |     |      |      |          |
+| ![ ](model1.png) | 1244/2660    |  0.05 sec    | 0.2 sec     | 0.3 sec  | 9 sec     |  15 min   |
+| ![ ](model2.png) | 2582/5420    |  0.2 sec     | 0.6 sec     | 2.6 sec  | 35 min    |  -        |
+| ![ ](model3.png) | 21556/442191 |  2 sec       | 7.5 sec     | 9 sec    | -         |  -        |
+
+
+*/
+} /* namespace CGAL */

Polygonal_surface_reconstruction/doc/Polygonal_surface_reconstruction/dependencies

@@ -0,0 +1,8 @@
+Manual
+Kernel_23
+STL_Extension
+Circulator
+Alpha_shapes_2
+Surface_mesh
+Point_set_processing_3
+Solver_interface

Polygonal_surface_reconstruction/doc/Polygonal_surface_reconstruction/examples.txt

@@ -0,0 +1,5 @@
+/*!
+\example Polygonal_surface_reconstruction/polyfit_example_without_input_planes.cpp
+\example Polygonal_surface_reconstruction/polyfit_example_user_provided_planes.cpp
+\example Polygonal_surface_reconstruction/polyfit_example_model_complexty_control.cpp
+*/

Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/CMakeLists.txt

@@ -0,0 +1,62 @@
+# Created by the script cgal_create_CMakeLists
+# This is the CMake script for compiling a set of CGAL applications.
+
+project( Polygonal_surface_reconstruction )
+
+
+cmake_minimum_required(VERSION 2.8.11)
+
+# CGAL and its components
+find_package( CGAL QUIET COMPONENTS  )
+
+if ( NOT CGAL_FOUND )
+
+  message(STATUS "This project requires the CGAL library, and will not be compiled.")
+  return()  
+
+endif()
+
+# include helper file
+include( ${CGAL_USE_FILE} )
+
+
+# Boost and its components
+find_package( Boost REQUIRED )
+
+if ( NOT Boost_FOUND )
+
+  message(STATUS "This project requires the Boost library, and will not be compiled.")
+
+  return()  
+
+endif()
+
+find_package( SCIP )
+
+if ( NOT SCIP_FOUND )
+
+  message(STATUS "This project requires the SCIP library, and will not be compiled.")
+
+  return()  
+
+endif()
+
+
+# include for local directory
+
+# include for local package
+include_directories( BEFORE ../../include )
+include_directories( BEFORE ${SCIP_INCLUDE_DIRS} )
+
+# Creating entries for all C++ files with "main" routine
+# ##########################################################
+
+include( CGAL_CreateSingleSourceCGALProgram )
+
+create_single_source_cgal_program( "polyfit_example_without_input_planes.cpp" )
+create_single_source_cgal_program( "polyfit_example_user_provided_planes.cpp" )
+create_single_source_cgal_program( "polyfit_example_model_complexty_control.cpp" )
+
+target_link_libraries( polyfit_example_without_input_planes PRIVATE ${SCIP_LIBRARIES} )
+target_link_libraries( polyfit_example_user_provided_planes PRIVATE ${SCIP_LIBRARIES} )
+target_link_libraries( polyfit_example_model_complexty_control PRIVATE ${SCIP_LIBRARIES} )

Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_model_complexty_control.cpp

@@ -0,0 +1,133 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/IO/read_ply_points.h>
+#include <CGAL/IO/Writer_OFF.h>
+#include <CGAL/property_map.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polygonal_surface_reconstruction.h>
+#include <CGAL/SCIP_mixed_integer_program_traits.h>
+#include <CGAL/Timer.h>
+
+#include <fstream>
+
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel		Kernel;
+typedef Kernel::Point_3											Point;
+typedef Kernel::Vector_3										Vector;
+typedef	CGAL::Polygonal_surface_reconstruction<Kernel>			Polygonal_surface_reconstruction;
+typedef CGAL::Surface_mesh<Point>								Surface_mesh;
+
+typedef CGAL::SCIP_mixed_integer_program_traits<double>			MIP_Solver;
+
+// Point with normal, and plane index
+typedef boost::tuple<Point, Vector, int>						PNI;
+typedef CGAL::Nth_of_tuple_property_map<0, PNI>					Point_map;
+typedef CGAL::Nth_of_tuple_property_map<1, PNI>					Normal_map;
+typedef CGAL::Nth_of_tuple_property_map<2, PNI>					Plane_index_map;
+
+/*
+* The following example shows how to control the model complexity by
+* increasing the influence of the model complexity term.
+* In this example, the intermediate results from plane extraction and 
+* candidate generation are cached and reused.
+*/
+
+int main()
+{
+    const std::string& input_file("data/building.ply");
+	std::ifstream input_stream(input_file.c_str());
+
+	std::vector<PNI> points; // store points
+
+	std::cout << "Loading point cloud: " << input_file << "...";
+	CGAL::Timer t;
+	t.start();
+
+	if (!input_stream ||
+		!CGAL::read_ply_points_with_properties(
+			input_stream,
+			std::back_inserter(points),
+			CGAL::make_ply_point_reader(Point_map()),
+			CGAL::make_ply_normal_reader(Normal_map()),
+			std::make_pair(Plane_index_map(), CGAL::PLY_property<int>("segment_index"))))
+	{
+		std::cerr << "Error: cannot read file " << input_file << std::endl;
+		return EXIT_FAILURE;
+	}
+	else
+		std::cout << " Done. " << points.size() << " points. Time: " << t.time() << " sec." << std::endl;
+
+	//////////////////////////////////////////////////////////////////////////
+
+	std::cout << "Generating candidate faces...";
+	t.reset();
+
+	Polygonal_surface_reconstruction algo(
+		points,
+		Point_map(),
+		Normal_map(),
+		Plane_index_map()
+	);
+
+	std::cout << " Done. Time: " << t.time() << " sec." << std::endl;
+
+	//////////////////////////////////////////////////////////////////////////
+
+	// Reconstruction with complexity control
+	
+	// Model 1: more detail
+	Surface_mesh model;
+
+	std::cout << "Reconstructing with complexity 0.2...";
+	t.reset();
+	if (!algo.reconstruct<MIP_Solver>(model, 0.43, 0.27, 0.2)) {
+		std::cerr << " Failed: " << algo.error_message() << std::endl;
+		return EXIT_FAILURE;
+	}
+	else {
+       const std::string& output_file = "data/building_result-0.2.off";
+       std::ofstream output_stream(output_file.c_str());
+       if (output_stream && CGAL::write_off(output_stream, model))
+			std::cout << " Done. Saved to " << output_file << ". Time: " << t.time() << " sec." << std::endl;
+		else {
+           std::cerr << " Failed saving file." << std::endl;
+           return EXIT_FAILURE;
+       }
+	}
+
+	// Model 2: a little less detail
+	std::cout << "Reconstructing with complexity 0.6...";
+	t.reset();
+	if (!algo.reconstruct<MIP_Solver>(model, 0.43, 0.27, 0.6)) {
+		std::cerr << " Failed: " << algo.error_message() << std::endl;
+		return EXIT_FAILURE;
+	}
+	else {
+       const std::string& output_file = "data/building_result-0.6.off";
+       std::ofstream output_stream(output_file.c_str());
+       if (output_stream && CGAL::write_off(output_stream, model))
+			std::cout << " Done. Saved to " << output_file << ". Time: " << t.time() << " sec." << std::endl;
+		else {
+           std::cerr << " Failed saving file." << std::endl;
+           return EXIT_FAILURE;
+       }
+	}
+
+	// Model 3: more less detail
+	std::cout << "Reconstructing with complexity 1.5...";
+	t.reset();
+	if (!algo.reconstruct<MIP_Solver>(model, 0.43, 0.27, 1.5)) {
+		std::cerr << " Failed: " << algo.error_message() << std::endl;
+		return EXIT_FAILURE;
+	}
+	else {
+		const std::string& output_file = "data/building_result-1.5.off";
+		std::ofstream output_stream(output_file.c_str());
+		if (output_stream && CGAL::write_off(output_stream, model))
+			std::cout << " Done. Saved to " << output_file << ". Time: " << t.time() << " sec." << std::endl;
+		else {
+			std::cerr << " Failed saving file." << std::endl;
+			return EXIT_FAILURE;
+		}
+	}
+	return EXIT_SUCCESS;
+}

Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_user_provided_planes.cpp

@@ -0,0 +1,96 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/IO/Writer_OFF.h>
+#include <CGAL/IO/read_ply_points.h>
+#include <CGAL/property_map.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polygonal_surface_reconstruction.h>
+#include <CGAL/SCIP_mixed_integer_program_traits.h>
+#include <CGAL/Timer.h>
+
+#include <fstream>
+
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel		Kernel;
+typedef Kernel::Point_3											Point;
+typedef Kernel::Vector_3										Vector;
+typedef	CGAL::Polygonal_surface_reconstruction<Kernel>			Polygonal_surface_reconstruction;
+typedef CGAL::Surface_mesh<Point>								Surface_mesh;
+
+// Point with normal, and plane index
+typedef boost::tuple<Point, Vector, int>						PNI;
+typedef CGAL::Nth_of_tuple_property_map<0, PNI>					Point_map;
+typedef CGAL::Nth_of_tuple_property_map<1, PNI>					Normal_map;
+typedef CGAL::Nth_of_tuple_property_map<2, PNI>					Plane_index_map;
+
+typedef CGAL::SCIP_mixed_integer_program_traits<double>			MIP_Solver;
+
+/*
+* The following example shows the reconstruction using user-provided
+* planar segments stored in PLY format. In the PLY format, a property
+* named "segment_index" stores the plane index for each point (-1 if 
+* the point is not assigned to a plane). 
+*/
+
+int main()
+{
+    const std::string& input_file("data/ball.ply");
+	std::ifstream input_stream(input_file.c_str());
+
+	std::vector<PNI> points; // store points
+
+	std::cout << "Loading point cloud: " << input_file << "...";
+	CGAL::Timer t;
+	t.start();
+
+	if (!input_stream ||
+		!CGAL::read_ply_points_with_properties(
+			input_stream,
+			std::back_inserter(points),
+			CGAL::make_ply_point_reader(Point_map()),
+			CGAL::make_ply_normal_reader(Normal_map()),
+			std::make_pair(Plane_index_map(), CGAL::PLY_property<int>("segment_index"))))
+	{
+		std::cerr << "Error: cannot read file " << input_file << std::endl;
+		return EXIT_FAILURE;
+	}
+	else
+		std::cout << " Done. " << points.size() << " points. Time: " << t.time() << " sec." << std::endl;
+
+	//////////////////////////////////////////////////////////////////////////
+
+	std::cout << "Generating candidate faces...";
+	t.reset();
+
+	Polygonal_surface_reconstruction algo(
+		points, 
+		Point_map(), 
+		Normal_map(),
+		Plane_index_map()
+	);
+
+	std::cout << " Done. Time: " << t.time() << " sec." << std::endl;
+
+	//////////////////////////////////////////////////////////////////////////
+
+	Surface_mesh model;
+
+	std::cout << "Reconstructing...";
+	t.reset();
+
+	if (!algo.reconstruct<MIP_Solver>(model)) {
+		std::cerr << " Failed: " << algo.error_message() << std::endl;
+		return EXIT_FAILURE;
+	}
+
+	// Saves the mesh model
+    const std::string& output_file("data/ball_result.off");
+    std::ofstream output_stream(output_file.c_str());
+    if (output_stream && CGAL::write_off(output_stream, model))
+		std::cout << " Done. Saved to " << output_file << ". Time: " << t.time() << " sec." << std::endl;
+	else {
+        std::cerr << " Failed saving file." << std::endl;
+        return EXIT_FAILURE;
+    }
+
+	return EXIT_SUCCESS;
+}

Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_without_input_planes.cpp

@@ -0,0 +1,136 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/IO/read_xyz_points.h>
+#include <CGAL/IO/Writer_OFF.h>
+#include <CGAL/property_map.h>
+#include <CGAL/Shape_detection_3.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polygonal_surface_reconstruction.h>
+#include <CGAL/SCIP_mixed_integer_program_traits.h>
+#include <CGAL/Timer.h>
+
+#include <fstream>
+
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel			Kernel;
+
+typedef Kernel::Point_3												Point;
+typedef Kernel::Vector_3											Vector;
+
+// Point with normal, and plane index
+typedef boost::tuple<Point, Vector, int>							PNI;
+typedef std::vector<PNI>											Point_vector;
+typedef CGAL::Nth_of_tuple_property_map<0, PNI>						Point_map;
+typedef CGAL::Nth_of_tuple_property_map<1, PNI>						Normal_map;
+typedef CGAL::Nth_of_tuple_property_map<2, PNI>						Plane_index_map;
+
+typedef CGAL::Shape_detection_3::Shape_detection_traits<Kernel, Point_vector, Point_map, Normal_map>	Traits;
+typedef CGAL::Shape_detection_3::Efficient_RANSAC<Traits>			Efficient_ransac;
+typedef CGAL::Shape_detection_3::Plane<Traits>						Plane;
+typedef CGAL::Shape_detection_3::Point_to_shape_index_map<Traits>	Point_to_shape_index_map;
+
+typedef	CGAL::Polygonal_surface_reconstruction<Kernel>				Polygonal_surface_reconstruction;
+typedef CGAL::Surface_mesh<Point>									Surface_mesh;
+
+typedef CGAL::SCIP_mixed_integer_program_traits<double>				MIP_Solver;
+
+/*
+* This example first extracts planes from the input point cloud
+* (using RANSAC with default parameters) and then reconstructs 
+* the surface model from the planes.
+*/
+
+int main()
+{
+	Point_vector points;
+
+	// Loads point set from a file. 
+	const std::string input_file("data/cube.pwn");
+    std::ifstream input_stream(input_file.c_str());
+	if (input_stream.fail()) {
+		std::cerr << "failed open file \'" <<input_file << "\'" << std::endl;
+		return EXIT_FAILURE;
+	}
+	std::cout << "Loading point cloud: " << input_file << "...";
+
+	CGAL::Timer t;
+	t.start();
+    if (!input_stream ||
+		!CGAL::read_xyz_points(input_stream,
+			std::back_inserter(points),
+			CGAL::parameters::point_map(Point_map()).normal_map(Normal_map())))
+	{
+		std::cerr << "Error: cannot read file " << input_file << std::endl;
+		return EXIT_FAILURE;
+	}
+	else
+		std::cout << " Done. " << points.size() << " points. Time: " << t.time() << " sec." << std::endl;
+
+	// Shape detection
+	Efficient_ransac ransac;
+	ransac.set_input(points);
+	ransac.add_shape_factory<Plane>();
+
+	std::cout << "Extracting planes...";
+	t.reset();
+	ransac.detect();
+
+	Efficient_ransac::Plane_range planes = ransac.planes();
+	std::size_t num_planes = planes.size();
+
+	std::cout << " Done. " << num_planes << " planes extracted. Time: " << t.time() << " sec." << std::endl;
+
+	// Stores the plane index of each point as the third element of the tuple.
+	Point_to_shape_index_map shape_index_map(points, planes);
+	for (std::size_t i = 0; i < points.size(); ++i) {
+		// Uses the get function from the property map that accesses the 3rd element of the tuple.
+		int plane_index = get(shape_index_map, i);
+		points[i].get<2>() = plane_index;
+	}
+
+	//////////////////////////////////////////////////////////////////////////
+
+	std::cout << "Generating candidate faces...";
+	t.reset();
+
+	Polygonal_surface_reconstruction algo(
+		points,
+		Point_map(),
+		Normal_map(),
+		Plane_index_map()
+	);
+
+	std::cout << " Done. Time: " << t.time() << " sec." << std::endl;
+
+	//////////////////////////////////////////////////////////////////////////
+
+	Surface_mesh model;
+
+	std::cout << "Reconstructing...";
+	t.reset();
+
+	if (!algo.reconstruct<MIP_Solver>(model)) {
+		std::cerr << " Failed: " << algo.error_message() << std::endl;
+		return EXIT_FAILURE;
+	}
+
+	const std::string& output_file("data/cube_result.off");
+	std::ofstream output_stream(output_file.c_str());
+	if (output_stream && CGAL::write_off(output_stream, model))
+		std::cout << " Done. Saved to " << output_file << ". Time: " << t.time() << " sec." << std::endl;
+	else {
+		std::cerr << " Failed saving file." << std::endl;
+		return EXIT_FAILURE;
+	}
+
+	//////////////////////////////////////////////////////////////////////////
+
+	// Also stores the candidate faces as a surface mesh to a file
+	Surface_mesh candidate_faces;
+	algo.output_candidate_faces(candidate_faces);
+	const std::string& candidate_faces_file("data/cube_candidate_faces.off");
+	std::ofstream candidate_stream(candidate_faces_file.c_str());
+	if (candidate_stream && CGAL::write_off(candidate_stream, candidate_faces))
+		std::cout << "Candidate faces saved to " << candidate_faces_file << "." << std::endl;
+
+	return EXIT_SUCCESS;
+}

Polygonal_surface_reconstruction/include/CGAL/Polygonal_surface_reconstruction.h

@@ -0,0 +1,480 @@
+// Copyright (c) 2018  Liangliang Nan
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// 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.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0+
+//
+// Author(s) : Liangliang Nan
+
+#ifndef CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_H
+#define CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_H
+
+#include <CGAL/license/Polygonal_surface_reconstruction.h>
+
+#include <CGAL/bounding_box.h>
+#include <CGAL/property_map.h>
+#include <CGAL/internal/hypothesis.h>
+#include <CGAL/internal/compute_confidences.h>
+#include <CGAL/internal/point_set_with_planes.h>
+
+#include <unordered_map>
+
+/*!
+\file Polygonal_surface_reconstruction.h
+*/
+
+namespace CGAL {
+
+	/*!
+	\ingroup PkgPolygonalSurfaceReconstruction
+
+	\brief
+
+	Implementation of the Polygonal Surface Reconstruction method.
+
+	Given a set of 3D points with unoriented normals sampled
+	from the outer boundary of a piecewise planar object, the Polygonal Surface
+	Reconstruction method \cgalCite{nan2017polyfit} outputs a simplified and
+	watertight surface mesh interpolating the input point set.
+
+	The method first generates a set of face candidates by intersecting the planar
+	primitives. Then an optimal subset of the candidate faces is selected through
+	optimization under hard constraints that enforce the final model to be manifold
+	and watertight.
+
+	The reconstruction assumes the planar segmentation of the point cloud is
+	provided in the input.
+
+	\tparam GeomTraits a geometric traits class, model of Kernel
+	*/
+	template <class GeomTraits>
+	class Polygonal_surface_reconstruction
+	{
+	public:
+
+		/// \name Types
+
+		typedef typename GeomTraits::FT				FT;			///< number type.
+		typedef typename GeomTraits::Point_3		Point;		///< point type.
+		typedef typename GeomTraits::Vector_3		Vector;		///< vector type.
+		typedef typename GeomTraits::Plane_3		Plane;		///< plane type.
+
+	private:
+
+		// Polygon mesh storing the candidate faces
+		typedef CGAL::Surface_mesh<Point>				Polygon_mesh;
+
+		typedef typename Polygon_mesh::Face_index		Face_descriptor;
+		typedef typename Polygon_mesh::Vertex_index		Vertex_descriptor;
+		typedef typename Polygon_mesh::Edge_index		Edge_descriptor;
+		typedef typename Polygon_mesh::Halfedge_index	Halfedge_descriptor;
+
+		// Public methods
+	public:
+
+		/// \name Creation
+
+		/*!
+		Creates a Polygonal Surface Reconstruction object.
+		After construction, candidate faces are generated and point/face confidence values are
+		computed, allowing to reuse them in the subsequent reconstruction step with different parameters.
+
+		\tparam PointRange is the range of input points, model of `ConstRange`.
+		\tparam PointMap is a model of `ReadablePropertyMap` with value	type `GeomTraits::Point_3`.
+		\tparam NormalMap is a model of `ReadablePropertyMap` with value type `GeomTraits::Vector_3`.
+		\tparam IndexMap is a model of `ReadablePropertyMap` with value	type `int`.
+
+		\param points range of input points.
+		\param point_map property map: value_type of `typename PointRange::const_iterator` -> `Point_3`
+		\param normal_map property map: value_type of `typename PointRange::const_iterator` -> `Vector_3`
+		\param index_map property map: value_type of `typename PointRange::const_iterator` -> `int`,
+		denoting the index of the plane it belongs to (-1 if the point is not assigned to a plane)
+		*/
+		template <
+			typename PointRange,
+			typename PointMap,
+			typename NormalMap,
+			typename IndexMap
+		>
+			Polygonal_surface_reconstruction(
+				const PointRange& points,
+				PointMap point_map,
+				NormalMap normal_map,
+				IndexMap index_map
+			);
+
+
+		/// \name Operations
+
+		/** Reconstructs a watertight polygonal mesh model.
+
+		\tparam MixedIntegerProgramTraits a model of `MixedIntegerProgramTraits`
+
+		\tparam PolygonMesh a model of `MutableFaceGraph`
+
+		\return `true` if the reconstruction succeeded, `false` otherwise.
+		*/
+		template <typename MixedIntegerProgramTraits, typename PolygonMesh>
+		bool reconstruct(
+			PolygonMesh& output_mesh,		///< the final reconstruction result
+			double wt_fitting = 0.43,		///< weight for the data fitting term.
+			double wt_coverage = 0.27,		///< weight for the point coverage term.
+			double wt_complexity = 0.30		///< weight for the model complexity term.
+		);
+
+		/*! Gives the user the possibility to access the intermediate candidate faces
+			(i.e., the faces induced by the intersection of the supporting planes).
+		\tparam PolygonMesh a model of `MutableFaceGraph`.
+		*/
+		template <typename PolygonMesh>
+		void output_candidate_faces(PolygonMesh& candidate_faces) const;
+
+		/// Gets the error message (if reconstruction failed).
+		const std::string& error_message() const { return error_message_; }
+
+		// Data members.
+	private:
+		internal::Hypothesis<GeomTraits> hypothesis_;
+
+		// The generated candidate faces stored as a polygon mesh
+		Polygon_mesh	candidate_faces_;
+
+		std::string		error_message_;
+
+	private: // Copying is not allowed
+		Polygonal_surface_reconstruction(const Polygonal_surface_reconstruction& psr);
+
+	}; // end of Polygonal_surface_reconstruction
+
+
+	//////////////////////////////////////////////////////////////////////////s
+
+	// implementations
+
+	template <class GeomTraits>
+
+	template <
+		typename PointRange,
+		typename PointMap,
+		typename NormalMap,
+		typename IndexMap
+	>
+		Polygonal_surface_reconstruction<GeomTraits>::Polygonal_surface_reconstruction(
+			const PointRange& points,
+			PointMap point_map,
+			NormalMap normal_map,
+			IndexMap index_map
+		) : error_message_("")
+	{
+		if (points.empty()) {
+			error_message_ = "empty input points";
+			return;
+		}
+
+		typedef internal::Planar_segment<GeomTraits>			Planar_segment;
+		typedef internal::Point_set_with_planes<GeomTraits>		Point_set_with_planes;
+
+		Point_set_with_planes point_set(points, point_map, normal_map, index_map);
+
+		const std::vector< Planar_segment* >& planar_segments = point_set.planar_segments();
+		if (planar_segments.size() < 4) {
+			error_message_ = "at least 4 planes required to reconstruct a closed surface mesh (only "
+				+ std::to_string(planar_segments.size()) + " provided)";
+			return;
+		}
+
+		hypothesis_.generate(point_set, candidate_faces_);
+
+		typedef internal::Candidate_confidences<GeomTraits>		Candidate_confidences;
+		Candidate_confidences conf;
+		conf.compute(point_set, candidate_faces_);
+	}
+
+
+	template <class GeomTraits>
+	template <typename PolygonMesh>
+	void Polygonal_surface_reconstruction<GeomTraits>::output_candidate_faces(PolygonMesh& candidate_faces) const {
+		candidate_faces.clear();	// make sure it is empty.
+		CGAL::copy_face_graph(candidate_faces_, candidate_faces);
+	}
+
+
+	template <class GeomTraits>
+	template <typename MixedIntegerProgramTraits, typename PolygonMesh>
+	bool Polygonal_surface_reconstruction<GeomTraits>::reconstruct(
+		PolygonMesh& output_mesh,
+		double wt_fitting /* = 0.43 */,
+		double wt_coverage /* = 0.27 */,
+		double wt_complexity /* = 0.30 */)
+	{
+		if (!error_message_.empty()) { // an error has occurred in the constructor
+			return false;
+		}
+
+		if (candidate_faces_.num_faces() < 4) {
+			error_message_ = "at least 4 candidate faces required to reconstruct a closed surface mesh (only "
+				+ std::to_string(candidate_faces_.num_faces()) + " computed)";
+			return false;
+		}
+
+		typedef typename internal::Hypothesis<GeomTraits>::Adjacency Adjacency;
+		const Adjacency& adjacency = hypothesis_.extract_adjacency(candidate_faces_);
+
+		// Internal data structure
+		Polygon_mesh target_mesh = candidate_faces_;
+
+		// The number of supporting points of each face
+		typename Polygon_mesh::template Property_map<Face_descriptor, std::size_t> face_num_supporting_points =
+			target_mesh.template add_property_map<Face_descriptor, std::size_t>("f:num_supporting_points").first;
+
+		// The area of each face
+		typename Polygon_mesh::template Property_map<Face_descriptor, FT> face_areas =
+			target_mesh.template add_property_map<Face_descriptor, FT>("f:face_area").first;
+
+		// The point covered area of each face
+		typename Polygon_mesh::template Property_map<Face_descriptor, FT> face_covered_areas =
+			target_mesh.template add_property_map<Face_descriptor, FT>("f:covered_area").first;
+
+		// The supporting plane of each face
+		typename Polygon_mesh::template Property_map<Face_descriptor, const Plane*> face_supporting_planes =
+			target_mesh.template add_property_map<Face_descriptor, const Plane*>("f:supp_plane").first;
+
+		// Gives each face an index
+		typename Polygon_mesh::template Property_map<Face_descriptor, std::size_t> face_indices =
+			target_mesh.template add_property_map<Face_descriptor, std::size_t>("f:index").first;
+
+		double total_points = 0.0;
+		std::size_t idx = 0;
+		BOOST_FOREACH(Face_descriptor f, target_mesh.faces()) {
+			total_points += face_num_supporting_points[f];
+			face_indices[f] = idx;
+			++idx;
+		}
+
+
+		typedef MixedIntegerProgramTraits								MIP_Solver;
+		typedef typename MixedIntegerProgramTraits::Variable			Variable;
+		typedef typename MixedIntegerProgramTraits::Linear_objective	Linear_objective;
+		typedef typename MixedIntegerProgramTraits::Linear_constraint	Linear_constraint;
+
+		MIP_Solver solver;
+
+		// Adds variables
+
+		// Binary variables:
+		// x[0] ... x[num_faces - 1] : binary labels of all the input faces
+		// x[num_faces] ... x[num_faces + num_edges - 1] : binary labels of all the intersecting edges (remain or not)
+		// x[num_faces + num_edges] ... x[num_faces + num_edges + num_edges] : binary labels of corner edges (sharp edge of not)
+
+		std::size_t num_faces = target_mesh.number_of_faces();
+		std::size_t num_edges(0);
+
+		typedef typename internal::Hypothesis<GeomTraits>::Intersection	Intersection;
+
+		std::unordered_map<const Intersection*, std::size_t> edge_usage_status;	// keep or remove an intersecting edges
+		for (std::size_t i = 0; i < adjacency.size(); ++i) {
+			const Intersection& fan = adjacency[i];
+			if (fan.size() == 4) {
+				std::size_t var_idx = num_faces + num_edges;
+				edge_usage_status[&fan] = var_idx;
+				++num_edges;
+			}
+		}
+
+		std::size_t total_variables = num_faces + num_edges + num_edges;
+
+		const std::vector<Variable*>& variables = solver.create_variables(total_variables);
+		for (std::size_t i = 0; i < total_variables; ++i) {
+			Variable* v = variables[i];
+			v->set_variable_type(Variable::BINARY);
+		}
+
+		// Adds objective
+
+		const typename Polygon_mesh::template Property_map<Vertex_descriptor, Point>& coords = target_mesh.points();
+		std::vector<Point> vertices(target_mesh.number_of_vertices());
+		idx = 0;
+		BOOST_FOREACH(Vertex_descriptor v, target_mesh.vertices()) {
+			vertices[idx] = coords[v];
+			++idx;
+		}
+
+		typedef typename CGAL::Iso_cuboid_3<GeomTraits>		Box;
+
+		const Box& box = CGAL::bounding_box(vertices.begin(), vertices.end());
+		FT dx = box.xmax() - box.xmin();
+		FT dy = box.ymax() - box.ymin();
+		FT dz = box.zmax() - box.zmin();
+		FT box_area = FT(2.0) * (dx * dy + dy * dz + dz * dx);
+
+		// Chooses a better scale: all actual values multiplied by total number of points
+		double coeff_data_fitting = wt_fitting;
+		double coeff_coverage = total_points * wt_coverage / box_area;
+		double coeff_complexity = total_points * wt_complexity / double(adjacency.size());
+
+		Linear_objective * objective = solver.create_objective(Linear_objective::MINIMIZE);
+
+		std::unordered_map<const Intersection*, std::size_t> edge_sharp_status;	// the edge is sharp or not
+		std::size_t num_sharp_edges = 0;
+		for (std::size_t i = 0; i < adjacency.size(); ++i) {
+			const Intersection& fan = adjacency[i];
+			if (fan.size() == 4) {
+				std::size_t var_idx = num_faces + num_edges + num_sharp_edges;
+				edge_sharp_status[&fan] = var_idx;
+
+				// Accumulates model complexity term
+				objective->add_coefficient(variables[var_idx], coeff_complexity);
+				++num_sharp_edges;
+			}
+		}
+		CGAL_assertion(num_edges == num_sharp_edges);
+
+		BOOST_FOREACH(Face_descriptor f, target_mesh.faces()) {
+			std::size_t var_idx = face_indices[f];
+
+			// Accumulates data fitting term
+			std::size_t num = face_num_supporting_points[f];
+			objective->add_coefficient(variables[var_idx], -coeff_data_fitting * num);
+
+			// Accumulates model coverage term
+			double uncovered_area = (face_areas[f] - face_covered_areas[f]);
+			objective->add_coefficient(variables[var_idx], coeff_coverage * uncovered_area);
+		}
+
+		// Adds constraints: the number of faces associated with an edge must be either 2 or 0
+		std::size_t var_edge_used_idx = 0;
+		for (std::size_t i = 0; i < adjacency.size(); ++i) {
+			Linear_constraint* c = solver.create_constraint(0.0, 0.0);
+			const Intersection& fan = adjacency[i];
+			for (std::size_t j = 0; j < fan.size(); ++j) {
+				Face_descriptor f = target_mesh.face(fan[j]);
+				std::size_t var_idx = face_indices[f];
+				c->add_coefficient(variables[var_idx], 1.0);
+			}
+
+			if (fan.size() == 4) {
+				std::size_t var_idx = num_faces + var_edge_used_idx;
+				c->add_coefficient(variables[var_idx], -2.0);  //
+				++var_edge_used_idx;
+			}
+			else { // boundary edge
+				   // will be set to 0 (i.e., we don't allow open surface)
+			}
+		}
+
+		// Adds constraints: for the sharp edges. The explanation of posing this constraint can be found here:
+		// https://user-images.githubusercontent.com/15526536/30185644-12085a9c-942b-11e7-831d-290dd2a4d50c.png
+		double M = 1.0;
+		for (std::size_t i = 0; i < adjacency.size(); ++i) {
+			const Intersection& fan = adjacency[i];
+			if (fan.size() != 4)
+				continue;
+
+			// If an edge is sharp, the edge must be selected first:
+			// X[var_edge_usage_idx] >= X[var_edge_sharp_idx]
+			Linear_constraint* c = solver.create_constraint(0.0);
+			std::size_t var_edge_usage_idx = edge_usage_status[&fan];
+			c->add_coefficient(variables[var_edge_usage_idx], 1.0);
+			std::size_t var_edge_sharp_idx = edge_sharp_status[&fan];
+			c->add_coefficient(variables[var_edge_sharp_idx], -1.0);
+
+			for (std::size_t j = 0; j < fan.size(); ++j) {
+				Face_descriptor f1 = target_mesh.face(fan[j]);
+				const Plane* plane1 = face_supporting_planes[f1];
+				std::size_t fid1 = face_indices[f1];
+				for (std::size_t k = j + 1; k < fan.size(); ++k) {
+					Face_descriptor f2 = target_mesh.face(fan[k]);
+					const Plane* plane2 = face_supporting_planes[f2];
+					std::size_t fid2 = face_indices[f2];
+
+					if (plane1 != plane2) {
+						// The constraint is:
+						//X[var_edge_sharp_idx] + M * (3 - (X[fid1] + X[fid2] + X[var_edge_usage_idx])) >= 1
+						// which equals to
+						//X[var_edge_sharp_idx] - M * X[fid1] - M * X[fid2] - M * X[var_edge_usage_idx] >= 1 - 3M
+						c = solver.create_constraint(1.0 - 3.0 * M);
+						c->add_coefficient(variables[var_edge_sharp_idx], 1.0);
+						c->add_coefficient(variables[fid1], -M);
+						c->add_coefficient(variables[fid2], -M);
+						c->add_coefficient(variables[var_edge_usage_idx], -M);
+					}
+				}
+			}
+		}
+
+		// Optimization
+
+		if (solver.solve()) {
+
+			// Marks results
+			const std::vector<double>& X = solver.solution();
+
+			std::vector<Face_descriptor> to_delete;
+			std::size_t f_idx(0);
+			BOOST_FOREACH(Face_descriptor f, target_mesh.faces()) {
+				if (static_cast<int>(std::round(X[f_idx])) == 0)
+					to_delete.push_back(f);
+				++f_idx;
+			}
+
+			for (std::size_t i = 0; i < to_delete.size(); ++i) {
+				Face_descriptor f = to_delete[i];
+				Halfedge_descriptor h = target_mesh.halfedge(f);
+				Euler::remove_face(h, target_mesh);
+			}
+
+			// Marks the sharp edges
+			typename Polygon_mesh::template Property_map<Edge_descriptor, bool> edge_is_sharp =
+				target_mesh.template add_property_map<Edge_descriptor, bool>("e:sharp_edges").first;
+			BOOST_FOREACH(Edge_descriptor e, target_mesh.edges())
+				edge_is_sharp[e] = false;
+
+			for (std::size_t i = 0; i < adjacency.size(); ++i) {
+				const Intersection& fan = adjacency[i];
+				if (fan.size() != 4)
+					continue;
+
+				std::size_t idx_sharp_var = edge_sharp_status[&fan];
+				if (static_cast<int>(X[idx_sharp_var]) == 1) {
+					for (std::size_t j = 0; j < fan.size(); ++j) {
+						Halfedge_descriptor h = fan[j];
+						Face_descriptor f = target_mesh.face(h);
+						if (f != Polygon_mesh::null_face()) { // some faces may be deleted
+							std::size_t fid = face_indices[f];
+							if (static_cast<int>(std::round(X[fid])) == 1) {
+								Edge_descriptor e = target_mesh.edge(h);
+								edge_is_sharp[e] = true;
+								break;
+							}
+						}
+					}
+				}
+			}
+
+			// Converts from internal data structure to the required `PolygonMesh`.
+			output_mesh.clear();	// make sure it is empty.
+			CGAL::copy_face_graph(target_mesh, output_mesh);
+		}
+		else {
+			error_message_ = "solving the binary program failed";
+			return false;
+		}
+
+		return true;
+	}
+
+} //namespace CGAL
+
+#endif // CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_H

Polygonal_surface_reconstruction/include/CGAL/internal/alpha_shape_mesh.h

@@ -0,0 +1,241 @@
+// Copyright(c) 2018  Liangliang Nan
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// 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.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0+
+//
+// Author(s) : Liangliang Nan
+
+#ifndef CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_ALPHA_SHAPE_MESH_H
+#define CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_ALPHA_SHAPE_MESH_H
+
+#include <CGAL/license/Polygonal_surface_reconstruction.h>
+
+#include <CGAL/Alpha_shape_2.h>
+#include <CGAL/Delaunay_triangulation_2.h>
+#include <CGAL/Triangulation_hierarchy_2.h>
+#include <CGAL/Surface_mesh.h>
+
+
+// warn the user if he/she is using an unsupported version of CGAL
+#if CGAL_VERSION_NR < 1041100000
+#error CGAL 4.11 or above is required (due to the breaking change in CGAL 4.11). Please update your code.
+#endif
+
+/*!
+\file alpha_shape_mesh.h
+*/
+
+namespace CGAL {
+
+	namespace internal {
+
+		/// \cond SKIP_IN_MANUA
+
+		template <typename Ht>
+		class Alpha_shape;
+
+		/// \endcond
+
+		/// \cond SKIP_IN_MANUA
+
+		/* A vertex class with an additional member representing its index */
+		template < class Gt, class VB = CGAL::Triangulation_hierarchy_vertex_base_2<Gt> >
+		class AS_vertex_base : public VB
+		{
+		public:
+			typedef VB								Base;
+			typedef typename VB::Vertex_handle      Vertex_handle;
+			typedef typename VB::Face_handle        Face_handle;
+			typedef typename VB::Point              Point;
+
+			template < typename TDS2 >
+			struct Rebind_TDS {
+				typedef typename VB::template Rebind_TDS<TDS2>::Other	VB2;
+				typedef AS_vertex_base<Gt, VB2>                         Other;
+			};
+
+		public:
+			AS_vertex_base() : Base(), index_(-1) {}
+			AS_vertex_base(const Point & p) : Base(p), index_(-1) {}
+			AS_vertex_base(const Point & p, Face_handle f) : Base(f, p), index_(-1) {}
+			AS_vertex_base(Face_handle f) : Base(f), index_(-1) {}
+
+			void set_index(int idx) { index_ = idx; }
+			int  index() const { return index_; }
+
+		private:
+			int index_;
+		};
+
+
+		template <typename Ht>
+		class Alpha_shape : public Alpha_shape_2<Ht>
+		{
+		public:
+			typedef Alpha_shape_2<Ht>						Parent_class;
+			typedef typename Ht::Point_2					Point2;
+			typedef typename Parent_class::Vertex_handle	Vertex_handle;
+
+		public:
+			// constructs alpha shapes from the input points
+			template <typename InputIterator>
+			Alpha_shape(InputIterator first, InputIterator beyond);
+		};
+
+		/// \endcond
+
+
+		/**
+		*	An Alpha Shape Mesh approximates the point covered region by a mesh representation.
+		*/
+
+		template <typename Kernel>
+		class Alpha_shape_mesh
+		{
+			typedef typename Kernel::FT						FT;
+			typedef typename Kernel::Point_2				Point2;
+			typedef typename Kernel::Point_3				Point3;
+			typedef typename Kernel::Plane_3				Plane3;
+			typedef CGAL::Surface_mesh<Point3>				Mesh3;
+			typedef typename Mesh3::Vertex_index			Vertex_descriptor;
+
+			typedef CGAL::Alpha_shape_vertex_base_2<Kernel>			Avb;
+			typedef AS_vertex_base<Avb>								Av;
+			typedef CGAL::Triangulation_face_base_2<Kernel>			Tf;
+			typedef CGAL::Alpha_shape_face_base_2<Kernel, Tf>		Af;
+			typedef CGAL::Triangulation_default_data_structure_2<Kernel, Av, Af> Tds;
+			typedef CGAL::Delaunay_triangulation_2<Kernel, Tds>		Dt;
+			typedef CGAL::Triangulation_hierarchy_2<Dt>				Ht;
+
+		public:
+			/// Given a set of 3D points lying on 'plane', constructs alpha shapes from the
+			/// the projection of the points onto 'plane'
+			template <typename InputIterator>
+			Alpha_shape_mesh(InputIterator first, InputIterator beyond, const Plane3& plane);
+
+			~Alpha_shape_mesh() { delete alpha_shape_; }
+
+			/// Extracts the 3D mesh representation of the alpha shapes
+			bool extract_mesh(FT alpha_value, Mesh3& mesh);
+
+		private:
+			Alpha_shape<Ht>*			alpha_shape_;
+			std::vector<const Point3*>  original_points_;
+		};
+
+
+		//////////////////////////////////////////////////////////////////////////
+
+		// implementation
+
+
+		template <typename Traits>
+		template <typename InputIterator>
+		Alpha_shape<Traits>::Alpha_shape(InputIterator first, InputIterator beyond) {
+			InputIterator it = first;
+			for (int id = 0; it != beyond; ++it, ++id) {
+				const Point2& p = *it;
+				Vertex_handle vh = Traits::insert(p);
+				if (vh->index() == -1)
+					vh->set_index(id);
+				else {
+					// p was not inserted (there might be a duplicated point)
+				}
+			}
+
+			if (Parent_class::dimension() == 2) {
+				// Computes the associated _interval_face_map
+				Parent_class::initialize_interval_face_map();
+
+				// Computes the associated _interval_edge_map
+				Parent_class::initialize_interval_edge_map();
+
+				// Computes the associated _interval_vertex_map
+				Parent_class::initialize_interval_vertex_map();
+
+				// Merges the two maps
+				Parent_class::initialize_alpha_spectrum();
+			}
+		}
+
+
+		template <typename Kernel>
+		template <typename InputIterator>
+		Alpha_shape_mesh<Kernel>::Alpha_shape_mesh(InputIterator first, InputIterator beyond, const Plane3& plane) {
+			original_points_.clear();
+
+			std::vector<Point2> pts;
+			for (InputIterator it = first; it != beyond; ++it) {
+				const Point3& p = *it;
+				const Point2& q = plane.to_2d(p);
+				pts.push_back(q);
+				original_points_.push_back(&p);
+			}
+			alpha_shape_ = new Alpha_shape<Ht>(pts.begin(), pts.end());
+		}
+
+
+		template <typename Kernel>
+		bool Alpha_shape_mesh<Kernel>::extract_mesh(FT alpha_value, Mesh3& mesh) {
+			alpha_shape_->set_alpha(alpha_value);
+
+			typedef std::vector<std::size_t> Triangle;
+			std::vector<Triangle>	faces;
+
+			typedef Alpha_shape<Ht>	Alpha_shape;
+
+			typename Alpha_shape::Finite_faces_iterator fit = alpha_shape_->finite_faces_begin();
+			for (; fit != alpha_shape_->finite_faces_end(); ++fit) {
+				if (alpha_shape_->classify(fit) == Alpha_shape::INTERIOR) {
+					Triangle tri;
+					for (int i = 0; i < 3; ++i) {
+						typename Alpha_shape::Vertex_handle vh = fit->vertex(i);
+						int idx = vh->index();
+						tri.push_back(idx);
+					}
+					faces.push_back(tri);
+				}
+			}
+
+			if (faces.empty())
+				return false;
+
+			mesh.clear();
+
+			std::vector<Vertex_descriptor> descriptors(original_points_.size());
+			for (std::size_t i = 0; i < original_points_.size(); ++i) {
+				const Point3* p = original_points_[i];
+				descriptors[i] = mesh.add_vertex(*p);
+			}
+
+			for (std::size_t i = 0; i < faces.size(); ++i) {
+				std::vector<Vertex_descriptor> face;
+				const Triangle& tri = faces[i];
+				for (std::size_t j = 0; j < tri.size(); ++j) {
+					std::size_t idx = tri[j];
+					face.push_back(descriptors[idx]);
+				}
+				mesh.add_face(face);;
+			}
+
+			return true;
+		}
+
+	} //namespace internal
+
+} //namespace CGAL
+
+#endif	// CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_ALPHA_SHAPE_MESH_H

Polygonal_surface_reconstruction/include/CGAL/internal/compute_confidences.h

@@ -0,0 +1,265 @@
+// Copyright (c) 2018  Liangliang Nan
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// 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.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0+
+//
+// Author(s) : Liangliang Nan
+
+#ifndef CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_CANDIDATE_CONFIDENCES_H
+#define CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_CANDIDATE_CONFIDENCES_H
+
+#include <CGAL/license/Polygonal_surface_reconstruction.h>
+
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/compute_average_spacing.h>
+#include <CGAL/Polygon_2.h>
+#include <CGAL/assertions.h>
+#include <CGAL/internal/point_set_with_planes.h>
+#include <CGAL/internal/alpha_shape_mesh.h>
+#include <CGAL/internal/hypothesis.h>
+#include <CGAL/internal/parameters.h>
+
+
+// Concurrency
+#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_tag Concurrency_tag;
+#else
+typedef CGAL::Sequential_tag Concurrency_tag;
+#endif
+
+
+namespace CGAL {
+
+
+	namespace internal {
+
+		/**
+		*	Computes the confidences of the candidate faces.
+		*/
+
+		template <typename Kernel>
+		class Candidate_confidences
+		{
+		private:
+			typedef typename Kernel::FT						FT;
+			typedef typename Kernel::Point_3				Point;
+			typedef typename Kernel::Point_2				Point2;
+			typedef typename Kernel::Vector_3				Vector;
+			typedef typename Kernel::Line_3					Line;
+			typedef typename Kernel::Segment_3				Segment;
+			typedef typename Kernel::Plane_3				Plane;
+			typedef CGAL::Polygon_2<Kernel>					Polygon;
+			typedef Planar_segment<Kernel>					Planar_segment;
+			typedef Point_set_with_planes<Kernel>			Point_set;
+			typedef CGAL::Surface_mesh<Point>				Polygon_mesh;
+			typedef typename Polygon_mesh::Face_index		Face_descriptor;
+			typedef typename Polygon_mesh::Edge_index		Edge_descriptor;
+			typedef typename Polygon_mesh::Vertex_index		Vertex_descriptor;
+			typedef typename Polygon_mesh::Halfedge_index	Halfedge_descriptor;
+
+		public:
+			Candidate_confidences() {}
+			~Candidate_confidences() {}
+
+			/// Computes the confidence values for each face
+			/// - supporting point number:	stored as property 'f:num_supporting_points'
+			/// - face area:				stored as property 'f:face_area'
+			/// - covered area:				stored as property 'f:covered_area'
+			void compute(const Point_set& point_set, Polygon_mesh& mesh);
+
+		private:
+			// Returns the indices of the supporting point for 'face'
+			std::vector<std::size_t> supporting_points(Face_descriptor face, const Polygon_mesh& mesh, const Point_set& point_set);
+
+			inline FT triangle_area(const Point& p1, const Point& p2, const Point& p3) const {
+				const Vector& orth = CGAL::cross_product(Vector(p1, p2), Vector(p1, p3));
+				return std::sqrt(orth.squared_length()) * FT(0.5);
+			}
+
+			FT face_area(Face_descriptor f, const Polygon_mesh& mesh) const {
+				FT result(0);
+
+				const typename Polygon_mesh::template Property_map<Vertex_descriptor, Point>& coords = mesh.points();
+
+				Halfedge_around_face_circulator<Polygon_mesh> cir(mesh.halfedge(f), mesh), done(cir);
+				Halfedge_descriptor p_hd = *cir;
+				Vertex_descriptor p_vd = mesh.target(p_hd);
+				const Point& p = coords[p_vd];
+				++cir;
+
+				do {
+					Halfedge_descriptor q_hd = *cir;
+					Vertex_descriptor q_vd = mesh.target(q_hd);
+					const Point& q = coords[q_vd];
+
+					Halfedge_descriptor r_hd = mesh.next(q_hd);
+					Vertex_descriptor r_vd = mesh.target(r_hd);
+					const Point& r = coords[r_vd];
+
+					result += triangle_area(p, q, r);
+
+					++cir;
+				} while (cir != done);
+
+				return result;
+			}
+		};
+
+
+		//////////////////////////////////////////////////////////////////////////
+
+		// implementation
+
+		template <typename Kernel>
+		std::vector<std::size_t> Candidate_confidences<Kernel>::supporting_points(Face_descriptor face, const Polygon_mesh& mesh, const Point_set& point_set) {
+			std::vector<std::size_t> indices;
+
+			if (face == Polygon_mesh::null_face())
+				return indices;
+
+			// The supporting planar segment of each face
+			typename Polygon_mesh::template Property_map<Face_descriptor, Planar_segment*> face_supporting_segments =
+				mesh.template property_map<Face_descriptor, Planar_segment*>("f:supp_segment").first;
+
+			Planar_segment* segment = face_supporting_segments[face];
+			if (segment == nullptr)
+				return indices;
+
+			// The supporting plane of each face
+			typename Polygon_mesh::template Property_map<Face_descriptor, const Plane*> face_supporting_planes =
+				mesh.template property_map<Face_descriptor, const Plane*>("f:supp_plane").first;
+			// We do everything by projecting the point onto the face's supporting plane
+			const Plane* supporting_plane = face_supporting_planes[face];
+			CGAL_assertion(supporting_plane == segment->supporting_plane());
+
+			Polygon plg; // The projection of the face onto it supporting plane
+			const typename Polygon_mesh::template Property_map<Vertex_descriptor, Point>& coords = mesh.points();
+			Halfedge_around_face_circulator<Polygon_mesh> cir(mesh.halfedge(face), mesh), done(cir);
+			do {
+				Halfedge_descriptor hd = *cir;
+				Vertex_descriptor vd = mesh.target(hd);
+				const Point& p = coords[vd];
+				const Point2& q = supporting_plane->to_2d(p);
+
+				// Removes duplicated vertices
+				// The last point in the polygon
+				if (!plg.is_empty()) {
+					const Point2& r = plg[plg.size() - 1];
+					if (CGAL::squared_distance(q, r) < CGAL::snap_squared_distance_threshold<FT>())
+						continue;
+				}
+				plg.push_back(q);
+
+				++cir;
+			} while (cir != done);
+
+			if (plg.size() < 3 || !plg.is_simple())
+				return indices;
+
+			const typename Point_set::Point_map& points = point_set.point_map();
+			for (std::size_t i = 0; i < segment->size(); ++i) {
+				std::size_t idx = segment->at(i);
+				const Point& p = points[idx];
+				if (plg.bounded_side(supporting_plane->to_2d(p)) == CGAL::ON_BOUNDED_SIDE)
+					indices.push_back(idx);
+			}
+
+			return indices;
+		}
+
+		template <typename Kernel>
+		void Candidate_confidences<Kernel>::compute(const Point_set& point_set, Polygon_mesh& mesh) {
+			const unsigned int K = 6;
+
+			const typename Point_set::Point_map& points = point_set.point_map();
+			FT avg_spacing = compute_average_spacing<Concurrency_tag>(points, K);
+
+			// The number of supporting points of each face
+			typename Polygon_mesh::template Property_map<Face_descriptor, std::size_t> face_num_supporting_points =
+				mesh.template add_property_map<Face_descriptor, std::size_t>("f:num_supporting_points").first;
+
+			// The area of each face
+			typename Polygon_mesh::template Property_map<Face_descriptor, FT> face_areas =
+				mesh.template add_property_map<Face_descriptor, FT>("f:face_area").first;
+
+			// The point covered area of each face
+			typename Polygon_mesh::template Property_map<Face_descriptor, FT> face_covered_areas =
+				mesh.template add_property_map<Face_descriptor, FT>("f:covered_area").first;
+
+			// The supporting plane of each face
+			typename Polygon_mesh::template Property_map<Face_descriptor, const Plane*> face_supporting_planes =
+				mesh.template property_map<Face_descriptor, const Plane*>("f:supp_plane").first;
+
+			FT degenerate_face_area_threshold = CGAL::snap_squared_distance_threshold<FT>() * CGAL::snap_squared_distance_threshold<FT>();
+
+			BOOST_FOREACH(Face_descriptor f, mesh.faces()) {
+				const Plane* supporting_plane = face_supporting_planes[f];
+				// Face area
+				FT area = face_area(f, mesh);
+				face_areas[f] = area;
+
+				if (area > degenerate_face_area_threshold) {
+					const std::vector<std::size_t>& indices = supporting_points(f, mesh, point_set);
+					face_num_supporting_points[f] = indices.size();
+
+					std::vector<Point> pts;
+					for (std::size_t i = 0; i < indices.size(); ++i) {
+						std::size_t idx = indices[i];
+						const Point& p = points[idx];
+						pts.push_back(p);
+					}
+
+					FT covered_area(0);
+					Alpha_shape_mesh<Kernel> alpha_mesh(pts.begin(), pts.end(), *supporting_plane);
+					Polygon_mesh covering_mesh;
+					FT radius = avg_spacing * FT(5.0);
+					if (alpha_mesh.extract_mesh(radius * radius, covering_mesh)) {
+						// We cannot use the area of the 3D faces, because the alpha shape mesh is
+						// not perfectly planar
+						const typename Polygon_mesh::template Property_map<Vertex_descriptor, Point>& coords = covering_mesh.points();
+						BOOST_FOREACH(Face_descriptor face, covering_mesh.faces()) {
+							// We have to use the projected version
+							Polygon plg; // the projection of the face onto it supporting plane
+							Halfedge_around_face_circulator<Polygon_mesh> cir(covering_mesh.halfedge(face), covering_mesh), done(cir);
+							do {
+								Halfedge_descriptor hd = *cir;
+								Vertex_descriptor vd = covering_mesh.target(hd);
+								const Point& p = coords[vd];
+								const Point2& q = supporting_plane->to_2d(p);
+								plg.push_back(q);
+								++cir;
+							} while (cir != done);
+							covered_area += std::abs(plg.area());
+						}
+					}
+
+					face_covered_areas[f] = covered_area;
+					if (covered_area > area)
+						face_covered_areas[f] = area;
+				}
+				else { // For tiny faces, we can simple assign zero supporting points
+					face_num_supporting_points[f] = 0;
+					face_covered_areas[f] = FT(0.0);
+				}
+			}
+		}
+
+	} //namespace internal
+
+} //namespace CGAL
+
+
+#endif // CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_CANDIDATE_CONFIDENCES_H

Polygonal_surface_reconstruction/include/CGAL/internal/parameters.h

@@ -0,0 +1,42 @@
+// Copyright (c) 2018  Liangliang Nan
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// 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.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0+
+//
+// Author(s) : Liangliang Nan
+
+#ifndef CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_PARAMETERS_H
+#define CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_PARAMETERS_H
+
+#include <CGAL/license/Polygonal_surface_reconstruction.h>
+
+namespace CGAL {
+
+	/// When an intersecting point (at an edge, computed from a plane and an edge)
+	/// is very close to an existing vertex (i.e., an end point of an edge), we
+	/// snap the intersecting point to the existing vertex. This way we can avoid
+	/// many thin faces.
+	/// \note Value really doesn't matter as long as it is small (default is 1e-10).
+	///       So this parameter is not intended to be changed by the user.
+	template <class FT>
+	FT snap_squared_distance_threshold() {
+		return FT(1e-10);
+	}
+
+
+} //namespace CGAL
+
+#endif // CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_PARAMETERS_H

Polygonal_surface_reconstruction/include/CGAL/internal/point_set_with_planes.h

@@ -0,0 +1,192 @@
+// Copyright (c) 2018  Liangliang Nan
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// 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.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0+
+//
+// Author(s) : Liangliang Nan
+
+#ifndef CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_POINT_SET_WITH_PLANES_H
+#define CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_POINT_SET_WITH_PLANES_H
+
+#include <CGAL/license/Polygonal_surface_reconstruction.h>
+
+#include <CGAL/Point_set_3.h>
+#include <CGAL/property_map.h>
+#include <CGAL/linear_least_squares_fitting_3.h>
+
+#include <vector>
+
+
+/*!
+\file Point_set_with_planes.h
+*/
+
+namespace CGAL {
+
+
+	namespace internal {
+
+		// forward declaration
+		template <typename Kernel>
+		class Point_set_with_planes;
+
+
+		/**
+		*	A group of points (represented by their indices) belonging to a planar segment in a point set.
+		*/
+		template <typename Kernel>
+		class Planar_segment : public std::vector<std::size_t>
+		{
+		public:
+			typedef typename Kernel::Point_3            Point;
+			typedef typename Kernel::Plane_3			Plane;
+			typedef Point_set_with_planes<Kernel>		Point_set;
+
+		public:
+
+			// \param point_set the point set that owns this planar segment.
+			Planar_segment(Point_set* point_set) : point_set_(point_set), supporting_plane_(nullptr) {}
+			~Planar_segment() {}
+
+			Point_set* point_set() { return point_set_; }
+			void set_point_set(Point_set* point_set) { point_set_ = point_set; }
+
+			// Fits and returns the supporting plane of this planar segment
+			Plane* fit_supporting_plane() {
+				const typename Point_set::Point_map& points = point_set_->point_map();
+				std::list<Point> pts;
+				for (std::size_t i = 0; i < size(); ++i) {
+					std::size_t idx = at(i);
+					pts.push_back(points[idx]);
+				}
+
+				if (supporting_plane_)
+					delete supporting_plane_;
+				supporting_plane_ = new Plane;
+				CGAL::linear_least_squares_fitting_3(pts.begin(), pts.end(), *supporting_plane_, CGAL::Dimension_tag<0>());
+				return supporting_plane_;
+			}
+
+			// Returns the supporting plane of this planar segment.
+			// \note Returned plane is valid only if fit_supporting_plane() has been called.
+			Plane* supporting_plane() const { return supporting_plane_; }
+
+		private:
+			Point_set * point_set_;
+			Plane *		supporting_plane_; // The hypothesis generator owns this plane and manages the memory
+		};
+
+
+		/**
+		*	An enriched point set that stores the extracted planar segments
+		*/
+		template <typename Kernel>
+		class Point_set_with_planes : public Point_set_3<typename Kernel::Point_3>
+		{
+		public:
+
+			typedef Point_set_3<typename Kernel::Point_3>	Base_class;
+			typedef Point_set_with_planes<Kernel>			This_class;
+
+			typedef typename Kernel::FT						FT;
+			typedef typename Kernel::Point_3				Point;
+			typedef typename Kernel::Vector_3				Vector;
+			typedef typename Kernel::Plane_3				Plane;
+			typedef Planar_segment<Kernel>					Planar_segment;
+
+		public:
+
+			/*!
+			\tparam PointRange The range of input points.
+
+			\tparam PointMap is a model of `ReadablePropertyMap` with value
+			type `Point_3<Kernel > `.
+
+			\tparam NormalMap is a model of `ReadablePropertyMap` with value
+			type `Vector_3<Kernel > `.
+
+			\tparam PlaneIndexMap is a model of `ReadablePropertyMap` with value
+			type `int `.
+
+			\param point_range the input point range.
+
+			\param point_map property map: value_type of `InputIterator` -> Point_3.
+
+			\param normal_map property map: value_type of `InputIterator` -> Vector_3.
+
+			\param plane_index_map property map: value_type of `InputIterator` -> int,
+			denoting the index of the plane it belongs to.
+			*/
+			template <
+				typename PointRange,
+				typename PointMap,
+				typename NormalMap,
+				typename PlaneIndexMap
+			>
+				Point_set_with_planes(
+					const PointRange& points,
+					PointMap point_map,
+					NormalMap normal_map,
+					PlaneIndexMap plane_index_map
+				)
+			{
+				Base_class::resize(points.size());
+				Base_class::add_normal_map();
+
+				// Gets to know the number of plane from the plane indices
+				int max_plane_index = 0;
+				for (typename PointRange::const_iterator it = points.begin(); it != points.end(); ++it) {
+					int plane_index = get(plane_index_map, *it);
+					if (plane_index > max_plane_index)
+						max_plane_index = plane_index;
+				}
+				std::size_t num_plane = max_plane_index + 1; // the first one has index 0
+
+				for (std::size_t i = 0; i < num_plane; ++i)
+					planar_segments_.push_back(new Planar_segment(this));
+
+				std::size_t idx = 0;
+				for (typename PointRange::const_iterator it = points.begin(); it != points.end(); ++it) {
+					Base_class::m_points[idx] = get(point_map, *it);
+					Base_class::m_normals[idx] = get(normal_map, *it);
+					int plane_index = get(plane_index_map, *it);
+					if (plane_index != -1) {
+						Planar_segment* ps = planar_segments_[plane_index];
+						ps->push_back(idx);
+					}
+
+					++idx;
+				}
+			}
+
+			~Point_set_with_planes() {
+				for (std::size_t i = 0; i < planar_segments_.size(); ++i)
+					delete planar_segments_[i];
+			}
+
+			std::vector< Planar_segment* >& planar_segments() { return planar_segments_; }
+			const std::vector< Planar_segment* >& planar_segments() const { return planar_segments_; }
+
+		private:
+			std::vector< Planar_segment* > planar_segments_;
+		};
+
+	} //namespace internal
+
+} //namespace CGAL
+
+
+#endif // CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_POINT_SET_WITH_PLANES_H

Polygonal_surface_reconstruction/package_info/Polygonal_surface_reconstruction/copyright.txt

@@ -0,0 +1 @@
+Liangliang Nan

Polygonal_surface_reconstruction/package_info/Polygonal_surface_reconstruction/description.txt

@@ -0,0 +1,5 @@
+PolyFit implements the hypothesis and selection based surface reconstruction method described in the following paper:
+
+Liangliang Nan and Peter Wonka. 
+PolyFit: Polygonal Surface Reconstruction from Point Clouds. 
+ICCV 2017

Polygonal_surface_reconstruction/package_info/Polygonal_surface_reconstruction/license.txt

@@ -0,0 +1 @@
+GPL (v3 or later)

Polygonal_surface_reconstruction/package_info/Polygonal_surface_reconstruction/maintainer

@@ -0,0 +1 @@
+Liangliang Nan <liangliang.nan@gmail.com>

Polygonal_surface_reconstruction/test/Polygonal_surface_reconstruction/CMakeLists.txt

@@ -0,0 +1,64 @@
+# Created by the script cgal_create_CMakeLists
+# This is the CMake script for compiling a set of CGAL applications.
+
+project( Polygonal_surface_reconstruction_Tests )
+
+
+cmake_minimum_required(VERSION 2.8.11)
+
+# CGAL and its components
+find_package( CGAL QUIET COMPONENTS  )
+
+if ( NOT CGAL_FOUND )
+
+  message(STATUS "This project requires the CGAL library, and will not be compiled.")
+  return()  
+
+endif()
+
+# include helper file
+include( ${CGAL_USE_FILE} )
+
+
+# Boost and its components
+find_package( Boost REQUIRED )
+
+if ( NOT Boost_FOUND )
+  message(STATUS "This project requires the Boost library, and will not be compiled.")
+  return()  
+endif()
+
+find_package( SCIP )
+if ( NOT SCIP_FOUND )
+  message(STATUS "This project requires the SCIP library, and will not be compiled.")
+  return()  
+endif()
+
+
+# include for local directory
+
+# include for local package
+include_directories( BEFORE ../../include )
+include_directories( BEFORE ${SCIP_INCLUDE_DIRS} )
+
+
+# Creating entries for all C++ files with "main" routine
+# ##########################################################
+
+include( CGAL_CreateSingleSourceCGALProgram )
+
+create_single_source_cgal_program( "polygonal_surface_reconstruction_test.cpp" )
+
+target_link_libraries( polygonal_surface_reconstruction_test PRIVATE ${SCIP_LIBRARIES} )
+
+
+find_package( GLPK QUIET)
+if ( NOT GLPK_FOUND )
+  target_compile_definitions(polygonal_surface_reconstruction_test PRIVATE SKIP_TEST_USING_GLPK)
+endif()
+
+
+if ( GLPK_FOUND )
+	include_directories( BEFORE ${GLPK_INCLUDE_DIRS} )
+	target_link_libraries( polygonal_surface_reconstruction_test PRIVATE ${GLPK_LIBRARIES} )
+endif()

Polygonal_surface_reconstruction/test/Polygonal_surface_reconstruction/polygonal_surface_reconstruction_test.cpp

@@ -0,0 +1,106 @@
+// polygonal_surface_reconstruction_test.cpp
+//
+// Test the Polygonal Surface Reconstruction method with different 
+//    - kernels(Simple_cartesian, EPICK)
+//    - solvers(GLPK, SCIP)
+//    - use/ignore provided planar segmentation
+//    - input file formats (pwn, ply). For ply format, a property "segment_index"
+//      must be present storing the plane index for each point(-1 if the point is
+//      not assigned to a plane).
+
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#ifndef SKIP_TEST_USING_GLPK
+#include <CGAL/GLPK_mixed_integer_program_traits.h>
+#endif
+#include <CGAL/SCIP_mixed_integer_program_traits.h>
+
+#include "polygonal_surface_reconstruction_test_framework.h"
+
+
+// kernels:
+typedef CGAL::Simple_cartesian<double>								Cartesian;
+typedef CGAL::Exact_predicates_inexact_constructions_kernel			Epick;
+
+// solvers:
+#ifndef SKIP_TEST_USING_GLPK
+typedef CGAL::GLPK_mixed_integer_program_traits<double>				GLPK_Solver;
+#endif
+typedef CGAL::SCIP_mixed_integer_program_traits<double>				SCIP_Solver;
+
+
+int main(int argc, char * argv[])
+{
+	std::cerr << "Testing the Polygonal Surface Reconstruction method...\n\n";
+
+#if 0
+    argc = 2;
+	argv[1] = "data/icosahedron.pwn";
+ 	//argv[1] = "data/house.ply";
+	//argv[1] = "data/chair.ply";
+	//argv[1] = "data/square.ply";
+#endif
+
+	// usage
+	if (argc - 1 == 0) {
+		std::cerr << "For the input point cloud, reconstruct a water-tight polygonal surface.\n";
+		std::cerr << "\n";
+		std::cerr << "Usage: " << argv[0] << " point_cloud_file" << std::endl;
+		std::cerr << "Input file formats are \'pwn\' and \'ply\'. No output.\n\n";
+		return EXIT_FAILURE;
+	}
+
+	char* input_file = argv[1];
+
+	//---------------------------------------------------------------------
+
+	std::cerr << "--- Using Simple cartesian kernel\n\n";
+
+	//---------------------------------------------------------------------
+
+#ifndef SKIP_TEST_USING_GLPK
+	std::cerr << "\n\t---- Using GLPK solver\n\n";
+
+	std::cerr << "\t\t---- using provided planes\n";
+	reconstruct<Cartesian, GLPK_Solver>(input_file, false);
+
+	std::cerr << "\n\t\t---- re-extract planes\n";
+	reconstruct<Cartesian, GLPK_Solver>(input_file, true);
+#endif
+
+	std::cerr << "\n\t---- Using SCIP solver\n\n";
+
+	std::cerr << "\t\t---- using provided planes\n";
+	reconstruct<Cartesian, SCIP_Solver>(input_file, false);
+
+	std::cerr << "\n\t\t---- re-extract planes\n\n";
+	reconstruct<Cartesian, SCIP_Solver>(input_file, true);
+
+
+	//---------------------------------------------------------------------
+
+	std::cerr << "\n--- Using Epick kernel\n\n";
+
+	//---------------------------------------------------------------------
+
+#ifndef SKIP_TEST_USING_GLPK
+	std::cerr << "\t---- Using GLPK solver\n\n";
+
+	std::cerr << "\t\t---- using provided planes\n";
+	reconstruct<Epick, GLPK_Solver>(input_file, false);
+
+	std::cerr << "\n\t\t---- re-extract planes\n";
+	reconstruct<Epick, GLPK_Solver>(input_file, true);
+#endif 
+
+	std::cerr << "\n\t---- Using SCIP solver\n\n";
+
+	std::cerr << "\t\t---- using provided planes\n";
+	reconstruct<Epick, SCIP_Solver>(input_file, false);
+
+	std::cerr << "\n\t\t---- re-extract planes\n";
+	reconstruct<Epick, SCIP_Solver>(input_file, true);
+}

Polygonal_surface_reconstruction/test/Polygonal_surface_reconstruction/polygonal_surface_reconstruction_test_framework.h

@@ -0,0 +1,167 @@
+#ifndef POLYFIT_TEST_FRAMEWORK_H_
+#define POLYFIT_TEST_FRAMEWORK_H_
+
+#include <CGAL/IO/read_xyz_points.h>
+#include <CGAL/IO/read_ply_points.h>
+#include <CGAL/property_map.h>
+#include <CGAL/Shape_detection_3.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polygonal_surface_reconstruction.h>
+#include <CGAL/Timer.h>
+
+#include <string>
+
+
+// This function enables to run the Polygonal Surface Reconstruction algorithm with different 
+//    - kernels(Simple_cartesian, EPICK)
+//    - solvers(GLPK, SCIP)
+//    - use/ignore provided planar segmentation
+//    - input file formats (pwn, ply). For ply format, a property named "segment_index" 
+//		stores the plane index for each point(-1 if the point is not assigned to a plane).
+
+template <typename Kernel, typename MIP_Solver>
+int reconstruct(const std::string& input_file, bool force_extract_planes)
+{
+        typedef typename Kernel::Point_3				Point;
+        typedef typename Kernel::Vector_3				Vector;
+        typedef	CGAL::Polygonal_surface_reconstruction<Kernel>		Polygonal_surface_reconstruction;
+        typedef CGAL::Surface_mesh<Point>				Surface_mesh;
+
+	// Point with normal, and plane index
+        typedef boost::tuple<Point, Vector, int>			PNI;
+        typedef std::vector<PNI>					Point_vector;
+        typedef CGAL::Nth_of_tuple_property_map<0, PNI>                 Point_map;
+        typedef CGAL::Nth_of_tuple_property_map<1, PNI>			Normal_map;
+        typedef CGAL::Nth_of_tuple_property_map<2, PNI>			Plane_index_map;
+
+	typedef CGAL::Shape_detection_3::Shape_detection_traits<Kernel, Point_vector, Point_map, Normal_map>	Traits;
+        typedef CGAL::Shape_detection_3::Efficient_RANSAC<Traits>		Efficient_ransac;
+        typedef CGAL::Shape_detection_3::Plane<Traits>				Plane;
+	typedef CGAL::Shape_detection_3::Point_to_shape_index_map<Traits>	Point_to_shape_index_map;
+
+
+	Point_vector points;
+
+	// Loads point set from a file. 
+	std::ifstream input_stream(input_file.c_str());
+        if (!input_stream) {
+            std::cerr << " Error: cannot read file " << input_file << std::endl;
+            return EXIT_FAILURE;
+        }
+	std::cout << "\t\t\tLoading point cloud: " << input_file << "...";
+
+	std::string extension = input_file.substr(input_file.find_last_of('.'));
+
+	bool input_has_planes = false;
+
+	CGAL::Timer t;
+	t.start();
+	if (extension == ".pwn") {
+                if (!CGAL::read_xyz_points(
+                            input_stream,
+                            std::back_inserter(points),
+                            CGAL::parameters::point_map(Point_map()).normal_map(Normal_map())))
+		{
+                        std::cerr << " Error: cannot read file " << input_file << std::endl;
+			return EXIT_FAILURE;
+		}
+		else
+			std::cout << " Done. " << points.size() << " points. Time: " << t.time() << " sec." << std::endl;
+	}
+	else if (extension == ".ply") {
+                if (!CGAL::read_ply_points_with_properties(
+				input_stream,
+				std::back_inserter(points),
+				CGAL::make_ply_point_reader(Point_map()),
+				CGAL::make_ply_normal_reader(Normal_map()),
+				std::make_pair(Plane_index_map(), CGAL::PLY_property<int>("segment_index"))))
+		{
+                        std::cerr << " Error: cannot read file " << input_file << std::endl;
+			return EXIT_FAILURE;
+		}
+		else
+			std::cout << " Done. " << points.size() << " points. Time: " << t.time() << " sec." << std::endl;
+
+		int max_plane_index = 0;
+		for (std::size_t i = 0; i < points.size(); ++i) {
+                        int plane_index = points[i].template get<2>();
+			if (plane_index > max_plane_index)
+				max_plane_index = plane_index;
+		}
+
+		int num_planes = max_plane_index + 1;
+		// if fewer than 4 planes provided, we consider it as not provided.
+		if (num_planes > 3) {
+			input_has_planes = true;
+			std::cout << "\t\t\tInput has " << num_planes << " planes" << std::endl;
+		}
+	}
+
+	if (!input_has_planes || force_extract_planes) {
+		if (input_has_planes && force_extract_planes)
+			std::cout << "\t\t\tExtracting planes (forced)...";
+		else if (!input_has_planes)
+			std::cout << "\t\t\tExtracting planes...";
+
+		// Shape detection
+		Efficient_ransac ransac;
+		ransac.set_input(points);
+                ransac.template add_shape_factory<Plane>();
+
+		t.reset();
+		ransac.detect();
+
+                typename Efficient_ransac::Plane_range planes = ransac.planes();
+		std::size_t num_planes = planes.size();
+
+		std::cout << " Done. " << num_planes << " planes extracted. Time: " << t.time() << " sec." << std::endl;
+
+		// store the plane index of each point as the third element of the tuple.
+		Point_to_shape_index_map shape_index_map(points, planes);
+		for (std::size_t i = 0; i < points.size(); ++i) {
+			// Use the get function from the property map that accesses the 3rd element of the tuple.
+			int plane_index = get(shape_index_map, i);
+                        points[i].template get<2>() = plane_index;
+		}
+	}
+
+	//////////////////////////////////////////////////////////////////////////
+
+	std::cout << "\t\t\tGenerating candidate faces...";
+	t.reset();
+
+	Polygonal_surface_reconstruction algo(
+		points,
+		Point_map(),
+		Normal_map(),
+		Plane_index_map()
+	);
+
+	std::cout << " Done. Time: " << t.time() << " sec." << std::endl;
+
+	//////////////////////////////////////////////////////////////////////////
+
+	Surface_mesh model;
+
+	std::cout << "\t\t\tReconstructing...";
+	t.reset();
+
+        if (!algo.template reconstruct<MIP_Solver>(model)) {
+		std::cerr << " Failed: " << algo.error_message() << std::endl;
+		return EXIT_FAILURE;
+	}
+	std::cout << " Done. Time: " << t.time() << " sec." << std::endl;
+
+	if (model.is_valid()) {
+		std::cout << "\t\t\tReconstructed model has " << model.number_of_faces() << " faces" << std::endl;
+		return EXIT_FAILURE;
+	}
+	else {
+		std::cout << "\t\t\tReconstructed model is not valid. Reconstruction maybe failed?" << std::endl;
+		return EXIT_SUCCESS;
+	}
+
+}
+
+
+#endif	// POLYFIT_TEST_FRAMEWORK_H_

Solver_interface/doc/Solver_interface/Concepts/MixedIntegerProgramTraits.h

@@ -0,0 +1,394 @@
+
+
+class MixedIntegerProgramTraits
+
+
+/*!
+\cgalConcept
+
+`MixedIntegerProgramVariable` is a concept of a variable in
+a Mixed Integer Programming (MIP) problem.
+
+\cgalHasModel `CGAL::Variable<FT>`
+
+*/
+template <typename FT>
+class MixedIntegerProgramVariable
+{
+public:
+	/// \name Types
+	/// @{
+
+	/*!
+
+	*/
+	typedef unspecified_type FT;
+
+	/*!
+	A variable can be continuous, integer, or binary
+	*/
+	enum Variable_type { CONTINUOUS, INTEGER, BINARY };
+
+	/// @}
+
+	/// \name Creation 
+	/// @{
+
+	/*!
+	Constructs a variable initialized with the pointer of the solver it belongs to,
+	the variable type, lower bound, upper bound, name, and index.
+	*/
+	MixedIntegerProgramVariable(MixedIntegerProgramTraits* solver, Variable_type type, FT lb =, FT ub, const std::string& name, int idx);
+
+	/// \name Operations 
+	/// @{
+
+	/// Returns the variable type
+	Variable_type variable_type() const;
+
+	/// Sets/Changes the variable type
+	void set_variable_type(Variable_type t);
+
+	/*!
+	Returns the name of the variable.
+	*/
+	const std::string& name() const;
+
+	/*!
+	Sets the name of the variable.
+	*/
+	void set_name(const std::string& n);
+
+	/*!
+	Returns the index of the variable.
+	*/
+	int  index() const;
+
+	/*!
+	Sets the index of the variable.
+	*/
+	void set_index(int idx);
+
+	/// Returns the solver that owns this variable
+	const MixedIntegerProgramTraits* solver() const;
+	MixedIntegerProgramTraits* solver();
+
+	/// Sets the lower bound
+	void set_lower_bound(FT lb);
+
+	/// Sets the upper bound
+	void set_upper_bound(FT ub);
+
+	/// Sets both lower and upper bounds
+	void set_bounds(FT lb, FT ub);
+
+	/// Gets the lower bound
+	FT lower_bound() const;
+
+	/// Gets the upper bound
+	FT upper_bound() const;
+
+	/// Gets both lower and upper bounds
+	void get_bounds(FT& lb, FT& ub) const;
+
+	/// Gets the infinity threshold (e.g., 1e20). 
+	/// Values greater than this value are considered as infinity.
+	static FT infinity();
+
+	/// Returns the value of the variable in the current solution.
+	/// \note (1) Valid only if the program was successfully solved.
+	///       (2) If the variable is integer and rounded == true, then the 
+	///           value will be rounded to the nearest integer.
+	FT solution_value(bool rounded = false) const;
+
+	/// Sets the solution value (should be called internally by the solver).
+	void set_solution_value(FT value);
+
+	/// @}
+
+}; /* end MixedIntegerProgramVariable */
+
+
+
+   /*!
+
+   \cgalConcept
+
+   `MixedIntegerProgramLinearConstraint` is a concept of a linear
+   constraint in a Mixed Integer Programming (MIP) problem.
+
+   \cgalHasModel `CGAL::Linear_constraint<FT>`
+   */
+template <typename FT>
+class MixedIntegerProgramLinearConstraint
+{
+public:
+	/// \name Creation 
+	/// @{
+
+	/*!
+	Constructs a linear constraint, initialized with the solver it belongs to,
+	the lower bound, upper bound, name, and index.
+	*/
+	MixedIntegerProgramLinearConstraint(MixedIntegerProgramTraits* solver, FT lb, FT ub, const std::string& name, int idx);
+
+	/// \name Operations 
+	/// @{
+
+	/*!
+	Returns the name of the constraint.
+	*/
+	const std::string& name() const;
+
+	/*!
+	Sets the name of the constraint.
+	*/
+	void set_name(const std::string& n);
+
+	/*!
+	Returns the index of the constraint.
+	*/
+	int  index() const;
+
+	/*!
+	Sets the index of the constraint.
+	*/
+	void set_index(int idx);
+
+	/// Returns the solver that owns this constraint
+	const MixedIntegerProgramTraits* solver() const;
+	MixedIntegerProgramTraits* solver();
+
+	/// Sets the lower bound
+	void set_lower_bound(FT lb);
+
+	/// Sets the upper bound
+	void set_upper_bound(FT ub);
+
+	/// Sets both lower and upper bounds
+	void set_bounds(FT lb, FT ub);
+
+	/// Gets the lower bound
+	FT lower_bound() const;
+
+	/// Gets the upper bound
+	FT upper_bound() const;
+
+	/// Gets both lower and upper bounds
+	void get_bounds(FT& lb, FT& ub) const;
+
+	/// Gets the infinity threshold (e.g., 1e20). 
+	/// Values greater than this value are considered as infinity.
+	static FT infinity();
+
+	/// Sets the coefficients of the constraint. 
+	void  set_coefficients(const std::unordered_map<const MixedIntegerProgramVariable*, FT>& coeffs);
+
+	/// Adds a coefficient to a variable of the constraint. 
+	void  add_coefficient(const MixedIntegerProgramVariable* var, FT coeff);
+
+	/// Returns the coefficients of the constraint. 
+	const std::unordered_map<const MixedIntegerProgramVariable*, FT>& coefficients() const;
+
+	/// Gets the coefficient of the variable in this constraint. 
+	FT get_coefficient(const MixedIntegerProgramVariable* var) const;
+
+	/// Sets the constant term.
+	void set_offset(FT value);
+
+	/// Gets the constant term.
+	FT offset() const;
+
+	/// Clears all variables and sets the constant term to zero.
+	/// Useful to reuse the object to define a new linear constraint.
+	void clear();
+
+	/// @}
+
+}; /* end MixedIntegerProgramLinearConstraint */
+
+
+
+   /*!
+
+   \cgalConcept
+
+   `MixedIntegerProgramLinearObjective` is a concept of the linear
+   objective function in a Mixed Integer Programming (MIP) problem.
+
+   \cgalHasModel `CGAL::Linear_objective<FT>`
+   */
+template <typename FT>
+class MixedIntegerProgramLinearObjective
+{
+public:
+	/// \name Types
+	/// @{
+
+	/// The objective sense (i.e., optimization direction)
+	enum Sense { MINIMIZE, MAXIMIZE, UNDEFINED };
+
+	/// @}
+
+	/// \name Creation 
+	/// @{
+
+	/*!
+	Constructs a linear objective, initialized with the solver it belongs to
+	and the objective sense.
+	*/
+	MixedIntegerProgramLinearObjective(MixedIntegerProgramTraits* solver, Sense sense);
+
+	/// \name Operations 
+	/// @{
+
+	/// Sets the objective sense.
+	void  set_sense(Sense sense);
+
+	/// Gets the objective sense.
+	Sense sense() const;
+
+	/// Sets the coefficients of the constraint. 
+	void  set_coefficients(const std::unordered_map<const MixedIntegerProgramVariable*, FT>& coeffs);
+
+	/// Adds a coefficient to a variable of the constraint. 
+	void  add_coefficient(const MixedIntegerProgramVariable* var, FT coeff);
+
+	/// Returns the coefficients of the constraint. 
+	const std::unordered_map<const MixedIntegerProgramVariable*, FT>& coefficients() const;
+
+	/// Gets the coefficient of the variable in this constraint. 
+	FT get_coefficient(const MixedIntegerProgramVariable* var) const;
+
+	/// Sets the constant term.
+	void set_offset(FT value);
+
+	/// Gets the constant term.
+	FT offset() const;
+
+	/// Clears the objective (i.e., removes all variables, resets the 
+	/// objective sense to UNDEFINED). Useful to reuse the object 
+	/// to define a new linear objective.
+	void clear();
+
+	/// @}
+
+}; /* end MixedIntegerProgramLinearObjective */
+
+
+
+/*!
+\ingroup PkgSolverConcepts
+\cgalConcept
+
+@brief Concept describing the set of requirements for (constrained or unconstrained)
+Mixed Integer Programming (MIP) problems. A model of this concept stores the integer
+variables, linear objective, and linear constraints (if any) and provides a method
+to solve the problem.
+
+\cgalHasModel `CGAL::Mixed_integer_program_traits<T>`
+\cgalHasModel `CGAL::GLPK_mixed_integer_program_traits<T>`
+\cgalHasModel `CGAL::SCIP_mixed_integer_program_traits<T>`
+*/
+
+template <typename FT>
+class MixedIntegerProgramTraits
+{
+public:
+
+	/// \name Creation
+	/// @{
+
+	/*!
+	Default constructor.
+	*/
+	MixedIntegerProgramTraits();
+
+	/// @}
+
+	/// \name Operations
+	/// @{
+
+	/// Creates a single variable, adds it to the solver, and returns its pointer.
+	/// \note Memory is managed by the solver and will be automatically released.
+	MixedIntegerProgramVariable* create_variable(Variable_type type, FT lb, FT ub, const std::string& name);
+
+	/// Creates a set of variables, adds them to the solver, and returns their pointers.
+	/// \note (1) Variables will be given default names, e.g., x0, x1...; 	
+	///		  (2) Memory is managed by the solver and will be automatically released when the
+	///		  solver is destroyed.
+	std::vector<MixedIntegerProgramVariable*> create_variables(std::size_t n);
+
+	/// Creates a single linear constraint, adds it to the solver, and returns the pointer.
+	/// \note Memory is managed by the solver and will be automatically released when the
+	///		  solver is destroyed.
+	MixedIntegerProgramLinearConstraint* create_constraint(FT lb, FT ub, const std::string& name);
+
+	/// Creates a set of linear constraints, adds them to the solver, and returns their pointers.	
+	/// \note (1) Constraints will be given default names, e.g., c0, c1...
+	///		  (2) Memory is managed by the solver and will be automatically released when the
+	///		  solver is destroyed.
+	std::vector<MixedIntegerProgramLinearConstraint*> create_constraints(std::size_t n);
+
+	/// Creates the objective function and returns the pointer.
+	/// \note Memory is managed by the solver and will be automatically released when the
+	///		  solver is destroyed.
+	MixedIntegerProgramLinearObjective* create_objective(Sense sense);
+
+	/// Returns the number of variables
+	std::size_t number_of_variables() const;
+
+	/// Returns the variables
+	const std::vector<MixedIntegerProgramVariable*>& variables() const;
+	std::vector<MixedIntegerProgramVariable*>& variables();
+
+	/// Returns the number of constraints
+	std::size_t number_of_constraints() const;
+
+	/// Returns the constraints
+	const std::vector<MixedIntegerProgramLinearConstraint*>& constraints() const;
+	std::vector<MixedIntegerProgramLinearConstraint*>& constraints();
+
+	/// Returns the number of continuous variables
+	std::size_t number_of_continuous_variables() const;
+
+	/// Returns the number of integer variables
+	std::size_t number_of_integer_variables() const;
+
+	/// Returns the number of binary variables
+	std::size_t number_of_binary_variables() const;
+
+	/// Returns true if all variables are continuous
+	bool is_continuous() const;
+
+	/// Returns true if this is a mixed integer program
+	bool is_mixed_integer_program() const;
+
+	/// Returns true if this is an integer program
+	bool is_integer_program() const;
+
+	/// Returns true if binary program
+	bool is_binary_program() const;
+
+	/// Returns the objective
+	const MixedIntegerProgramLinearObjective * objective() const;
+	MixedIntegerProgramLinearObjective * objective();
+
+	/// Solves the program. Returns false if failed.
+	bool solve();
+
+	/// Returns the result. 
+	/// \note (1) Result is valid only if the solver succeeded.
+	///       (2) Each entry in the result corresponds to the variable with the
+	///			 same index in the program.
+	const std::vector<FT>& solution() const;
+
+	/// Returns the error message.
+	/// \note This function should be called after call to solve().
+	const std::string& error_message() const { return error_message_; }
+
+	/// Clears all variables, constraints, and the objective.
+	void clear();
+
+	/// @}
+}; /* end MixedIntegerProgramTraits */

Solver_interface/doc/Solver_interface/PackageDescription.txt

@@ -9,8 +9,8 @@
 \cgalPkgDescriptionBegin{CGAL and Solvers,PkgSolverInterface}
 \cgalPkgPicture{solver.png}
 \cgalPkgSummaryBegin
-\cgalPkgAuthors{Simon Giraudot, Pierre Alliez, Frédéric Cazals, Gaël Guennebaud, Bruno Lévy, Marc Pouget, and Laurent Saboret}
+\cgalPkgAuthors{Simon Giraudot, Pierre Alliez, Frédéric Cazals, Gaël Guennebaud, Bruno Lévy, Marc Pouget, Laurent Saboret, and Liangliang Nan}
-\cgalPkgDesc{This package provides concepts and models for solving linear systems with dense or sparse matrices.}
+\cgalPkgDesc{This package provides concepts and models for solving linear systems with dense or sparse matrices, Mixed Integer Porgramming (MIP) problems with or without constraints.}
 \cgalPkgManuals{Chapter_CGAL_and_Solvers,PkgSolverInterfaceRef}
 \cgalPkgSummaryEnd
 \cgalPkgShortInfoBegin
@@ -29,6 +29,7 @@
 - `SparseLinearAlgebraTraits_d`
 - `SparseLinearAlgebraWithFactorTraits_d`
 - `SvdTraits`
+- `MixedIntegerProgramTraits`
 
 \cgalCRPSection{Classes}
 
@@ -39,4 +40,13 @@
 - `CGAL::Eigen_sparse_matrix`
 - `CGAL::Eigen_sparse_symmetric_matrix`
 - `CGAL::Eigen_svd`
+- `CGAL::Lapack_vector`
+- `CGAL::Lapack_matrix`
+- `CGAL::Lapack_svd`
+- `CGAL::Mixed_integer_program_traits`
+- `CGAL::GLPK_mixed_integer_program_traits`
+- `CGAL::SCIP_mixed_integer_program_traits`
+- `CGAL::Variable`
+- `CGAL::Linear_constraint`
+- `CGAL::Linear_objective`
 */

Solver_interface/doc/Solver_interface/Solver_interface.txt

@@ -1,25 +1,31 @@
-namespace CGAL {
+namespace CGAL {
 /*!
 
 \mainpage User Manual 
 \anchor Chapter_CGAL_and_Solvers
 \anchor chapterSolvers
 \cgalAutoToc
-\authors Simon Giraudot, Pierre Alliez, Frédéric Cazals, Gaël Guennebaud, Bruno Lévy, Marc Pouget, and Laurent Saboret
+\authors Simon Giraudot, Pierre Alliez, Frédéric Cazals, Gaël Guennebaud, Bruno Lévy, Marc Pouget, Laurent Saboret, and Liangliang Nan
 
 Several \cgal packages have to solve linear systems with dense or
-sparse matrices. This package provides concepts and models for that
+sparse matrices, linear integer programs. This package provides 
-purpose.
+concepts and models for that purpose.
-
-We generally provide models using the \ref thirdpartyEigen
-library. Wrappers for the Eigen classes `Eigen_matrix` and
-`Eigen_vector` are also provided when needed.
 
+For linear systems, we generally provide models using the 
+\ref thirdpartyEigen library. Wrappers for the Eigen classes 
+`Eigen_matrix` and `Eigen_vector` are also provided when needed.
 It is straightforward to develop equivalent models for
 other solvers, for example those found in the <a
 href="https://software.intel.com/en-us/mkl">Intel Math Kernel
 Library (MKL)</a>.
 
+For mixed integer programs (either constrained or unconstrained),
+we provide models using \ref thirdpartySCIP and \ref thirdpartyGLPK
+libraries. It is also possible to derive new models from other
+high performance libraries, e.g., 
+<a href = "https://projects.coin-or.org/Cbc"> CBC </a>, 
+<a href = "http://www.gurobi.com/"> Gurobi </a>.
+
 
 \section SectionSolverDiagonalize Matrix Diagonalization
 
@@ -87,6 +93,26 @@ the solver:
 
 \cgalExample{Solver_interface/sparse_solvers.cpp}
 
+
+
+\section SectionMIPSolver Mixed Integer Program Solvers
+
+The concept `MixedIntegerProgramTraits` defines an interface for 
+formulating and solving (constrained or unconstrained) mixed integer 
+programs. It can also be used for general linear programs. 
+
+The field type is `double`. We provide two models of this concept: 
+`GLPK_mixed_integer_program_traits` using \ref thirdpartyGLPK and 
+`SCIP_mixed_integer_program_traits` using \ref thirdpartySCIP.
+
+Here is an example that shows how to formulate and solve a simple 
+mixed integer programs using this solver.
+
+\cgalExample{Solver_interface/mixed_integer_program.cpp}
+
+
+
+
 \section SolversHistory Implementation History
 This package is the result of the increasing needs for linear solvers
 in \cgal.  The first packages that introduced the solver concepts were
@@ -96,7 +122,10 @@ on \sc{Taucs}, \sc{LAPACK}, \sc{BLAS} and \sc{OpenNL}. Gaël Guennebaud
 then introduced new models using the \ref thirdpartyEigen library that
 became the only supported models by \cgal.  Later on the packages \ref
 PkgSurfaceMeshSkeletonization and \ref PkgSurfaceMeshDeformation
-extended the existing concepts.
+extended the existing concepts. Liangliang Nan introduced the concept
+`MixedIntegerProgramTraits` and two models for solving mixed integer
+programs when implementing the \ref PkgPolygonalSurfaceReconstruction
+package.
 
 Simon Giraudot was responsible for gathering all concepts and classes, and also wrote this user manual
 with the help of Andreas Fabri.

Solver_interface/doc/Solver_interface/dependencies

@@ -8,3 +8,4 @@ Surface_mesh_skeletonization
 Surface_mesh_deformation
 Jet_fitting_3
 Poisson_surface_reconstruction_3
+Polygonal_surface_reconstruction

Solver_interface/doc/Solver_interface/examples.txt

@@ -2,4 +2,5 @@
 \example Solver_interface/diagonalize_matrix.cpp	
 \example Solver_interface/singular_value_decomposition.cpp
 \example Solver_interface/sparse_solvers.cpp
+\example Solver_interface/mixed_integer_program.cpp
 */

Solver_interface/examples/Solver_interface/CMakeLists.txt

@@ -21,6 +21,8 @@ if ( CGAL_FOUND )
 
   create_single_source_cgal_program( "diagonalize_matrix.cpp" )
 
+  create_single_source_cgal_program( "mixed_integer_program.cpp" )
+
 else()
   
     message(STATUS "This program requires the CGAL library, and will not be compiled.")

Solver_interface/examples/Solver_interface/mixed_integer_program.cpp

@@ -0,0 +1,98 @@
+
+/* 
+* This example shows how to formulate and solve the following MIP problem:
+*
+*	Maximize
+*   Objective: x1 + 2 x2 + 3 x3 + x4
+*   Subject to
+*   	c1: - x1 +   x2 + x3 +  10 x4 <= 20
+*   	c2:   x1 - 3 x2 + x3          <= 30
+*   	c3:          x2      - 3.5 x4  =  0
+*   Bounds
+*   	0 <= x1 <= 40
+*   	2 <= x4 <= 3
+*   General
+*   	x4 is integer
+*
+*   	Expected results: x1=40; x2=10.5; x3=19.5; x4=3;
+*/
+
+
+#if 1	// Uses SCIP
+
+#include <CGAL/SCIP_mixed_integer_program_traits.h>
+typedef CGAL::SCIP_mixed_integer_program_traits<double>			MIP_Solver;
+
+#else	// Uses GLPK
+
+#include <CGAL/GLPK_mixed_integer_program_traits.h>
+typedef CGAL::GLPK_mixed_integer_program_traits<double>			MIP_Solver;
+
+#endif
+
+typedef typename MIP_Solver::Variable			Variable;
+typedef typename MIP_Solver::Linear_objective	Linear_objective;
+typedef typename MIP_Solver::Linear_constraint	Linear_constraint;
+
+#include <iostream>
+
+
+int main()
+{
+	MIP_Solver solver;
+
+	// Variable x1
+	Variable* x1 = solver.create_variable(Variable::CONTINUOUS, 0, 40, "x1");
+
+	// Variable x2
+	// You can create first (using default parameters) and then assign values.
+	Variable* x2 = solver.create_variable();
+	x2->set_name("x2");	// This is optional (a default will be given)
+
+	// Variable x3
+	Variable* x3 = solver.create_variable();	// Uses all default parameters
+	x3->set_name("x3");
+
+	// Variable x4
+	Variable* x4 = solver.create_variable(Variable::INTEGER, 2, 3, "x4");
+
+	// Objective. 
+	// Be careful this is "MAXIMIZE"
+	Linear_objective * obj = solver.create_objective(Linear_objective::MAXIMIZE);
+	obj->add_coefficient(x1, 1.0);
+	obj->add_coefficient(x2, 2.0);
+	obj->add_coefficient(x3, 3.0);
+	obj->add_coefficient(x4, 1.0);
+
+	// Constraint c1: -x1 + x2 + x3 + 10 x4 <= 20
+	Linear_constraint* c1 = solver.create_constraint(-Linear_constraint::infinity(), 20, "c1");
+	c1->add_coefficient(x1, -1);
+	c1->add_coefficient(x2, 1);
+	c1->add_coefficient(x3, 1);
+	c1->add_coefficient(x4, 10);
+
+	// Constraint c2: x1 - 3 x2 + x3 <= 30
+	Linear_constraint* c2 = solver.create_constraint(-Linear_constraint::infinity(), 30, "c2");
+	c2->add_coefficient(x1, 1);
+	c2->add_coefficient(x2, -3);
+	c2->add_coefficient(x3, 1);
+
+	// Constraint c3:  x2 - 3.5 x4 = 0
+	Linear_constraint* c3 = solver.create_constraint(0, 0, "c3");
+	c3->add_coefficient(x2, 1);
+	c3->add_coefficient(x4, -3.5);
+
+	// Solve
+	if (solver.solve()) {
+		std::cout << "result: " << std::endl;
+		const std::vector<double>& results = solver.solution();
+		for (std::size_t i = 0; i < results.size(); ++i) {
+			std::cout << "\tx" << i + 1 << ": " << results[i] << std::endl;
+		}
+		return EXIT_SUCCESS;
+	}
+	else {
+		std::cerr << "solving problem failed" << std::endl;
+		return EXIT_FAILURE;
+	}
+}

Solver_interface/include/CGAL/GLPK_mixed_integer_program_traits.h

@@ -0,0 +1,296 @@
+// Copyright (c) 2018  Liangliang Nan
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// 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.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0+
+//
+// Author(s) : Liangliang Nan
+
+#ifndef CGAL_GLPK_MIXED_INTEGER_PROGRAM_TRAITS_H
+#define CGAL_GLPK_MIXED_INTEGER_PROGRAM_TRAITS_H
+
+#include <CGAL/Mixed_integer_program_traits.h>
+
+#include "3rd_glpk/glpk.h"
+
+
+namespace CGAL {
+
+	/// \ingroup PkgSolver
+	///
+	/// This class provides an interface for formulating and solving 
+	/// mixed integer programs (constrained or unconstrained) using 
+	/// \ref thirdpartyGLPK. 
+	///
+	/// \ref thirdpartyGLPK must be available on the system.
+	/// \note For better performance, please consider using 
+	///		  `SCIP_mixed_integer_program_traits`, or derive a new 
+	///		  model from `Mixed_integer_program_traits`.
+	///
+	/// \cond SKIP_IN_MANUAL
+	/// \tparam FT Number type
+	/// \endcond
+	///
+	/// \cgalModels `MixedIntegerProgramTraits`
+	///
+	/// \sa `SCIP_mixed_integer_program_traits`
+
+	template <typename FT>
+	class GLPK_mixed_integer_program_traits : public Mixed_integer_program_traits<FT>
+	{
+		/// \cond SKIP_IN_MANUAL
+	public:
+		typedef Mixed_integer_program_traits<FT>		Base_class;
+		typedef typename Base_class::Variable			Variable;
+		typedef typename Base_class::Linear_constraint	Linear_constraint;
+		typedef typename Base_class::Linear_objective	Linear_objective;
+		typedef typename Linear_objective::Sense		Sense;
+		typedef typename Variable::Variable_type		Variable_type;
+
+	public:
+
+		/// Solves the program. Returns false if fails.
+		virtual bool solve();
+		/// \endcond
+	};
+
+	//////////////////////////////////////////////////////////////////////////
+
+	// implementation
+
+	/// \cond SKIP_IN_MANUAL
+
+	namespace internal {
+
+		// Infers "bound type" (required by GLPK) from the bounds values
+		template <typename FT>
+		int bound_type(FT lb, FT ub) {
+			typedef Variable<FT> Variable;
+
+			if (lb <= -Variable::infinity() && ub >= Variable::infinity())
+				return GLP_FR;		// free (unbounded) variable
+
+			else if (lb > -Variable::infinity() && ub >= Variable::infinity())
+				return GLP_LO;		// variable with lower bound
+
+			else if (lb <= -Variable::infinity() && ub < Variable::infinity())
+				return GLP_UP;		// variable with upper bound
+
+			else {// lb > -Variable::infinity() && ub < Variable::infinity()
+				if (lb == ub)
+					return GLP_FX;	// fixed variable
+				else
+					return GLP_DB;  // double-bounded variable
+			}
+		}
+	}
+	/// \endcond
+
+	template<typename FT>
+	bool GLPK_mixed_integer_program_traits<FT>::solve() {
+		Base_class::error_message_.clear();
+
+		glp_prob* lp = glp_create_prob();
+		if (!lp) {
+			Base_class::error_message_ = "failed creating GLPK program";
+			return false;
+		}
+
+		std::size_t num_integer_variables = 0;
+
+		// Creates variables
+
+		// This "static_cast<>" suppresses many annoying warnings: "conversion from 'size_t' to 'int', possible loss of data"
+		int num_variables = static_cast<int>(Base_class::variables_.size());
+
+		glp_add_cols(lp, num_variables);
+		for (int i = 0; i < num_variables; ++i) {
+			const Variable* var = Base_class::variables_[i];
+			glp_set_col_name(lp, i + 1, var->name().c_str());
+
+			if (var->variable_type() == Variable::INTEGER) {
+				glp_set_col_kind(lp, i + 1, GLP_IV);	// glpk uses 1-based arrays
+				++num_integer_variables;
+			}
+			else if (var->variable_type() == Variable::BINARY) {
+				glp_set_col_kind(lp, i + 1, GLP_BV);	// glpk uses 1-based arrays
+				++num_integer_variables;
+			}
+			else
+				glp_set_col_kind(lp, i + 1, GLP_CV);	// Continuous variable
+
+			FT lb, ub;
+			var->get_bounds(lb, ub);
+
+			int type = internal::bound_type(lb, ub);
+			glp_set_col_bnds(lp, i + 1, type, lb, ub);
+		}
+
+		// Adds constraints
+
+		// This "static_cast<>" suppresses many annoying warnings: "conversion from 'size_t' to 'int', possible loss of data"
+		int num_constraints = static_cast<int>(Base_class::constraints_.size());
+		glp_add_rows(lp, num_constraints);
+
+		for (int i = 0; i < num_constraints; ++i) {
+			const Linear_constraint* c = Base_class::constraints_[i];
+			const std::unordered_map<const Variable*, FT>& coeffs = c->coefficients();
+			typename std::unordered_map<const Variable*, FT>::const_iterator cur = coeffs.begin();
+
+			std::vector<int>	indices(coeffs.size() + 1, 0);		// glpk uses 1-based arrays
+			std::vector<FT> coefficients(coeffs.size() + 1, 0.0);	// glpk uses 1-based arrays
+			std::size_t idx = 1; // glpk uses 1-based arrays
+			for (; cur != coeffs.end(); ++cur) {
+				int var_idx = cur->first->index();
+				FT coeff = cur->second;
+
+				indices[idx] = var_idx + 1;	 // glpk uses 1-based arrays
+				coefficients[idx] = coeff;
+				++idx;
+			}
+
+			glp_set_mat_row(lp, i + 1, static_cast<int>(coeffs.size()), indices.data(), coefficients.data());
+
+			FT lb, ub;
+			c->get_bounds(lb, ub);
+
+			int type = internal::bound_type(lb, ub);
+			glp_set_row_bnds(lp, i + 1, type, lb, ub);
+
+			glp_set_row_name(lp, i + 1, c->name().c_str());
+		}
+
+		// Sets objective 
+
+		// Determines the coefficient of each variable in the objective function
+		const std::unordered_map<const Variable*, FT>& obj_coeffs = Base_class::objective_->coefficients();
+		typename std::unordered_map<const Variable*, FT>::const_iterator cur = obj_coeffs.begin();
+		for (; cur != obj_coeffs.end(); ++cur) {
+			int var_idx = cur->first->index();
+			FT coeff = cur->second;
+			glp_set_obj_coef(lp, var_idx + 1, coeff); // glpk uses 1-based arrays
+		}
+
+		// Sets objective function sense
+		bool minimize = (Base_class::objective_->sense() == Linear_objective::MINIMIZE);
+		glp_set_obj_dir(lp, minimize ? GLP_MIN : GLP_MAX);
+		int msg_level = GLP_MSG_ERR;
+		int status = -1;
+		if (num_integer_variables == 0) { // Continuous problem
+			glp_smcp parm;
+			glp_init_smcp(&parm);
+			parm.msg_lev = msg_level;
+			status = glp_simplex(lp, &parm);
+		}
+		else { // Solves as MIP problem
+			glp_iocp parm;
+			glp_init_iocp(&parm);
+			parm.msg_lev = msg_level;
+			parm.presolve = GLP_ON;
+			// The routine glp_intopt is a driver to the MIP solver based on the branch-and-cut method,
+			// which is a hybrid of branch-and-bound and cutting plane methods.
+			status = glp_intopt(lp, &parm);
+		}
+
+		switch (status) {
+		case 0: {
+			if (num_integer_variables == 0) { // continuous problem
+				Base_class::result_.resize(num_variables);
+				for (int i = 0; i < num_variables; ++i) {
+					FT x = glp_get_col_prim(lp, i + 1);	// glpk uses 1-based arrays
+					Variable* v = Base_class::variables_[i];
+					v->set_solution_value(x);
+					Base_class::result_[i] = x;
+				}
+			}
+			else { // MIP problem
+				Base_class::result_.resize(num_variables);
+				for (int i = 0; i < num_variables; ++i) {
+					FT x = glp_mip_col_val(lp, i + 1);		// glpk uses 1-based arrays
+
+					Variable* v = Base_class::variables_[i];
+					if (v->variable_type() != Variable::CONTINUOUS)
+						x = static_cast<int>(std::round(x));
+
+					v->set_solution_value(x);
+					Base_class::result_[i] = x;
+				}
+			}
+			break;
+		}
+
+		case GLP_EBOUND:
+			Base_class::error_message_ =
+				"Unable to start the search, because some FT-bounded variables have incorrect"
+				"bounds or some integer variables have non - integer(fractional) bounds.";
+			break;
+
+		case GLP_EROOT:
+			Base_class::error_message_ =
+				"Unable to start the search, because optimal basis for initial LP relaxation is not"
+				"provided. (This code may appear only if the presolver is disabled.)";
+			break;
+
+		case GLP_ENOPFS:
+			Base_class::error_message_ =
+				"Unable to start the search, because LP relaxation of the MIP problem instance has"
+				"no primal feasible solution. (This code may appear only if the presolver is enabled.)";
+			break;
+
+		case GLP_ENODFS:
+			Base_class::error_message_ =
+				"Unable to start the search, because LP relaxation of the MIP problem instance has"
+				"no dual feasible solution.In other word, this code means that if the LP relaxation"
+				"has at least one primal feasible solution, its optimal solution is unbounded, so if the"
+				"MIP problem has at least one integer feasible solution, its(integer) optimal solution"
+				"is also unbounded. (This code may appear only if the presolver is enabled.)";
+			break;
+
+		case GLP_EFAIL:
+			Base_class::error_message_ =
+				"The search was prematurely terminated due to the solver failure.";
+			break;
+
+		case GLP_EMIPGAP:
+			Base_class::error_message_ =
+				"The search was prematurely terminated, because the relative mip gap tolerance has been reached.";
+			break;
+
+		case GLP_ETMLIM:
+			Base_class::error_message_ =
+				"The search was prematurely terminated, because the time limit has been exceeded.";
+			break;
+
+		case GLP_ESTOP:
+			Base_class::error_message_ =
+				"The search was prematurely terminated by application. (This code may appear only"
+				"if the advanced solver interface is used.)";
+			break;
+
+		default:
+			Base_class::error_message_ =
+				"optimization was stopped with status code " + std::to_string(status);
+			break;
+		}
+
+		glp_delete_prob(lp);
+
+		return (status == 0);
+	}
+
+
+} // namespace CGAL
+
+#endif // CGAL_GLPK_MIXED_INTEGER_PROGRAM_TRAITS_H

Solver_interface/include/CGAL/Mixed_integer_program_traits.h

@@ -0,0 +1,720 @@
+// Copyright (c) 2018  Liangliang Nan
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// 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.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0+
+//
+// Author(s) : Liangliang Nan
+
+#ifndef CGAL_MIXED_INTEGER_PROGRAM_TRAITS_H
+#define CGAL_MIXED_INTEGER_PROGRAM_TRAITS_H
+
+#include <string>
+#include <vector>
+#include <unordered_map>
+#include <iostream>
+#include <iomanip>
+#include <sstream>
+
+namespace CGAL {
+
+	/// \cond SKIP_IN_MANUAL
+	// forward declaration
+	template <typename FT>
+	class Mixed_integer_program_traits;
+
+	/// The base class of solver element(e.g., Variable, Linear_constraint, and Linear_objective) in
+	/// a mixed integer program
+
+	template <typename FT>
+	class Solver_entry
+	{
+	public:
+		typedef	Mixed_integer_program_traits<FT>	Solver;
+
+	private:
+		/// A solver element (e.g., variable, constraint, objective) cannot belong to multiple solvers.
+		/// "solver" owns this entry.
+		Solver_entry(
+			Solver* solver,
+			const std::string& name = "",
+			int idx = 0
+		) : solver_(solver), name_(name), index_(idx) {}
+
+	public:
+		const std::string& name() const { return name_; }
+		void set_name(const std::string& n) { name_ = n; }
+
+		int  index() const { return index_; }
+		void set_index(int idx) { index_ = idx; }
+
+		/// The solver that owns this entry
+		const Solver* solver() const { return solver_; }
+		Solver* solver() { return solver_; }
+
+	private:
+		Solver*		solver_; // The solver that owns this entry
+		std::string	name_;
+		int			index_;
+
+		template <typename T> friend class Variable;
+		template <typename T> friend class Linear_expression;
+		template <typename T> friend class Mixed_integer_program_traits;
+	};
+
+
+	/// The base class of solver elements that have bound constraints.
+	template <typename FT>
+	class Bound
+	{
+	private:
+		Bound(FT lb = -infinity(), FT ub = +infinity());
+
+	public:
+		void set_bounds(FT lb, FT ub);
+		void set_lower_bound(FT lb) { lower_bound_ = lb; }
+		void set_upper_bound(FT ub) { upper_bound_ = ub; }
+
+		void get_bounds(FT& lb, FT& ub) const;
+		FT lower_bound() const { return lower_bound_; }
+		FT upper_bound() const { return upper_bound_; }
+
+		static FT infinity();
+
+	private:
+		FT		lower_bound_;
+		FT		upper_bound_;
+
+		static FT infinity_;
+
+		template <typename T> friend class Variable;
+		template <typename T> friend class Linear_constraint;
+		template <typename T> friend class Mixed_integer_program_traits;
+	};
+	/// \endcond
+
+	/// The variables of mixed integer programs.
+	///
+	/// \cgalModels `MixedIntegerProgramVariable`
+	template <typename FT>
+	class Variable : public Solver_entry<FT>, public Bound<FT>
+	{
+		/// \cond SKIP_IN_MANUAL
+	public:
+		enum Variable_type { CONTINUOUS, INTEGER, BINARY };
+
+		typedef	Bound<FT>							Bound;
+		typedef	Solver_entry<FT>					Solver_entry;
+		typedef	Mixed_integer_program_traits<FT>	Solver;
+
+	private:
+		/// A variable cannot belong to several solvers.
+		/// "solver" owns this variable.
+		Variable(
+			Solver* solver,
+			Variable_type type = CONTINUOUS,
+			FT lb = -Bound::infinity(),
+			FT ub = +Bound::infinity(),
+			const std::string& name = "",
+			int idx = 0
+		);
+
+	public:
+		Variable_type variable_type() const { return variable_type_; }
+		void set_variable_type(Variable_type t);
+
+		/// Returns the value of the variable in the current solution.
+		/// \note (1) Valid only if the program was successfully solved.
+		///       (2) If the variable is integer and rounded == true, then the 
+		///           value will be rounded to the nearest integer.
+		FT solution_value(bool rounded = false) const;
+
+		void set_solution_value(FT value) { solution_value_ = value; }
+
+	private:
+		Variable_type	variable_type_;
+		FT				solution_value_;
+
+		template <typename T> friend class Mixed_integer_program_traits;
+		/// \endcond
+	};
+
+
+	/// \cond SKIP_IN_MANUAL
+	/// The base class of Linear_constraint and Linear_objective.
+	template <typename FT>
+	class Linear_expression : public Solver_entry<FT>
+	{
+	public:
+		typedef Solver_entry<FT>                        Solver_entry;
+		typedef Variable<FT>				Variable;
+		typedef	Mixed_integer_program_traits<FT>	Solver;
+
+	private:
+		/// An expression cannot belong to several solvers.
+		/// "solver" owns this expression.
+		Linear_expression(
+			Solver* solver,
+			const std::string& name = "",
+			int idx = 0
+		);
+
+	public:
+		/// Adds a coefficient to a variable. 
+		void  add_coefficient(const Variable* var, FT coeff);
+
+		const std::unordered_map<const Variable*, FT>& coefficients() const { return coefficients_; }
+		void  set_coefficients(const std::unordered_map<const Variable*, FT>& coeffs) { coefficients_ = coeffs; }
+
+		FT get_coefficient(const Variable* var) const;
+
+		// The constant term
+		void set_offset(FT value) { offset_ = value; }
+		FT offset() const { return offset_; }
+
+		/// Evaluates the value of this expression at the solution found.
+		/// \note (1) valid only if the problem was successfully solved.
+		///       (2) if a variable is integer and rounded == true, then the 
+		///           variable value will be rounded to the nearest integer.
+		FT solution_value(bool rounded = false) const;
+
+		virtual void clear() { coefficients_.clear(); offset_ = 0.0; }
+
+	private:
+		std::unordered_map<const Variable*, FT>	coefficients_;
+		FT	offset_;
+
+		template <typename T> friend class Linear_constraint;
+		template <typename T> friend class Linear_objective;
+		template <typename T> friend class Mixed_integer_program_traits;
+	};
+	/// \endcond
+
+	/// The linear constraint.
+	/// \cgalModels `MixedIntegerProgramLinearConstraint`
+	template <typename FT>
+	class Linear_constraint : public Linear_expression<FT>, public Bound<FT>
+	{
+		/// \cond SKIP_IN_MANUAL
+	public:
+		typedef	Bound<FT>                               Bound;
+		typedef	Linear_expression<FT>                   Linear_expression;
+		typedef	Mixed_integer_program_traits<FT>	Solver;
+
+	private:
+		/// A constraint cannot belong to several solvers.
+		/// The "solver" owns this constraint.
+		Linear_constraint(
+			Solver* solver,
+			FT lb = -Bound::infinity(),
+			FT ub = +Bound::infinity(),
+			const std::string& name = "",
+			int idx = 0
+		);
+
+		virtual ~Linear_constraint() {}
+
+		template <typename T> friend class Mixed_integer_program_traits;
+		/// \endcond
+	};
+
+
+	/// The linear objective.
+	///
+	/// \cgalModels `MixedIntegerProgramLinearObjective`
+	///
+	template <typename FT>
+	class Linear_objective : public Linear_expression<FT>
+	{
+		/// \cond SKIP_IN_MANUAL
+	public:
+		typedef	Mixed_integer_program_traits<FT>            Solver;
+		typedef Linear_expression<FT>                       Linear_expression;
+		typedef typename Linear_expression::Solver_entry    Solver_entry;
+
+		enum Sense { MINIMIZE, MAXIMIZE, UNDEFINED };
+
+	private:
+		/// An objective cannot belong to several solvers.
+		/// "solver" owns this objective.
+		Linear_objective(Solver* solver, Sense sense);
+		virtual ~Linear_objective() {}
+
+	public:
+		void  set_sense(Sense sense) { sense_ = sense; }
+		Sense sense() const { return sense_; }
+
+		void clear();
+
+	private:
+		Sense sense_;
+
+		template <typename T> friend class Mixed_integer_program_traits;
+		/// \endcond
+	};
+
+	/// \ingroup PkgSolver
+	///
+	/// The class `Mixed_integer_program_traits` provides an interface for 
+	/// formulating and solving (constrained or unconstrained) mixed integer 
+	/// programs. It can also be used for general linear programs.
+	/// \note The solve() function is virtual and thus this class cannot be
+	///		  instantiated directly. Client code should use the inherited 
+	///       classes, i.e., `GLPK_mixed_integer_program_traits` or
+	///		  `SCIP_mixed_integer_program_traits`. Alternatively, use
+	///       `Mixed_integer_program_traits` as a base to derive a new model
+	///       (using e.g., <a href = "https://projects.coin-or.org/Cbc"> CBC </a>, 
+	///       <a href = "http://www.gurobi.com/"> Gurobi </a> for better 
+	///       performance).
+	///
+	/// \cond SKIP_IN_MANUAL
+	/// \tparam FT Number type
+	/// \endcond
+	///
+	/// \cgalModels `MixedIntegerProgramTraits`
+
+	template <typename FT>
+	class Mixed_integer_program_traits
+	{
+		/// \cond SKIP_IN_MANUAL
+	public:
+		typedef CGAL::Variable<FT>						Variable;
+		typedef CGAL::Linear_constraint<FT>				Linear_constraint;
+		typedef CGAL::Linear_objective<FT>				Linear_objective;
+		typedef typename Linear_objective::Sense		Sense;
+		typedef typename Variable::Variable_type		Variable_type;
+
+	public:
+		Mixed_integer_program_traits();
+		~Mixed_integer_program_traits();
+
+		/// Creates a single variable, add it to the solver, and returns its pointer.
+		/// \note (1) If name is empty or not provided, a default name (e.g., x0, x1...) 
+		///		  will be given.
+		///		  (2) Memory is managed by the solver and will be automatically released 
+		///		  when the solver is destroyed.
+		Variable* create_variable(
+			Variable_type type = Variable::CONTINUOUS,
+			FT lb = -Variable::infinity(),
+			FT ub = +Variable::infinity(),
+			const std::string& name = ""
+		);
+
+		/// Creates a set of variables and add them to the solver.
+		/// \note (1) Variables will be given default names, e.g., x0, x1...
+		///		  (2) Memory is managed by the solver and will be automatically released 
+		///		  when the solver is destroyed.
+		std::vector<Variable*> create_variables(std::size_t n);
+
+		/// Creates a single linear constraint, add it to the solver, and returns the pointer.
+		/// \note (1) If 'name' is empty or not provided, a default name (e.g., c0, c1...) will be given.
+		///		  (2) Memory is managed by the solver and will be automatically released when the
+		///		  solver is destroyed.
+		Linear_constraint* create_constraint(
+			FT lb = -Variable::infinity(),
+			FT ub = +Variable::infinity(),
+			const std::string& name = ""
+		);
+
+		/// Creates a set of linear constraints and add them to the solver.	
+		/// \note (1) Constraints with be given default names, e.g., c0, c1...
+		///		  (2) Memory is managed by the solver and will be automatically released 
+		///		  when the solver is destroyed.
+		std::vector<Linear_constraint*> create_constraints(std::size_t n);
+
+		/// Creates the objective function and returns the pointer.
+		/// \note Memory is managed by the solver and will be automatically released 
+		///		  when the solver is destroyed.
+		Linear_objective* create_objective(Sense sense = Linear_objective::MINIMIZE);
+
+		std::size_t number_of_variables() const { return variables_.size(); }
+		const std::vector<Variable*>& variables() const { return variables_; }
+		std::vector<Variable*>& variables() { return variables_; }
+
+		std::size_t number_of_constraints() const { return constraints_.size(); }
+		const std::vector<Linear_constraint*>& constraints() const { return constraints_; }
+		std::vector<Linear_constraint*>& constraints() { return constraints_; }
+
+		const Linear_objective* objective() const;
+		Linear_objective* objective();
+
+		std::size_t number_of_continuous_variables() const;
+		std::size_t number_of_integer_variables() const;
+		std::size_t number_of_binary_variables() const;
+
+		bool is_continuous() const;				// Returns true if all variables are continuous
+		bool is_mixed_integer_program() const;	// Returns true if mixed integer program
+		bool is_integer_program() const;		// Returns true if integer program
+		bool is_binary_program() const;			// Returns true if binary program
+
+		/// Clears all variables, constraints, and the objective.
+		void clear();
+
+		/// Solves the program. Returns false if fails.
+		virtual bool solve() = 0;
+
+		/// Returns the result. 
+		/// The result can also be retrieved using Variable::solution_value().
+		/// \note (1) Result is valid only if the solver succeeded.
+		///       (2) Each entry in the result corresponds to the variable with the
+		///			 same index in the program.
+		const std::vector<FT>& solution() const { return result_; }
+
+		/// Returns the error message.
+		/// \note This function should be called after call to solve().
+		const std::string& error_message() const { return error_message_; }
+
+	protected:
+		Linear_objective*				objective_;
+		std::vector<Variable*>			variables_;
+		std::vector<Linear_constraint*>	constraints_;
+
+		std::vector<FT>		result_;
+		std::string			error_message_;
+		/// \endcond
+	};
+
+	//////////////////////////////////////////////////////////////////////////
+
+	// implementation
+
+	template<typename FT>
+	FT Bound<FT>::infinity_ = 1e20;		// values larger than this value are considered infinity
+
+	template<typename FT>
+	Bound<FT>::Bound(FT lb /* = -infinity() */, FT ub /* = +infinity() */)
+		: lower_bound_(lb)
+		, upper_bound_(ub)
+	{
+	}
+
+	template<typename FT>
+	FT Bound<FT>::infinity() {
+		return infinity_;
+	}
+
+	template<typename FT>
+	void Bound<FT>::set_bounds(FT lb, FT ub) {
+		lower_bound_ = lb;
+		upper_bound_ = ub;
+	}
+
+	template<typename FT>
+	void Bound<FT>::get_bounds(FT& lb, FT& ub) const {
+		lb = lower_bound_;
+		ub = upper_bound_;
+	}
+
+
+	template<typename FT>
+	Variable<FT>::Variable(
+		Solver* solver,
+		Variable_type type /* = CONTINUOUS */,
+		FT lb /* = -infinity() */,
+		FT ub /* = +infinity() */,
+		const std::string& name /* = "" */,
+		int idx /* = 0*/
+	)
+		: Solver_entry(solver, name, idx)
+		, Bound(lb, ub)
+		, variable_type_(type)
+		, solution_value_(0.0)
+	{
+		if (type == BINARY)
+			Bound::set_bounds(0.0, 1.0);
+	}
+
+
+	template<typename FT>
+	void Variable<FT>::set_variable_type(Variable_type type) {
+		variable_type_ = type;
+		if (type == BINARY)
+			Bound::set_bounds(0.0, 1.0);
+	}
+
+
+	template<typename FT>
+	FT Variable<FT>::solution_value(bool rounded /* = false*/) const {
+		if (rounded && variable_type_ != CONTINUOUS)
+			return std::round(solution_value_);
+		else
+			return solution_value_;
+	}
+
+
+	//////////////////////////////////////////////////////////////////////////
+
+	template<typename FT>
+	Linear_expression<FT>::Linear_expression(Solver* solver, const std::string& name, int idx)
+		: Solver_entry(solver, name, idx)
+		, offset_(0.0)
+	{
+	}
+
+
+	template<typename FT>
+	void Linear_expression<FT>::add_coefficient(const Variable* var, FT coeff) {
+		if (coefficients_.find(var) == coefficients_.end())
+			coefficients_[var] = coeff;
+		else
+			coefficients_[var] += coeff;
+	}
+
+
+	template<typename FT>
+	FT Linear_expression<FT>::get_coefficient(const Variable* var) const {
+		typename std::unordered_map<const Variable*, FT>::const_iterator pos = coefficients_.find(var);
+		if (pos != coefficients_.end())
+			return pos->second;
+		else {
+			std::cerr << "linear expression does not own variable " << var->name() << " (" << var->index() << ")" << std::endl;
+			return 0.0;
+		}
+	}
+
+
+	template<typename FT>
+	FT Linear_expression<FT>::solution_value(bool rounded /* = false*/) const {
+		FT solution = offset_;
+
+		typename std::unordered_map<const Variable*, FT>::const_iterator it = coefficients_.begin();
+		for (; it != coefficients_.end(); ++it) {
+			const Variable* var = it->first;
+			FT coeff = it->second;
+			solution += var->solution_value(rounded) * coeff;
+		}
+		return solution;
+	}
+
+
+	template<typename FT>
+	void Linear_objective<FT>::clear() {
+		Linear_expression::clear();
+		Solver_entry::set_name("");
+		Solver_entry::set_index(0);
+	}
+
+
+	template<typename FT>
+	Linear_constraint<FT>::Linear_constraint(
+		Solver* solver,
+		FT lb /* = -infinity() */,
+		FT ub /* = +infinity() */,
+		const std::string& name/* = "" */,
+		int idx /* = 0*/
+	)
+		: Linear_expression(solver, name, idx)
+		, Bound(lb, ub)
+	{
+	}
+
+	template<typename FT>
+	Linear_objective<FT>::Linear_objective(Solver* solver, Sense sense)
+		: Linear_expression(solver)
+		, sense_(sense)
+	{
+	}
+
+	template<typename FT>
+	Mixed_integer_program_traits<FT>::Mixed_integer_program_traits() {
+		// Intentionally set the objective to UNDEFINED, so it will allow me to warn
+		// the user if he/she forgot to set the objective sense.
+		objective_ = new Linear_objective(this, Linear_objective::UNDEFINED);
+	}
+
+	template<typename FT>
+	Mixed_integer_program_traits<FT>::~Mixed_integer_program_traits() {
+		clear();
+		delete objective_;
+	}
+
+	template<typename FT>
+	void Mixed_integer_program_traits<FT>::clear() {
+		for (std::size_t i = 0; i < variables_.size(); ++i)
+			delete variables_[i];
+		variables_.clear();
+
+		for (std::size_t i = 0; i < constraints_.size(); ++i)
+			delete constraints_[i];
+		constraints_.clear();
+
+		objective_->clear();
+	}
+
+	/// \cond SKIP_IN_MANUAL
+	namespace internal {
+		/**
+		* Converts an integer v to a string of specified 'width' by
+		* filling with character 'fill'
+		*/
+		template <class Int>
+		inline std::string from_integer(Int v, int width, char fill) {
+			std::ostringstream string_stream;
+			string_stream << std::setfill(fill) << std::setw(width) << v;
+			return string_stream.str();
+		}
+	}
+	/// \endcond
+
+	template<typename FT>
+	typename Mixed_integer_program_traits<FT>::Variable* Mixed_integer_program_traits<FT>::create_variable(
+		Variable_type type /* = Variable::CONTINUOUS */,
+		FT lb /* = -Variable::infinity() */,
+		FT ub /* = +Variable::infinity() */,
+		const std::string& name /* = "" */)
+	{
+		Variable* v = new Variable(this, type, lb, ub);
+
+		std::size_t idx = variables_.size();
+		v->set_index(static_cast<int>(idx));
+
+		const std::string& fixed_name = name.empty() ? "x" + internal::from_integer(idx, 9, '0') : name;
+		v->set_name(fixed_name);
+
+		variables_.push_back(v);
+		return v;
+	}
+
+	template<typename FT>
+	std::vector<typename Mixed_integer_program_traits<FT>::Variable*> Mixed_integer_program_traits<FT>::create_variables(std::size_t n) {
+		std::vector<Variable*> variables;
+		for (std::size_t i = 0; i < n; ++i) {
+			Variable* v = create_variable();
+			variables.push_back(v);
+		}
+		return variables;
+	}
+
+	template<typename FT>
+	typename Mixed_integer_program_traits<FT>::Linear_constraint* Mixed_integer_program_traits<FT>::create_constraint(
+		FT lb /* = -Variable::infinity() */,
+		FT ub /* = +Variable::infinity() */,
+		const std::string& name /* = "" */)
+	{
+		Linear_constraint* c = new Linear_constraint(this, lb, ub);
+
+		std::size_t idx = constraints_.size();
+		c->set_index(static_cast<int>(idx));
+
+		const std::string& fixed_name = name.empty() ? "c" + internal::from_integer(idx, 9, '0') : name;
+		c->set_name(fixed_name);
+
+		constraints_.push_back(c);
+		return c;
+	}
+
+	template<typename FT>
+	std::vector<typename Mixed_integer_program_traits<FT>::Linear_constraint*> Mixed_integer_program_traits<FT>::create_constraints(std::size_t n) {
+		std::vector<Linear_constraint*> constraints;
+		for (std::size_t i = 0; i < n; ++i) {
+			Linear_constraint* v = create_constraint();
+			constraints.push_back(v);
+		}
+		return constraints;
+	}
+
+	template<typename FT>
+	typename Mixed_integer_program_traits<FT>::Linear_objective * Mixed_integer_program_traits<FT>::create_objective(Sense sense /* = Linear_objective ::MINIMIZE*/) {
+		if (objective_)
+			delete objective_;
+
+		objective_ = new Linear_objective(this, sense);
+		return objective_;
+	}
+
+
+	template<typename FT>
+	const typename Mixed_integer_program_traits<FT>::Linear_objective * Mixed_integer_program_traits<FT>::objective() const {
+		if (!objective_)
+			std::cerr << "please call \'create_objective()\' to create an objective first" << std::endl;
+
+		return objective_;
+	}
+
+	template<typename FT>
+	typename Mixed_integer_program_traits<FT>::Linear_objective * Mixed_integer_program_traits<FT>::objective() {
+		if (!objective_)
+			std::cerr << "please call \'create_objective()\' to create an objective first" << std::endl;
+
+		return objective_;
+	}
+
+	template<typename FT>
+	std::size_t Mixed_integer_program_traits<FT>::number_of_continuous_variables() const {
+		std::size_t num_continuous_var = 0;
+		for (std::size_t i = 0; i < variables_.size(); ++i) {
+			const Variable* v = variables_[i];
+			if (v->variable_type() == Variable::CONTINUOUS)
+				++num_continuous_var;
+		}
+		return num_continuous_var;
+	}
+
+	template<typename FT>
+	std::size_t Mixed_integer_program_traits<FT>::number_of_integer_variables() const {
+		std::size_t num_iteger_var = 0;
+		for (std::size_t i = 0; i < variables_.size(); ++i) {
+			const Variable* v = variables_[i];
+			if (v->variable_type() == Variable::INTEGER)
+				++num_iteger_var;
+		}
+		return num_iteger_var;
+	}
+
+	template<typename FT>
+	std::size_t Mixed_integer_program_traits<FT>::number_of_binary_variables() const {
+		std::size_t num_binary_var = 0;
+		for (std::size_t i = 0; i < variables_.size(); ++i) {
+			const Variable* v = variables_[i];
+			if (v->variable_type() == Variable::BINARY)
+				++num_binary_var;
+		}
+		return num_binary_var;
+	}
+
+	// Returns true if all variables are continuous
+	template<typename FT>
+	bool Mixed_integer_program_traits<FT>::is_continuous() const {
+		std::size_t num = number_of_continuous_variables();
+		return (num > 0) && (num == variables_.size());
+	}
+
+
+	// Returns true if this is a mixed integer program
+	template<typename FT>
+	bool Mixed_integer_program_traits<FT>::is_mixed_integer_program() const {
+		std::size_t num = number_of_continuous_variables();
+		return (num > 0) && (num < variables_.size());
+	}
+
+
+	// Returns true if inter program
+	template<typename FT>
+	bool Mixed_integer_program_traits<FT>::is_integer_program() const {
+		std::size_t num = number_of_integer_variables();
+		return (num > 0) && (num == variables_.size());
+	}
+
+
+	// Returns true if binary program
+	template<typename FT>
+	bool Mixed_integer_program_traits<FT>::is_binary_program() const {
+		std::size_t num = number_of_binary_variables();
+		return (num > 0) && (num == variables_.size());
+	}
+
+} // namespace CGAL
+
+#endif // CGAL_MIXED_INTEGER_PROGRAM_TRAITS_H

Solver_interface/include/CGAL/SCIP_mixed_integer_program_traits.h

@@ -0,0 +1,244 @@
+// Copyright (c) 2018  Liangliang Nan
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// 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.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0+
+//
+// Author(s) : Liangliang Nan
+
+#ifndef CGAL_SCIP_MIXED_INTEGER_PROGRAM_TRAITS_H
+#define CGAL_SCIP_MIXED_INTEGER_PROGRAM_TRAITS_H
+
+#include <CGAL/Mixed_integer_program_traits.h>
+
+#include "scip/scip.h"
+#include "scip/scipdefplugins.h"
+
+namespace CGAL {
+
+	/// \ingroup PkgSolver
+	///
+	/// This class provides an interface for formulating and solving 
+	/// mixed integer programs (constrained or unconstrained) using 
+	/// \ref thirdpartySCIP. 
+	///
+	/// \ref thirdpartySCIP must be available on the system.
+	///
+	/// \cond SKIP_IN_MANUAL
+	/// \tparam FT Number type
+	/// \endcond
+	///
+	/// \cgalModels `MixedIntegerProgramTraits`
+	///
+	/// \sa `GLPK_mixed_integer_program_traits`
+
+	template <typename FT>
+	class SCIP_mixed_integer_program_traits : public Mixed_integer_program_traits<FT>
+	{
+		/// \cond SKIP_IN_MANUAL
+	public:
+		typedef Mixed_integer_program_traits<FT>		Base_class;
+		typedef typename Base_class::Variable			Variable;
+		typedef typename Base_class::Linear_constraint	Linear_constraint;
+		typedef typename Base_class::Linear_objective	Linear_objective;
+		typedef typename Linear_objective::Sense		Sense;
+		typedef typename Variable::Variable_type		Variable_type;
+
+	public:
+
+		/// Solves the program. Returns false if fails.
+		virtual bool solve();
+		/// \endcond
+	};
+
+	//////////////////////////////////////////////////////////////////////////
+
+	// implementation
+
+	template<typename FT>
+	bool SCIP_mixed_integer_program_traits<FT>::solve()
+	{
+		Base_class::error_message_.clear();
+
+		Scip* scip = 0;
+		SCIP_CALL(SCIPcreate(&scip));
+		SCIP_CALL(SCIPincludeDefaultPlugins(scip));
+
+		// Disables scip output to stdout
+		SCIPmessagehdlrSetQuiet(SCIPgetMessagehdlr(scip), TRUE);
+
+		// Uses wall clock time because getting CPU user seconds
+		// involves calling times() which is very expensive
+		SCIP_CALL(SCIPsetIntParam(scip, "timing/clocktype", SCIP_CLOCKTYPE_WALL));
+
+		// Creates empty problem 
+		SCIP_CALL(SCIPcreateProbBasic(scip, "Solver_interface"));
+
+		// Creates variables
+		std::vector<SCIP_VAR*> scip_variables;
+		for (std::size_t i = 0; i < Base_class::variables_.size(); ++i) {
+			const Variable* var = Base_class::variables_[i];
+			SCIP_VAR* v = 0;
+
+			double lb, ub;
+			var->get_bounds(lb, ub);
+
+			switch (var->variable_type())
+			{
+			case Variable::CONTINUOUS:
+				SCIP_CALL(SCIPcreateVar(scip, &v, var->name().c_str(), lb, ub, 0.0, SCIP_VARTYPE_CONTINUOUS, TRUE, FALSE, 0, 0, 0, 0, 0));
+				break;
+			case Variable::INTEGER:
+				SCIP_CALL(SCIPcreateVar(scip, &v, var->name().c_str(), lb, ub, 0.0, SCIP_VARTYPE_INTEGER, TRUE, FALSE, 0, 0, 0, 0, 0));
+				break;
+			case Variable::BINARY:
+				SCIP_CALL(SCIPcreateVar(scip, &v, var->name().c_str(), 0, 1, 0.0, SCIP_VARTYPE_BINARY, TRUE, FALSE, 0, 0, 0, 0, 0));
+				break;
+			}
+			// Adds the SCIP_VAR object to the scip problem
+			SCIP_CALL(SCIPaddVar(scip, v));
+
+			// Stores the SCIP_VAR pointer for later access
+			scip_variables.push_back(v);
+		}
+
+		// Adds constraints
+
+		std::vector<SCIP_CONS*> scip_constraints;
+		for (std::size_t i = 0; i < Base_class::constraints_.size(); ++i) {
+			const Linear_constraint* c = Base_class::constraints_[i];
+			const std::unordered_map<const Variable*, double>& coeffs = c->coefficients();
+			typename std::unordered_map<const Variable*, double>::const_iterator cur = coeffs.begin();
+
+			std::vector<SCIP_VAR*>	cstr_variables(coeffs.size());
+			std::vector<double>		cstr_values(coeffs.size());
+			std::size_t idx = 0;
+			for (; cur != coeffs.end(); ++cur) {
+				std::size_t var_idx = cur->first->index();
+				double coeff = cur->second;
+				cstr_variables[idx] = scip_variables[var_idx];
+				cstr_values[idx] = coeff;
+				++idx;
+			}
+
+			// Creates SCIP_CONS object
+			SCIP_CONS* cons = 0;
+			const std::string& name = c->name();
+
+			double lb, ub;
+			c->get_bounds(lb, ub);
+
+			SCIP_CALL(SCIPcreateConsLinear(scip, &cons, name.c_str(), int(coeffs.size()), cstr_variables.data(), cstr_values.data(), lb, ub, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE));
+
+			// Adds the constraint to scip
+			SCIP_CALL(SCIPaddCons(scip, cons));			
+
+			// Stores the constraint for later on
+			scip_constraints.push_back(cons);
+		}
+
+		// Sets objective
+
+		// Determines the coefficient of each variable in the objective function
+		const std::unordered_map<const Variable*, double>& obj_coeffs = Base_class::objective_->coefficients();
+		typename std::unordered_map<const Variable*, double>::const_iterator cur = obj_coeffs.begin();
+		for (; cur != obj_coeffs.end(); ++cur) {
+			const Variable* var = cur->first;
+			double coeff = cur->second;
+			SCIP_CALL(SCIPchgVarObj(scip, scip_variables[var->index()], coeff));
+		}
+
+		// Sets the objective sense
+		bool minimize = (Base_class::objective_->sense() == Linear_objective::MINIMIZE);
+		SCIP_CALL(SCIPsetObjsense(scip, minimize ? SCIP_OBJSENSE_MINIMIZE : SCIP_OBJSENSE_MAXIMIZE));
+
+		// Turns presolve on (it's the SCIP default).
+		bool presolve = true;
+		if (presolve)
+			SCIP_CALL(SCIPsetIntParam(scip, "presolving/maxrounds", -1)); // maximal number of presolving rounds (-1: unlimited, 0: off)
+		else
+			SCIP_CALL(SCIPsetIntParam(scip, "presolving/maxrounds", 0));  // disable presolve
+
+		bool status = false;
+		// This tells scip to start the solution process
+		if (SCIPsolve(scip) == SCIP_OKAY) {
+			// Gets the best found solution from scip
+			SCIP_SOL* sol = SCIPgetBestSol(scip);
+			if (sol) {
+				// If optimal or feasible solution is found.
+				Base_class::result_.resize(Base_class::variables_.size());
+				for (std::size_t i = 0; i < Base_class::variables_.size(); ++i) {
+					FT x = SCIPgetSolVal(scip, sol, scip_variables[i]);
+
+					Variable* v = Base_class::variables_[i];
+					if (v->variable_type() != Variable::CONTINUOUS)
+						x = static_cast<int>(std::round(x));
+
+					v->set_solution_value(x);
+					Base_class::result_[i] = x;
+				}
+				status = true;
+			}
+		}
+
+		// Reports the status: optimal, infeasible, etc.
+		SCIP_STATUS scip_status = SCIPgetStatus(scip);
+		switch (scip_status) {
+		case SCIP_STATUS_OPTIMAL:
+			// Provides info only if fails.
+			break;
+		case SCIP_STATUS_GAPLIMIT:
+			// To be consistent with the other solvers.
+			// Provides info only if fails.
+			break;
+		case SCIP_STATUS_INFEASIBLE:
+			Base_class::error_message_ = "model was infeasible";
+			break;
+		case SCIP_STATUS_UNBOUNDED:
+			Base_class::error_message_ = "model was unbounded";
+			break;
+		case SCIP_STATUS_INFORUNBD:
+			Base_class::error_message_ = "model was either infeasible or unbounded";
+			break;
+		case SCIP_STATUS_TIMELIMIT:
+			Base_class::error_message_ = "aborted due to time limit";
+			break;
+		default:
+			Base_class::error_message_ = "aborted with status: " + std::to_string(scip_status);
+			break;
+		}
+
+		SCIP_CALL(SCIPresetParams(scip));
+
+		// Since the SCIPcreateVar captures all variables, we have to release them now
+		for (std::size_t i = 0; i < scip_variables.size(); ++i)
+			SCIP_CALL(SCIPreleaseVar(scip, &scip_variables[i]));
+		scip_variables.clear();
+
+		// The same for the constraints
+		for (std::size_t i = 0; i < scip_constraints.size(); ++i)
+			SCIP_CALL(SCIPreleaseCons(scip, &scip_constraints[i]));
+		scip_constraints.clear();
+
+		// After releasing all vars and cons we can free the scip problem.
+		// Remember this has always to be the last call to scip
+		SCIP_CALL(SCIPfree(&scip));
+
+		return status;
+	}
+
+} // namespace CGAL
+
+#endif // CGAL_SCIP_MIXED_INTEGER_PROGRAM_TRAITS_H

Solver_interface/package_info/Solver_interface/description.txt

@@ -1,2 +1,2 @@
-Solver_interface: contains classes to interface third party linear algebra solvers
+Solver_interface: contains classes to interface third party linear algebra solvers (such as Eigen for example),
-to CGAL (such as Eigen for example)
+mixed integer program solvers (such as SCIP for example) to CGAL.

Stream_support/include/CGAL/IO/PLY.h

@@ -24,6 +24,8 @@
 #include <CGAL/Kernel_traits.h>
 #include <CGAL/IO/io.h>
 
+#include <boost/property_map/property_map.hpp>
+
 #include <tuple>
 #include <iostream>
 #include <sstream>