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.

General

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:

  if (CMAKE_COMPILER_IS_GNUCXX)
    add_definitions(-std=c++11)
  endif()

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 http://doc.cgal.org/latest/Number_types/classCGAL_1_1Sqrt__extension.html) is used to represent such numbers. It supports the arithmetic operations '+', '-', '*', and '/'.

Contest

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 n² points on an n x n grid.

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.

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.

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.

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.

Q&A

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.

Answer

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 http://doc.cgal.org/latest/Number_types/classCGAL_1_1Sqrt__extension.html. 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?

Answer

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.

Answer

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.

Answer

cmake -DCGAL_DIR=/usr/local/lib/cgal_4.9/lib/CGAL .

export CGAL_DIR=/usr/local/lib/cgal_4.9/lib

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?

Answer

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 http://doc.cgal.org/latest/Arrangement_on_surface_2/classCGAL_1_1Arr__extended__dcel.html.

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:

vh->set_data(x);

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