verfiygrid: changed datatype of 2D point to Point.

src/Samplegridicon.cc

@@ -316,10 +316,10 @@ compute_child_from_bounds(CellIndex *cellindex2, long ncells2, double *grid_cent
   size_t *nbr_addr = &knnWeights.m_addr[0];
 
   int ncorner = 3;
-  double cell_corners_xyz[12];
-  double cell_corners_plane_projection[8];
   double center_point_xyz[3];
-  double center_point_plane_projection[2];
+  double cell_corners_xyz[12];
+  Point center_point_plane_projection;
+  Varray<Point> cell_corners_plane_projection(4);
 
   long *child2 = (long *) Malloc(MAX_CHILDS * ncells2 * sizeof(long));
   cellindex2->child = child2;

src/Verifygrid.cc

@@ -171,7 +171,7 @@ find_coordinate_to_ignore(const double *cell_corners_xyz)
 }
 
 static inline double
-is_point_left_of_edge(const double (&point1)[2], const double (&point2)[2], const double (&point)[2])
+is_point_left_of_edge(const Point &point1, const Point &point2, const Point &point)
 {
   /*
      Computes whether a point is left of the line through point1 and point2.
@@ -184,11 +184,11 @@ is_point_left_of_edge(const double (&point1)[2], const double (&point2)[2], cons
      This algorithm is by Dan Sunday (geomalgorithms.com) and is completely free for use and modification.
   */
 
-  return ((point2[0] - point1[0]) * (point[1] - point1[1]) - (point[0] - point1[0]) * (point2[1] - point1[1]));
+  return ((point2.x - point1.x) * (point.y - point1.y) - (point.x - point1.x) * (point2.y - point1.y));
 }
 
 int
-winding_numbers_algorithm(const double *cell_corners, int number_corners, const double (&point)[2])
+winding_numbers_algorithm(const Varray<Point> &cell_corners, int number_corners, const Point &point)
 {
   /*
      Computes whether a point is inside the bounds of a cell. This is the solution to the point in polygon problem.
@@ -198,27 +198,29 @@ winding_numbers_algorithm(const double *cell_corners, int number_corners, const
 
   int winding_number = 0;
 
-  for (int i = 0; i < number_corners - 1; i++)
+  for (int k = 0; k < number_corners - 1; k++)
     {
-      if (cell_corners[i * 2 + 1] <= point[1])
+      if (cell_corners[k].y <= point.y)
         {
-          if (cell_corners[(i + 1) * 2 + 1] > point[1])
+          if (cell_corners[k + 1].y > point.y)
             {
-              const double point1[2] = { cell_corners[i * 2 + 0], cell_corners[i * 2 + 1] };
-              const double point2[2] = { cell_corners[(i + 1) * 2 + 0], cell_corners[(i + 1) * 2 + 1] };
+              const auto &point1 = cell_corners[k];
+              const auto &point2 = cell_corners[k + 1];
 
               if (is_point_left_of_edge(point1, point2, point) > 0) winding_number++;
+              // printf("    1: %d %d %g %g %g %g %g %g\n", k, winding_number, point1.x, point1.y, point2.x, point2.y, point.x, point.y);
             }
         }
       else
         {
-          if (cell_corners[(i + 1) * 2 + 1] <= point[1])
+          if (cell_corners[k + 1].y <= point.y)
             {
-              const double point1[2] = { cell_corners[i * 2 + 0], cell_corners[i * 2 + 1] };
-              const double point2[2] = { cell_corners[(i + 1) * 2 + 0], cell_corners[(i + 1) * 2 + 1] };
+              const auto &point1 = cell_corners[k];
+              const auto &point2 = cell_corners[k + 1];
 
               if (is_point_left_of_edge(point1, point2, point) < 0) winding_number--;
-            }
+              //  printf("    2: %d %d %g %g %g %g %g %g\n", k, winding_number, point1.x, point1.y, point2.x, point2.y, point.x, point.y);
+           }
         }
     }
 
@@ -234,7 +236,7 @@ sign(double x)
 }
 
 static bool
-is_simple_polygon_convex(const Varray<double> &cell_corners, int number_corners)
+is_simple_polygon_convex(const Varray<Point> &cell_corners, int number_corners)
 {
   // Tests in which direction the polygon winds when walking along its edges. Does so for all edges of the polygon.
 
@@ -242,10 +244,8 @@ is_simple_polygon_convex(const Varray<double> &cell_corners, int number_corners)
 
   for (int i = 0; i < number_corners - 2; i++)
     {
-      double turns_to = (cell_corners[i * 2 + 0] - cell_corners[(i + 1) * 2 + 0])
-                            * (cell_corners[(i + 1) * 2 + 1] - cell_corners[(i + 2) * 2 + 1])
-                        - (cell_corners[i * 2 + 1] - cell_corners[(i + 1) * 2 + 1])
-                              * (cell_corners[(i + 1) * 2 + 0] - cell_corners[(i + 2) * 2 + 0]);
+      auto turns_to = (cell_corners[i].x - cell_corners[i + 1].x) * (cell_corners[i + 1].y - cell_corners[i + 2].y)
+                    - (cell_corners[i].y - cell_corners[i + 1].y) * (cell_corners[i + 1].x - cell_corners[i + 2].x);
 
       // In the first iteration the direction of winding of the entire polygon is set. Better not be 0.
 
@@ -264,8 +264,9 @@ is_simple_polygon_convex(const Varray<double> &cell_corners, int number_corners)
   return true;
 }
 
+
 double
-calculate_the_polygon_area(const double *cell_corners, int number_corners)
+calculate_the_polygon_area(const Varray<Point> &cell_corners, int number_corners)
 {
   /* This algorithm is based on the calculation from Wolfram Mathworld Polygon Area. It results in the area of planar
    * non-self-intersecting polygon. */
@@ -273,8 +274,7 @@ calculate_the_polygon_area(const double *cell_corners, int number_corners)
   double twice_the_polygon_area = 0.0;
 
   for (int i = 0; i < number_corners - 1; i++)
-    twice_the_polygon_area
-        += (cell_corners[i * 2 + 0] * cell_corners[(i + 1) * 2 + 1]) - (cell_corners[(i + 1) * 2 + 0] * cell_corners[i * 2 + 1]);
+    twice_the_polygon_area += (cell_corners[i].x * cell_corners[i + 1].y) - (cell_corners[i + 1].x * cell_corners[i].y);
 
   return twice_the_polygon_area / 2.0;
 }
@@ -397,49 +397,42 @@ void set_cell_corners_xyz(int ncorner, const double *cell_corners_lon, const dou
   cell_corners_xyz[ncorner * 3 + 2] = cell_corners_xyz[2];
 }
 
-void set_center_point_plane_projection(int coordinate_to_ignore, const double (&center_point_xyz)[3], double (&center_point_plane_projection)[2])
+void set_center_point_plane_projection(int coordinate_to_ignore, const double (&center_point_xyz)[3], Point &center_point_plane_projection)
 {
+  // clang-format off
   switch (coordinate_to_ignore)
     {
-    case 1:
-      center_point_plane_projection[0] = center_point_xyz[1];
-      center_point_plane_projection[1] = center_point_xyz[2];
-      break;
-    case 2:
-      center_point_plane_projection[0] = center_point_xyz[2];
-      center_point_plane_projection[1] = center_point_xyz[0];
-      break;
-    case 3:
-      center_point_plane_projection[0] = center_point_xyz[0];
-      center_point_plane_projection[1] = center_point_xyz[1];
-      break;
+    case 1: center_point_plane_projection = Point {center_point_xyz[1], center_point_xyz[2]}; break;
+    case 2: center_point_plane_projection = Point {center_point_xyz[2], center_point_xyz[0]}; break;
+    case 3: center_point_plane_projection = Point {center_point_xyz[0], center_point_xyz[1]}; break;
     }
+  // clang-format on
 }
 
 void set_cell_corners_plane_projection(int coordinate_to_ignore, int ncorner, const double *cell_corners_xyz,
-                                       double *cell_corners_plane_projection)
+                                       Varray<Point> &cell_corners_plane_projection)
 {
   switch (coordinate_to_ignore)
     {
     case 1:
       for (int corner_no = 0; corner_no <= ncorner; corner_no++)
         {
-          cell_corners_plane_projection[corner_no * 2 + 0] = cell_corners_xyz[corner_no * 3 + 1];
-          cell_corners_plane_projection[corner_no * 2 + 1] = cell_corners_xyz[corner_no * 3 + 2];
+          cell_corners_plane_projection[corner_no].x = cell_corners_xyz[corner_no * 3 + 1];
+          cell_corners_plane_projection[corner_no].y = cell_corners_xyz[corner_no * 3 + 2];
         }
       break;
     case 2:
       for (int corner_no = 0; corner_no <= ncorner; corner_no++)
         {
-          cell_corners_plane_projection[corner_no * 2 + 0] = cell_corners_xyz[corner_no * 3 + 2];
-          cell_corners_plane_projection[corner_no * 2 + 1] = cell_corners_xyz[corner_no * 3 + 0];
+          cell_corners_plane_projection[corner_no].x = cell_corners_xyz[corner_no * 3 + 2];
+          cell_corners_plane_projection[corner_no].y = cell_corners_xyz[corner_no * 3 + 0];
         }
       break;
     case 3:
       for (int corner_no = 0; corner_no <= ncorner; corner_no++)
         {
-          cell_corners_plane_projection[corner_no * 2 + 0] = cell_corners_xyz[corner_no * 3 + 0];
-          cell_corners_plane_projection[corner_no * 2 + 1] = cell_corners_xyz[corner_no * 3 + 1];
+          cell_corners_plane_projection[corner_no].x = cell_corners_xyz[corner_no * 3 + 0];
+          cell_corners_plane_projection[corner_no].y = cell_corners_xyz[corner_no * 3 + 1];
         }
       break;
     }
@@ -465,12 +458,11 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
 
      The results of the tests are printed on stdout.
   */
-
   double center_point_xyz[3];
   Varray<double> cell_corners_xyz_open_cell(3 * ncorner);
 
   double corner_coordinates[3];
-  double center_point_plane_projection[2];
+  Point center_point_plane_projection;
 
   size_t no_of_cells_with_duplicates = 0;
   size_t no_usable_cells = 0;
@@ -496,7 +488,7 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
   std::vector<bool> marked_duplicate_indices(ncorner);
   Varray<double> cell_corners_xyz_without_duplicates(3 * ncorner);
   Varray<double> cell_corners_xyz(3 * (ncorner + 1));
-  Varray<double> cell_corners_plane_projection(2 * (ncorner + 1));
+  Varray<Point> cell_corners_plane_projection(ncorner + 1);
   Varray<size_t> cells_with_duplicates;
   cells_with_duplicates.reserve(gridsize);
 
@@ -539,7 +531,7 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
           if (IS_EQUAL(cell_corners_xyz_open_cell[off + 0], cell_corners_xyz_open_cell[off2 + 0])
               && IS_EQUAL(cell_corners_xyz_open_cell[off + 1], cell_corners_xyz_open_cell[off2 + 1])
               && IS_EQUAL(cell_corners_xyz_open_cell[off + 2], cell_corners_xyz_open_cell[off2 + 2]))
-            actual_number_of_corners = actual_number_of_corners - 1;
+            actual_number_of_corners--;
           else
             break;
         }
@@ -654,7 +646,7 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
       if (cell_corners_xyz[coordinate_to_ignore - 1] < 0) invert_result = true;
 
       set_center_point_plane_projection(coordinate_to_ignore, center_point_xyz, center_point_plane_projection);
-      set_cell_corners_plane_projection(coordinate_to_ignore, actual_number_of_corners, cell_corners_xyz.data(), cell_corners_plane_projection.data());
+      set_cell_corners_plane_projection(coordinate_to_ignore, actual_number_of_corners, cell_corners_xyz.data(), cell_corners_plane_projection);
 
       // Checking for convexity of the cell.
 
@@ -671,7 +663,7 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
 
       // Checking the arrangement or direction of cell vertices.
 
-      const auto polygon_area = calculate_the_polygon_area(cell_corners_plane_projection.data(), actual_number_of_corners + 1);
+      const auto polygon_area = calculate_the_polygon_area(cell_corners_plane_projection, actual_number_of_corners + 1);
       auto is_clockwise = are_polygon_vertices_arranged_in_clockwise_order(polygon_area);
 
       /* If the direction of the vertices was flipped during the projection onto the two-dimensional plane, the previous
@@ -694,9 +686,23 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
         {
           no_counterclockwise_cells += 1;
         }
-
+      /*
+      printf("ncorner actual_number_of_corners %d %d\n", ncorner, actual_number_of_corners);
+  printf("cell_no, is_clockwise, polygon_area %zu %d, %g\n", cell_no, is_clockwise, polygon_area);
+  printf("bounds %g %g   %g %g  %g %g  %g %g\n",
+         grid_corner_lon[cell_no * ncorner], grid_corner_lat[cell_no * ncorner],
+         grid_corner_lon[cell_no * ncorner+1], grid_corner_lat[cell_no * ncorner+1],
+         grid_corner_lon[cell_no * ncorner+2], grid_corner_lat[cell_no * ncorner+2],
+         grid_corner_lon[cell_no * ncorner+3], grid_corner_lat[cell_no * ncorner+3]);
+  printf("plane %g %g   %g %g  %g %g  %g %g    %g %g\n",
+         cell_corners_plane_projection[0].x, cell_corners_plane_projection[0].y,
+         cell_corners_plane_projection[1].x, cell_corners_plane_projection[1].y,
+         cell_corners_plane_projection[2].x, cell_corners_plane_projection[2].y,
+         cell_corners_plane_projection[3].x, cell_corners_plane_projection[3].y,
+         center_point_plane_projection.x, center_point_plane_projection.y);
+      */
       // The winding numbers algorithm is used to test whether the presumed center point is within the bounds of the cell.
-      auto winding_number = winding_numbers_algorithm(cell_corners_plane_projection.data(), actual_number_of_corners + 1,
+      auto winding_number = winding_numbers_algorithm(cell_corners_plane_projection, actual_number_of_corners + 1,
                                                       center_point_plane_projection);
 
       // if ( winding_number == 0 ) printf("%d,", cell_no+1);

src/verifygrid.h

@@ -17,12 +17,14 @@
 #ifndef VERIFYGRID_H
 #define VERIFYGRID_H
 
+struct Point { double x, y; };
+
 void set_cell_corners_xyz(int ncorner, const double *cell_corners_lon, const double *cell_corners_lat, double *cell_corners_xyz);
-void set_center_point_plane_projection(int coordinate_to_ignore, const double (&center_point_xyz)[3], double (&center_point_plane_projection)[2]);
-void set_cell_corners_plane_projection(int coordinate_to_ignore, int ncorner, const double *cell_corners_xyz, double *cell_corners_plane_projection);
+void set_center_point_plane_projection(int coordinate_to_ignore, const double (&center_point_xyz)[3], Point &center_point_plane_projection);
+void set_cell_corners_plane_projection(int coordinate_to_ignore, int ncorner, const double *cell_corners_xyz, Varray<Point> &cell_corners_plane_projection);
 int find_coordinate_to_ignore(const double *cell_corners_xyz);
-double calculate_the_polygon_area(const double *cell_corners, int number_corners);
+double calculate_the_polygon_area(const Varray<Point> &cell_corners, int number_corners);
 bool are_polygon_vertices_arranged_in_clockwise_order(double cell_area);
-int winding_numbers_algorithm(const double *cell_corners, int number_corners, const double (&point)[2]);
+int winding_numbers_algorithm(const Varray<Point> &cell_corners, int number_corners, const Point &point);
 
 #endif /* VERIFYGRID_H */