From 51012f9f84f45102c9277f73b4062abf25ad7681 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Tue, 9 May 2017 21:41:47 -0400 Subject: [PATCH] Add distance to polygon function Fix bug Add acknowledgements Add test Improve comment --- include/aspect/utilities.h | 11 ++++ source/utilities.cc | 67 ++++++++++++++++++++++++ tests/prm_distance_polygon.cc | 39 ++++++++++++++ tests/prm_distance_polygon.prm | 1 + tests/prm_distance_polygon/screen-output | 10 ++++ 5 files changed, 128 insertions(+) create mode 100644 tests/prm_distance_polygon.cc create mode 100644 tests/prm_distance_polygon.prm create mode 100644 tests/prm_distance_polygon/screen-output diff --git a/include/aspect/utilities.h b/include/aspect/utilities.h index 5006ea9c67f..22f7b03d998 100644 --- a/include/aspect/utilities.h +++ b/include/aspect/utilities.h @@ -140,6 +140,17 @@ namespace aspect polygon_contains_point(const std::vector > &point_list, const dealii::Point<2> &point); + /** + * Given a 2d point and a list of points which form a polygon, compute the smallest + * distance of the point to the polygon. The sign is negative for points outside of + * the polygon and positive for points inside the polygon. Note that the polygon + * is required to be convex for this function. + */ + template + double + signed_distance_to_polygon(const std::vector > &point_list, + const dealii::Point<2> &point); + /** * Given a vector @p v in @p dim dimensional space, return a set * of (dim-1) vectors that are orthogonal to @p v and to each diff --git a/source/utilities.cc b/source/utilities.cc index 3f2beb0adb3..725ee178322 100644 --- a/source/utilities.cc +++ b/source/utilities.cc @@ -312,6 +312,70 @@ namespace aspect return (wn != 0); } + template + double + signed_distance_to_polygon(const std::vector > &point_list, + const dealii::Point<2> &point) + { + // If the point lies outside polygon, we give it a negative sign, + // inside a positive sign. + const double sign = polygon_contains_point(point_list, point) ? 1.0 : -1.0; + + /** + * This code is based on http://geomalgorithms.com/a02-_lines.html#Distance-to-Infinite-Line, + * and therefore requires the following copyright notice: + * + * Copyright 2000 softSurfer, 2012 Dan Sunday + * This code may be freely used and modified for any purpose + * providing that this copyright notice is included with it. + * SoftSurfer makes no warranty for this code, and cannot be held + * liable for any real or imagined damage resulting from its use. + * Users of this code must verify correctness for their application. + * + */ + + const unsigned int n_poly_points = point_list.size(); + AssertThrow(n_poly_points >= 3, ExcMessage("Not enough polygon points were specified.")); + + // Initialize a vector of distances for each point of the polygon with a very large distance + std::vector distances(n_poly_points, 1e23); + + // Create another polygon but with all points shifted 1 position to the right + std::vector > shifted_point_list(n_poly_points); + shifted_point_list[0] = point_list[n_poly_points-1]; + + for (unsigned int i = 0; i < n_poly_points-1; ++i) + shifted_point_list[i+1] = point_list[i]; + + for (unsigned int i = 0; i < n_poly_points; ++i) + { + // Create vector along the polygon line segment + Tensor<1,2> vector_segment = shifted_point_list[i] - point_list[i]; + // Create vector from point to the second segment point + Tensor<1,2> vector_point_segment = point - point_list[i]; + + // Compute dot products to get angles + const double c1 = vector_point_segment * vector_segment; + const double c2 = vector_segment * vector_segment; + + // point lies closer to not-shifted polygon point, but perpendicular base line lies outside segment + if (c1 <= 0.0) + distances[i] = (Tensor<1,2> (point_list[i] - point)).norm(); + // point lies closer to shifted polygon point, but perpendicular base line lies outside segment + else if (c2 <= c1) + distances[i] = (Tensor<1,2> (shifted_point_list[i] - point)).norm(); + // perpendicular base line lies on segment + else + { + const Point<2> point_on_segment = point_list[i] + (c1/c2) * vector_segment; + distances[i] = (Tensor<1,2> (point - point_on_segment)).norm(); + } + } + + // Return the minimum of the distances of the point to all polygon segments + return *std::min_element(distances.begin(),distances.end()) * sign; + } + template std_cxx11::array,dim-1> orthogonal_vectors (const Tensor<1,dim> &v) @@ -1874,6 +1938,9 @@ namespace aspect template bool polygon_contains_point<2>(const std::vector > &pointList, const dealii::Point<2> &point); template bool polygon_contains_point<3>(const std::vector > &pointList, const dealii::Point<2> &point); + template double signed_distance_to_polygon<2>(const std::vector > &pointList, const dealii::Point<2> &point); + template double signed_distance_to_polygon<3>(const std::vector > &pointList, const dealii::Point<2> &point); + template std_cxx11::array,1> orthogonal_vectors (const Tensor<1,2> &v); template std_cxx11::array,2> orthogonal_vectors (const Tensor<1,3> &v); diff --git a/tests/prm_distance_polygon.cc b/tests/prm_distance_polygon.cc new file mode 100644 index 00000000000..f282a3565ab --- /dev/null +++ b/tests/prm_distance_polygon.cc @@ -0,0 +1,39 @@ +#include +#include +#include + +#include +#include + +int f() +{ + using namespace aspect::Utilities; + + const int dim=3; + + // A square polygon + std::vector< Point<2> > polygon(4); + polygon[0] = Point<2>(0.0,0.0); + polygon[1] = Point<2>(1.0,0.0); + polygon[2] = Point<2>(1.0,1.0); + polygon[3] = Point<2>(0.0,1.0); + + // Some points inside and outside the polygon + Point<2> points[] = {Point<2>(0.5,-1), Point<2>(0.5,0.5), Point<2>(0.001,0.2), Point<2>(2.0,2.0), Point<2>(0.001,0.75)}; + + std::cout << "Testing distance to polygon function with the following parameters: (polygon) " << polygon[0] << ", " << polygon[1] << ", " << polygon[2] << ", " << polygon[3] + << ", (points) " + << points[0] << ", " << points[1] << ", " << points[2] << ", " << points[3] << ", " + << points[4] << std::endl; + + for (unsigned int i = 0; i < 5; i++) + { + std::cout << "Minimal distance of point " << points[i] << " to polygon = " << signed_distance_to_polygon(polygon,points[i]) << std::endl; + } + + exit(0); + return 42; +} +// run this function by initializing a global variable by it +int i = f(); + diff --git a/tests/prm_distance_polygon.prm b/tests/prm_distance_polygon.prm new file mode 100644 index 00000000000..583a5b0da38 --- /dev/null +++ b/tests/prm_distance_polygon.prm @@ -0,0 +1 @@ +set Additional shared libraries = tests/libprm_distance_polygon.so diff --git a/tests/prm_distance_polygon/screen-output b/tests/prm_distance_polygon/screen-output new file mode 100644 index 00000000000..747e17f3ff1 --- /dev/null +++ b/tests/prm_distance_polygon/screen-output @@ -0,0 +1,10 @@ +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- + +Loading shared library <./libprm_distance_polygon.so> +Testing distance to polygon function with the following parameters: (polygon) 0 0, 1 0, 1 1, 0 1, (points) 0.5 -1, 0.5 0.5, 0.001 0.2, 2 2, 0.001 0.75 +Minimal distance of point 0.5 -1 to polygon = -1 +Minimal distance of point 0.5 0.5 to polygon = 0.5 +Minimal distance of point 0.001 0.2 to polygon = 0.001 +Minimal distance of point 2 2 to polygon = -1.41421 +Minimal distance of point 0.001 0.75 to polygon = 0.001