Personal tools
You are here: Home Courses Computational Geometry Fall 2016/2017 Assignments Programming Assignment
« March 2021 »
Log in

Forgot your password?

Programming Assignment

Note: During the Last 24 hours I've made several changes to this webpage, mainly uploaded test cases. You can obtain these test cases and try out your program on them. I've created archive files of the input files for your convenience.

While editing this webpage I changed few things back and forth. I suggest that you revisit the material you have obtained. In particular I've slightly modified the exporter methods; see 2nd item in the Q&A section.


Your solution to this exercise should consist of at least one source file, CMakeLists.txt (input to cmake), and a pdf file called max_cover.pdf. The latter should provide some basic information about your program. It should also describe the algorithm implemented.

Running cmake on the provided CMakeLists.txt followed by compilation and linkage should produce an executable called max_cover.

If you organize all your code in a single source file, say max_cover.cpp, and place this file in a dedicated directory that does not contain anything else, you can simply run:

<CGAL>/Scripts/scripts/cgal_create_cmake_script example

to generate CMakeLists.txt, where <CGAL> is the directory where CGAL was installed from sources at. If you want to use C++11 features, edit the generated CMakeLists.txt file, and add the following lines:


Alternatively, you can use the CMakeLists.txt provided as an example below.

The program max_cover accepts two arguments on the command line, namely, the square of the radius as a rational number and the name of the file that contains the input points. By default, the file is called points.txt. A rational number is given either as a decimal number or as a numerator and denominator separated by '/' (e.g., 1/3), where the numerator and denominator are both unlimited integers. When given as a decimal number it is first converted to the nearest floating-point number, and then to a rational number represented by a numerator and a denominator. Observe that such numbers may have long bit lengths, and, as a consequence, arithmetic operations on them may be slow.

You can develop on Windows using VC10, VC11, VC12, or VC14, but the code must compile and run using g++ on Linux.

You must use a single thread.

Input, Output, and more

The input file is a text file (as its extension suggests). The format of the file that contains the disc center points follows:

n x0 y0 x1 y1 ... xn-1 yn-1

where n is the number of distinct points and each point is given by its x and y coordinates, which are rational numbers.

There will be a contest between all the different programs (including ours). To this end you are also required to add to your code a timer that measures the execution time-consumption.

Use a Boost timer.

  • Add the directive below at the top of your cpp source file:
    #include <boost/timer.hpp>
  • Introduce a timer after reading the input files:
    boost::timer timer;
  • At the end of the code that executes the algorithm (and before the start of the code that prints out the final result) measure the elapsed time:
    double secs = timer.elapsed();

An example of a program that accepts the input described above, constructs an arrangement induced by circles, of the given radius, centered at the input points can be found here.

An example of a CMakeLists.txt file that can be used to build the program above can be found here.

The output of the program should consists of three text lines.

In the following output quantities and coordinates are embedded in '<' and '>', and thus the '<' and '>' characters should not be printed out.

The first line should be:

"Execution time: <s>"

where <s> is the timer elapsed time in seconds.

The second line should list the number of points covered by the maximum covering disc:

"The maximum covering disc covers <k> points"

The third line should list the lexicographically smallest center of all discs with the same maximum number of covered points:

"The maximum covering disc is centered at <c>"

Use the exporter operator

std::ostream& operator<<(std::ostream& os, const Point& p)

defined in the example program max_cover.cpp above to export the center point to the standard output.

The coordinates of the points in the arrangement constructed in the example program max_cover.cpp are algebraic numbers that are not necessarily rational. However, they are a "light" vertion of algebraic numbers, which can be expressed as a + b * sqrt(c), where a, b, and c are rational numbers. The CGAL template CGAL::Sqrt_extension (see is used to represent such numbers. It supports the arithmetic operations '+', '-', '*', and '/'.


As part of the contest the performance of your programs will be compared on various inputs.

Make sure your program can handle at least several hundreds input points in various special positions.

The numbers under the groups (6 right columns) indicate the time consumption in seconds.

The first type of input consists of sets, such that each set contains n2 points on an n x n grid.

Square Radius
n Input Points Covered Points
Disc Center
Efi Ilia & Tzvika
Michael & Omri
Omer & Gal
Shahaf & Yair
1 1  points_grid_1.txt  1 ((0/1,-1/1,1/1),(0/1,0/1,0/1)) 0.000171
 0.000557 v 0.000411 0.000435 0.00026
1 3  points_grid_9.txt  5 ((0/1,1/1,1/1),(1/1,0/1,0/1))
 0.002767  0.003549 v 0.005224 0.001694 0.003309
1/4 3  points_grid_9.txt  2 ((0/1,0/1,0/1),(1/2,0/1,0/1)) 0.000531 0.002722  x 0.000745 0.00065
9 3  points_grid_9.txt  9 ((2/1,-1/4,128/1),(1/1,0/1,128/1)) 0.004545 0.008600 v 0.005172 0.005743 0.006004
1 4  points_grid_16.txt  5 ((0/1,1/1,1/1),(1/1,0/1,0/1)) 0.004097 0.005276  x 0.005451 0.003147
16 4  points_grid_16.txt  16 ((3/1,-1/6,495/1),(3/2,0/1,495/1)) 0.007967 0.030261  v 0.011857 0.019906
12  points_grid_144.txt  5 ((0/1,1/1,1/1),(1/1,0/1,0/1))
0.021015 0.044852  x  0.03297 0.026561
12  points_grid_144.txt  144 ((11/1,-1/22,55055/1),(11/2,0/1,55055/1)) 0.553085 3.687130  v 0.941872 2.42152 2.77307
20  points_grid_400.txt  5 ((0/1,1/1,1/1),(1/1,0/1,0/1))
0.044532 0.120384  x 0.101102 0.06623 0.124325
20  points_grid_400.txt  400 ((19/1,-1/38,447279/1),(19/2,0/1,447279/1))
5.21071 31.6068  v 8.93174 36.5862
31  points_grid_961.txt  5 ((0/1,1/1,1/1),(1/1,0/1,0/1)) 0.109369 0.318584  x 0.232092 0.155417
31  points_grid_961.txt  961  ((30/1,-1/60,2.6496e+06/1),(15/1,0/1,2.6496e+06/1)) 36.6941 251.273  v 70.7738 158.519

The archive file points_grid.tar.gz contains all the files in the table above and a short script that runs max_cover on them with the appropriate radius square.

The second type of input consists of sets, such that each set contains n non-antipodal points on a circle of radius 1.

Square Radius n
Input Points Covered Points Disc Center Efi Ilia & Tzvika Michael & Omri Mor Omer & Gal
Shahaf & Yair
1 3 points_circle_3.txt
3 ((-32/65,-2/7,12544/4225),(4/65,1/28,12544/4225))
 0.006508 v  0.002417  0.002502  0.00139
1 4 points_circle_4.txt
4 ((0/1,0/1,0/1),(0/1,0/1,0/1))
 0.001012  0.003825  x  0.003446  0.003854  0.006174
1 8 points_circle_8.txt
8 ((0/1,0/1,0/1),(0/1,0/1,0/1))  0.002746  0.017415  x  0.005636  0.01075  0.009475
1 16 points_circle_16.txt
16 ((0/1,0/1,0/1),(0/1,0/1,0/1))  0.010952  0.056314  x  0.016905  0.028548  0.027855
1 32 points_circle_32.txt
32 ((0/1,0/1,0/1),(0/1,0/1,0/1))  0.026821 0.203225
 x  0.046018  0.099752  0.114437
1 64 points_circle_64.txt
64 ((0/1,0/1,0/1),(0/1,0/1,0/1))  0.117463  0.784588  x  0.181588  0.396457  0.453168
1 128 points_circle_128.txt
128 ((0/1,0/1,0/1),(0/1,0/1,0/1))  0.489025  3.17453  x  0.769714  1.60429 1.87088
1 256 points_circle_256.txt
256 ((0/1,0/1,0/1),(0/1,0/1,0/1))  2.2189  12.8307  x  3.30405  6.58824  7.65752
1 512 points_circle_512.txt
512 ((0/1,0/1,0/1),(0/1,0/1,0/1))  10.5987  60.1585  x  18.7911  28.7101  33.5745
1 1024 points_circle_1024.txt
1024 ((0/1,0/1,0/1),(0/1,0/1,0/1))  59.3485  307.348  x  93.2172  137.297  156.172

The archive file points_circle.tar.gz contains all the files in the table above and a short script that runs max_cover on them.

The second type of input consists of sets, such that each set contains n antipodal pairs of points on a circle of radius 1.

Square Radius n Input Points Covered Points Disc Center Efi
Ilia & Tzvika Michael & Omri Mor Omer & Gal Shahaf & Yair
1 3 points_circle_ap_6.txt 6 ((0/1,0/1,0/1),(0/1,0/1,0/1)) 0.003574
 x 0.006889
0.005411 0.008799
1 4
points_circle_ap_8.txt 8 ((0/1,0/1,0/1),(0/1,0/1,0/1)) 0.009011
 x 0.005773
1 8
points_circle_ap_16.txt 16 ((0/1,0/1,0/1),(0/1,0/1,0/1)) 0.01434
 x 0.017916
1 16 points_circle_ap_32.txt 32 ((0/1,0/1,0/1),(0/1,0/1,0/1)) 0.034638
 x 0.065083
1 32
points_circle_ap_64.txt 64 ((0/1,0/1,0/1),(0/1,0/1,0/1)) 0.088728
 x 0.253143
1 64 points_circle_ap_128.txt 128 ((0/1,0/1,0/1),(0/1,0/1,0/1)) 0.394491
 x 1.08274
1 128
points_circle_ap_256.txt 256 ((0/1,0/1,0/1),(0/1,0/1,0/1)) 1.56419
 x 4.31515
1 256 points_circle_ap_512.txt 512 ((0/1,0/1,0/1),(0/1,0/1,0/1)) 7.6899
 x 19.7682
1 512
points_circle_ap_1024.txt 1024 ((0/1,0/1,0/1),(0/1,0/1,0/1)) 43.9682
 x 103.305
1 1024
points_circle_ap_2048.txt 2048 ((0/1,0/1,0/1),(0/1,0/1,0/1)) 342.062

 x  491.836 724.608

The archive file points_circle_ap.tar.gz contains all the files in the table above and a short script that runs max_cover on them.

The forth type consists of a mixture of the above. That is, each set consists of nm subsets of points, where each subset contains k points on a circle of radius 1 centered on a grid point of an n x m grid. Each subset of k=a+2b points contains a non-antipodal points and 2b pairs of antipodal points.

Square Radius n m a b
Input Points Covered Points Disc Center Efi Ilia & Tzvika Michael & Omri Mor Omer & Gal Shahaf & Yair
1 3 3
points_mix_63.txt 18 ((553/425,-633/3034,277611/180625),(1033/850,-247/1517,277611/180625)) 0.064244
0.384814 v
0.100866 0.177599 0.19505
2 3
points_mix_63.txt 28 ((81/34,-1/5,15/1),(25/17,-1/10,15/1)) 0.087212
0.619766 v 0.143738 0.290456 0.31478
4 3
points_mix_63.txt 41 ((1531/850,-529/4621,1.00692e+07/180625),(729/425,-919/9242,1.00692e+07/180625)) 0.119665
v 0.196896 0.401929 0.475066
points_mix_63.txt 50 ((72/65,-63/778,6224/4225),(56/65,97/1556,6224/4225)) 0.125936
0.978487 v 0.226521 0.465974 0.540257
points_mix_63.txt 57
((1382/493,-355/1542,3.57012e+07/243049),(385/986,13/771,3.57012e+07/243049)) 0.142349
1.05042 v 0.248979

The archive file points_mix.tar.gz contains all the files in the table above and a short script that runs max_cover on them.

Installation and Assistance

You can obtain the latest CGAL public release from here.
The CGAL installation manual is available here.
Further detailed instructions for installing CGAL on Windows is available here. This is a bit outdated though.
Ipe is a drawing editor for creating figures in PDF or (encapsulated) Postscript format.


1. Question

I'm using an arrangement of type CGAL::Arrangement_2<Traits> where Traits is an object of type CGAL::Arr_circle_segment_traits_2<Kernel> and the kernel is of some rational kernel type. I'd like to find out whether a given point is in the closed disc bounded by a circle centered at another given point, where the points may be either of type Kernel::Point_2 or Traits::Point_2.


A point of type Kernel::Point_2 has rational coordinates. (The type of the coordinates is Kernel::FT.) An object of type Kernel::Bounded_side_2 has an operator that accepts a circle and a point and determines whether the point is inside the circle, outside, or on the circle boundary, where the point and circle center-point are both of type Kernel::Point_2.

A point in the arrangement may have algebraic coordinates that are not rational. Thus, the type Kernel::FT cannot be used to represent the coordinates of such a point in an exact manner. Instead, a special type is used---it is an instance of CGAL::Sqrt_extension; see In particular, the type of such points is Traits::Point_2 and the type of its coordinates is Traits::CoordNT.

Assume that the center of the circle of of type Traits::Point_2 and the query point is of type Kernel::Point_2, the function below can be used to achieve the task

bool is_on(const Traits::Point_2& c, const Kernel::Point_2& p, const Kernel::FT& square_radius)
  auto xd = c.x() - p.x();
  auto yd = c.y() - p.y();
  return (xd*xd + yd*yd) <= square_radius;

The function above can be converted into a template, parametrized with two point types. These types can be substituted with either Kernel::Point_2 and Traits:Point_2, respectively, or Traits:Point_2 and Kernel::Point_2, respectively, when the template is instantiated. (I believe that it is possible to handle the case where the types of both, the circle center-point and the query point, is Traits:Point_2. However, it requires some modifications applied to the code above.)

If you only want to compute the intersection of two circles, you can use a circular kernel, e.g., an object of type CGAL::Exact_circular_kernel_2. The intersection point is of type Kernel::Circular_arc_point_2 and its coordinates are of type, you guessed it, an instance of CGAL::Sqrt_extension. The circle kernel has a nested type called Has_on_2, which can be used to determine whether a query point is inside a circle, outside, or on its boundary. Unfortunately, the circle center-point cannot be of type Kernel::Circular_arc_point_2. However, the template above can be used to do the job as well (as long as the query point is of type Kernel::Point_2).

Disclaimer: I haven't used any of the above in my implementation.

2. Question

How should I print the center point?


Use the exporters bellow. As indicated above, they are also defined in the sample program.

/*! Export a coordinate of type square-root extension to the standard output.
std::ostream& operator<<(std::ostream& os, const CoordNT& c)
  CGAL::Rational_traits<Number_type> rt;
  std::cout << "("
            << rt.numerator(c.a0()) << "/" << rt.denominator(c.a0()) << ","
            << rt.numerator(c.a1()) << "/" << rt.denominator(c.a1()) << ","
            << rt.numerator(c.root()) << "/" << rt.denominator(c.root()) << ")";
  return os;

/*! Export a point of two coordinates of type square-root extension to the
 * standard output.
std::ostream& operator<<(std::ostream& os, const Point& p)
  std::cout << "(" << p.x() << "," << p.y() << ")";
  return os;

3. Question

When I use those functions, I get the correct points, but with a different representation.


We are using exact arithmetic, so 1/2 is equivalent to 2/4, and a Sqrt_extension (0/1,0/1,0/1) (which is 0) is equivalent to (0/1, -1/1,1/1). There is a function in CGAL that can be used to compare two such numbers. Better yet, you can use the functor Traits::Equal_2, where Traits is an instance of CGAL::Arr_circle_segment_traits_2, to determine whether 2 points with coordinates of type CGAL::Sqrt_extension are equal, or Kernel::Equal_2, where Kernel is an instance of CGAL::Exact_circular_kernel_2. In general, in order to determine whether n0 = a0 + b0 * sqrt(c0) is equivalent to n1 = a1 + b1 * sqrt(c1), you would need to shuffle around the components and apply a power-of-2 twice.

4. Question

The example program fails to compile on Nova.


The latest version of CGAL (namely 4.9) has been installed under /usr/local/lib/cgal_4.9/.
You can either run
cmake -DCGAL_DIR=/usr/local/lib/cgal_4.9/lib/CGAL .
or set the environment variable
export CGAL_DIR=/usr/local/lib/cgal_4.9/lib
and then run cmake . as usual.

If your program will eventually compile and run on Windows, but fail on Nova, we will accept it. (We will probably ask you to apply some minor changes and re-submit without any penalty.)

5. Question

How can I map arrangement cells to indices?


It is possible to extend the Vertex, halfedge, or face cells of the arrangement data structure. As a mater of fact, when you extend them, you extend the corresponding cells of the underlying halfedge data-structure. If you want to extend all the 3 types, you can use the Arr_extended_dcel template; see

Say you want to extend the vertex, halfedge, and face with size_t, bool, and size_t, respectively. Then, you can define the following types:

#include <CGAL/Arr_extended_dcel.h>

typedef CGAL::Arr_extended_dcel<Traits, size_t, bool, size_t>   Dcel;
typedef CGAL::Arrangement_2<Traits, Dcel>                       Arrangement;

You can set and get the data stored using the member functions .data() and .set_data() of the appropriate cell type, respectively.

For example, the statement bellow sets the data of the vertex pointed by a vertex handle vh:


where x is of the type that was used to extend the vertex (size_t in the example above).

You can always create an external map from vertex handles to indices. This is, naturally, less efficient time-wise, but saves space...

Efi Fogel hompage

Document Actions