Merge pull request #167 from mgastner/fix-concentric-circle-float-issue

Adjust Grid Bug Fix, Other Code Improvements

.gitignore

@@ -38,6 +38,9 @@ cartogram
 # Ignore test logs
 *.txt
 
+# But, don't ignore CMakeLists.txt
+!CMakeLists.txt
+
 # Boilerplate code for .gitignore
 
 # Prerequisites
@@ -73,13 +76,18 @@ cartogram
 *.out
 *.app
 
+# Ignore backup files
+*.backup
+
 # Ignore CMakeFiles directory
 CMakeFiles
+cmake
 cmake-build-*/
 bin
 cmake_install.cmake
 CTestTestfile.cmake
 CMakeCache.txt
 Makefile
 .idea
-.Rproj.user
+.Rproj.user
+.vscode-ctags

.vscode/tasks.json

@@ -0,0 +1,16 @@
+{
+	"version": "2.0.0",
+	"tasks": [
+		{
+			"type": "cmake",
+			"label": "CMake: install",
+			"command": "install",
+			"problemMatcher": [],
+			"detail": "CMake template install task",
+			"group": {
+				"kind": "build",
+				"isDefault": true
+			}
+		}
+	]
+}

CMakeLists.txt

@@ -4,6 +4,7 @@ project(cartogram LANGUAGES CXX)
 # ========== Project Setup ==========
 set(CMAKE_CXX_STANDARD 20)
 set(THREADS_PREFER_PTHREAD_FLAG ON)
+set(CMAKE_COLOR_DIAGNOSTICS ON)
 set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
 
 # Default build type
@@ -24,21 +25,6 @@ find_package(PkgConfig REQUIRED)
 pkg_search_module(fftw REQUIRED fftw3 IMPORTED_TARGET)
 pkg_search_module(CAIRO REQUIRED CAIRO IMPORTED_TARGET)
 
-# ========== Compiler Setup ==========
-if(APPLE)
-  set(LLVM_BASE_PATH "/usr/local/opt/llvm@17/bin/")
-
-  if(${CMAKE_HOST_SYSTEM_PROCESSOR} STREQUAL "arm64")
-    set(LLVM_BASE_PATH "/opt/homebrew/opt/llvm@17/bin/")
-  endif()
-
-  set(CMAKE_CXX_COMPILER "${LLVM_BASE_PATH}clang++")
-  set(CMAKE_C_COMPILER "${LLVM_BASE_PATH}clang")
-elseif(UNIX)
-  set(CMAKE_CXX_COMPILER "g++-13")
-  set(CMAKE_C_COMPILER "gcc-13")
-endif()
-
 # ========== Source Files ==========
 file(GLOB_RECURSE CARTOGRAM_SOURCES "src/*.cpp")
 add_executable(cartogram ${CARTOGRAM_SOURCES})

README.md

@@ -94,12 +94,12 @@ The CSV file should be in the following format:
 
 ### Testing
 
-If you'd like to contribute to the project, please run our tests after you make any changes. 
+If you'd like to contribute to the project, please run our tests after you make any changes.
 
 To run the unit tests, execute the following command:
 
     ctest --verbose
-    
+
 To learn more about the tests, you may go to the `cartogram-cpp/tests` directory and read the `README.md` file.
 
 Additionally, you may go to the `cartogram-cpp/tests` directory and run the following command:

include/cartogram_info.hpp

@@ -37,7 +37,7 @@ class CartogramInfo
   void read_geojson(const std::string &, bool, std::string &);
   std::map<std::string, InsetState> &ref_to_inset_states();
   void replace_missing_and_zero_target_areas();
-  void set_map_name(const std::string &);
+  std::string set_map_name(const std::string &);
   void shift_insets_to_target_position();
   void write_geojson(const std::string &, const std::string &, bool);
 };

include/cgal_typedef.hpp

@@ -32,7 +32,7 @@ typedef CGAL::Orthtree<CGAL::Orthtree_traits_2<Scd>, std::vector<Point>>
   Quadtree;
 
 // Delaunay triangulation
-typedef CGAL::Delaunay_triangulation_2<Scd> Delaunay;
+typedef CGAL::Constrained_Delaunay_triangulation_2<Scd> Delaunay;
 typedef Delaunay::Line_face_circulator Line_face_circulator;
 typedef Delaunay::Face_handle Face_handle;
 

include/constants.hpp

@@ -4,14 +4,15 @@
 #include <limits>
 #include <numbers>
 
-constexpr unsigned int default_long_grid_length = 128;
-constexpr unsigned int max_allowed_grid_length = 4096;
-constexpr unsigned int default_grid_factor = 4;
+constexpr unsigned int default_long_grid_length = 256;
+constexpr unsigned int max_allowed_autoscale_grid_length = 1024;
+constexpr unsigned int default_grid_factor = 2;
 constexpr double dbl_epsilon = std::numeric_limits<double>::epsilon();
 constexpr double dbl_inf = std::numeric_limits<double>::infinity();
 constexpr double dbl_resolution = 1e-8;
 constexpr unsigned int max_integrations = 100;
 constexpr double max_permitted_area_error = 0.01;
+constexpr double max_permitted_area_drift = 0.05;
 constexpr double padding_unless_world = 1.5;
 constexpr double pi = std::numbers::pi;
 constexpr double earth_surface_area = 510.1e6;

include/inset_state.hpp

@@ -12,6 +12,14 @@
 struct max_area_error_info {
   double value;
   std::string geo_div;
+
+  // Conversion to std::tuple to enable structured bindings
+  // this allows:
+  // auto [value, geo_div] = max_area_error();
+  operator std::tuple<double, std::string>() const
+  {
+    return std::make_tuple(value, geo_div);
+  }
 };
 
 struct proj_qd {  // quadtree-delaunay projection
@@ -23,9 +31,14 @@ class InsetState
 {
 private:
   std::unordered_map<std::string, double> area_errors_;
+
   std::unordered_set<Point> unique_quadtree_corners_;
   proj_qd proj_qd_;
-  std::vector<proj_qd> proj_sequence_;
+
+  // Failed constraints
+  std::vector<Segment> failed_constraints_dt_projected_;
+  std::vector<Segment> failed_constraints_;
+  Delaunay og_dt_;
 
   // Bounding boxes of Quadtree cells
   std::vector<Bbox> quadtree_bboxes_;
@@ -35,9 +48,7 @@ class InsetState
   std::unordered_map<std::string, Color> colors_;
 
   // Unblurred density mean, min, max
-  double dens_min_;
-  double dens_mean_;
-  double dens_max_;
+  double dens_min_, dens_mean_, dens_max_, exterior_density_;
 
   // Scaling factor to convert equal-area-projection unit to lx*ly unit.
   double latt_const_;
@@ -49,9 +60,16 @@ class InsetState
   // Geographic divisions in this inset
   std::vector<GeoDiv> geo_divs_;
 
+  // Create a map from GeoDiv ID to index in geo_divs_
+  std::map<std::string, size_t> geo_divs_id_to_index_;
+
   // Copy of original data
   std::vector<GeoDiv> geo_divs_original_;
 
+  // Copy to original data to be transform to keep the same points in the final
+  // cartogram
+  std::vector<GeoDiv> geo_divs_original_transformed_;
+
   // Chosen diagonal for each grid cell
   boost::multi_array<int, 2> grid_diagonals_;
 
@@ -67,6 +85,7 @@ class InsetState
   std::unordered_map<std::string, bool> is_input_target_area_missing_;
   std::unordered_map<std::string, std::string> labels_;
   unsigned int lx_{}, ly_{};  // Lattice dimensions
+  unsigned int n_fails_during_flatten_density_;
   unsigned int n_finished_integrations_;
   std::string pos_;  // Position of inset ("C", "T" etc.)
   boost::multi_array<Point, 2> proj_;  // Cartogram projection
@@ -91,18 +110,20 @@ class InsetState
   void apply_smyth_craster_projection();
 
   // Calculate difference between initial area and current area
-  double area_drift() const;
+  double area_expansion_factor() const;
   double area_error_at(const std::string &) const;
   void auto_color();  // Automatically color GeoDivs
   Bbox bbox(bool = false) const;
   void blur_density(double, bool);
   double blur_width() const;
-  void check_topology();
+  void check_topology() const;
   int chosen_diag(const Point v[4], unsigned int &, bool = false) const;
   Color color_at(const std::string &) const;
   bool color_found(const std::string &) const;
   bool colors_empty() const;
   unsigned int colors_size() const;
+  bool continue_integrating() const;
+  void create_and_store_quadtree_cell_corners();
   void create_contiguity_graph(unsigned int);
   void create_delaunay_t();
   void densify_geo_divs();
@@ -116,10 +137,11 @@ class InsetState
   void fill_grid_diagonals(bool = false);
 
   // Density functions
-  void fill_with_density(bool);  // Fill map with density, using scanlines
+  void fill_with_density_rays(bool);  // Fill map with density, using scanlines
+  void fill_with_density_clip(bool);  // Fill map with density, using clipping
   void flatten_density();  // Flatten said density with integration
   void flatten_ellipse_density();
-  void flatten_density_with_node_vertices();
+  bool flatten_density_with_node_vertices();
 
   const std::vector<GeoDiv> &geo_divs() const;
   std::vector<std::vector<Color>> grid_cell_colors(unsigned int cell_width);
@@ -138,24 +160,29 @@ class InsetState
     unsigned int y,
     unsigned int cell_width);
   double grid_cell_area_km(const unsigned int i, const unsigned int j);
-  void holes_inside_polygons();
+  void holes_inside_polygons() const;
   double grid_cell_target_area_per_km(
     const unsigned int i,
     const unsigned int j,
     const double total_target_area,
     const double total_inset_area);
   Bbox get_bbox_bar(const double bar_width, const double bar_height);
 
+  GeoDiv &get_geo_div(const std::string &);
+
   std::pair<double, unsigned int> get_km_legend_length();
   std::pair<double, unsigned int> get_visual_variable_legend_length();
 
+  void increment_n_fails_during_flatten_density();
   void increment_integration();
   double initial_area() const;
   double initial_target_area() const;
   void initialize_cum_proj();
   void initialize_identity_proj();
   void insert_color(const std::string &, const Color &);
   void insert_color(const std::string &, std::string &);
+  bool insert_constraint_safely(const Point &p1, const Point &p2);
+  bool insert_constraint_safely_to_dt(Delaunay &dt, const Point &p1, const Point &p2);
   void insert_label(const std::string &, const std::string &);
   void insert_target_area(const std::string &, double);
   void insert_whether_input_target_area_is_missing(const std::string &, bool);
@@ -166,6 +193,7 @@ class InsetState
     char,
     unsigned int) const;
   bool is_input_target_area_missing(const std::string &) const;
+  void is_simple(const char* caller_func) const;
   std::string label_at(const std::string &) const;
   double latt_const() const;
   unsigned int lx() const;
@@ -179,6 +207,7 @@ class InsetState
   std::pair<Point, Point> max_and_min_grid_cell_area_index(
     unsigned int cell_width);
   unsigned int n_finished_integrations() const;
+  unsigned int n_fails_during_flatten_density() const;
   unsigned int n_geo_divs() const;
   unsigned long n_points() const;
   unsigned int n_rings() const;
@@ -189,9 +218,8 @@ class InsetState
   Point projected_point(const Point &, bool = false) const;
   Point projected_point_with_triangulation(const Point &, bool = false) const;
   void project_with_cum_proj();
-  void project_with_delaunay_t();
+  void project_with_delaunay_t(bool);
   void project_with_triangulation();
-  void project_with_proj_sequence();
   void push_back(const GeoDiv &);
   FTReal2d &ref_to_fluxx_init();
   FTReal2d &ref_to_fluxy_init();
@@ -202,7 +230,6 @@ class InsetState
   void replace_target_area(const std::string &, double);
   void rescale_map(unsigned int, bool);
   void revert_smyth_craster_projection();
-  void rings_are_simple();
   void set_area_errors();
   void set_grid_dimensions(unsigned int, unsigned int);
   void set_geo_divs(std::vector<GeoDiv> new_geo_divs);
@@ -225,22 +252,22 @@ class InsetState
   std::array<Point, 3> untransformed_triangle(const Point &, bool = false)
     const;
   void trim_grid_heatmap(cairo_t *cr, double padding);
+  void update_delaunay_t();
 
   // Cairo functions
   void write_cairo_map(
     const std::string &,
     bool,
-    std::unordered_map<Point, Point> = std::unordered_map<Point, Point>());
+    const std::unordered_map<Point, Vector> = std::unordered_map<Point, Vector>()) const;
   void write_cairo_polygons_to_svg(
     const std::string &,
     bool,
     bool,
     bool,
-    std::unordered_map<Point, Point> &);
+    const std::unordered_map<Point, Vector> &) const;
 
-  void write_delaunay_triangles(const std::string &);
+  void write_delaunay_triangles(const std::string &, const bool);
   void write_grid_heatmap_data(const std::string filename);
-
   void write_grid_heatmap_image(
     const std::string filename,
     const bool plot_equal_area_map,
@@ -254,7 +281,7 @@ class InsetState
     cairo_t *cr,
     const bool fill_polygons,
     const bool colors,
-    const bool plot_equal_area_map);
+    const bool plot_equal_area_map) const;
   void write_labels_on_surface(cairo_t *cr);
   void write_density_image(
     const std::string filename,
@@ -266,4 +293,4 @@ class InsetState
   void write_quadtree(const std::string &);
 };
 
-#endif // INSET_STATE_HPP_
+#endif  // INSET_STATE_HPP_

include/intersection.hpp

@@ -14,6 +14,7 @@ class intersection
 public:
   double target_density{};  // GeoDiv's target_density
   std::string geo_div_id;  // GeoDIv's ID
+  size_t pwh_idx;  // Polygon ID
   bool ray_enters{};  // Does the ray enter (true) a GeoDiv or exit (false)?
 
   // Overloading "<" operator, similar to above
@@ -40,6 +41,7 @@ void add_intersections(
   double,
   double,
   const std::string &,
+  size_t,
   char);
 
 #endif // INTERSECTION_HPP_

include/parse_arguments.hpp

@@ -24,6 +24,7 @@ argparse::ArgumentParser parsed_arguments(
   bool &plot_polygons,
   bool &remove_tiny_polygons,
   double &minimum_polygon_area,
-  bool &plot_quadtree);
+  bool &plot_quadtree,
+  bool &rays);
 
 #endif // PARSE_ARGUMENTS_HPP_

include/progress_tracker.hpp

@@ -26,18 +26,7 @@ class ProgressTracker
 private:
   double total_geo_divs_;  // Total number of GeoDivs to monitor progress
   double progress_;  // Progress measured on a scale from 0 (start) to 1 (end)
-  indicators::ProgressBar bar_{
-    indicators::option::BarWidth{75},
-    indicators::option::Start{"["},
-    indicators::option::Fill{"■"},
-    indicators::option::Lead{"■"},
-    indicators::option::Remainder{"-"},
-    indicators::option::End{"]"},
-    indicators::option::FontStyles{
-      std::vector<indicators::FontStyle>{indicators::FontStyle::bold}},
-    indicators::option::ShowPercentage{true},
-    indicators::option::ShowElapsedTime{true},
-    indicators::option::Stream{std::cerr}};
+  indicators::ProgressBar bar_;
 };
 
 #endif  // PROGRESS_TRACKER_H

include/round_point.hpp

@@ -11,5 +11,6 @@ double rounded_to_bicimal(
   const unsigned int,
   const unsigned int);
 Point rounded_point(const Point &, const unsigned int, const unsigned int);
+Point rounded_point(const Point &, const unsigned int);
 
 #endif // ROUND_POINT_HPP_