diff --git a/pyroomacoustics/libroom_src/geometry.cpp b/pyroomacoustics/libroom_src/geometry.cpp index 1c612278..26bdb098 100644 --- a/pyroomacoustics/libroom_src/geometry.cpp +++ b/pyroomacoustics/libroom_src/geometry.cpp @@ -233,20 +233,55 @@ Eigen::Vector3f cross(Eigen::Vector3f v1, Eigen::Vector3f v2) return v1.cross(v2); } +bool on_segment(const Eigen::Vector2f &c1, const Eigen::Vector2f &c2, const Eigen::Vector2f &p) { + // Given three collinear points c1, c2, p, the function checks if + // point p lies on line segment c1 <-> c2 + + float x_down, x_up, y_down, y_up; + x_down = fminf(c1.coeff(0), c2.coeff(0)); + x_up = fmaxf(c1.coeff(0), c2.coeff(0)); + y_down = fminf(c1.coeff(1), c2.coeff(1)); + y_up = fmaxf(c1.coeff(1), c2.coeff(1)); + return (x_down <= p.coeff(0) && p.coeff(0) <= x_up && y_down <= p.coeff(1) && p.coeff(1) <= y_up); +} -int is_inside_2d_polygon(const Eigen::Vector2f &p, +inline int is_left(const Eigen::Vector2f &P0, const Eigen::Vector2f &P1, const Eigen::Vector2f &P2) +{ + // isLeft(): tests if a point is Left|On|Right of an infinite line. + // on the line is defined as within epsilon + // Input: three points P0, P1, and P2 + // Return: >0 for P2 left of the line through P0 and P1 + // =0 for P2 on the line + // <0 for P2 right of the line + // See: Algorithm 1 "Area of Triangles and Polygons" + float test_value = ( + (P1.coeff(0) - P0.coeff(0)) * (P2.coeff(1) - P0.coeff(1)) + - (P2.coeff(0) - P0.coeff(0)) * (P1.coeff(1) - P0.coeff(1)) + ); + + if (fabsf(test_value) < libroom_eps) + return 0; + else if (test_value > 0) + return 1; + else + return -1; +} + +int is_inside_2d_polygon(const Eigen::Vector2f &p_test, const Eigen::Matrix &corners) { /* - Checks if a given point is inside a given polygon in 2D. + Checks if a given point is inside a given polygon in 2D using the winding + number method. This function checks if a point (defined by its coordinates) is inside a polygon (defined by an array of coordinates of its corners) by counting - the number of intersections between the borders and a segment linking - the given point with a computed point outside the polygon. - A boolean is also returned to indicate if a point is on a border of the - polygon (the point is still considered inside), which can be useful for - limit cases computations. + the number of left intersections between the edges and a half-line starting from + the test point. + + This is Dave Sunday's winding number algorithm described here: + http://profs.ic.uff.br/~anselmo/cursos/CGI/slidesNovos/Inclusion%20of%20a%20Point%20in%20a%20Polygon.pdf + It was modified to explicitely check for points on the edges of the polygon. p: (array size 2) coordinates of the point corners: (array size 2xN, N>2) coordinates of the corners of the polygon @@ -257,67 +292,35 @@ int is_inside_2d_polygon(const Eigen::Vector2f &p, 0 : the point is inside 1 : the point is on the boundary */ - - bool is_inside = false; // initialize point not in the polygon - int c1c2p, c1c2p0, pp0c1, pp0c2; int n_corners = corners.cols(); - - // find a point outside the polygon - int i_min; - corners.row(0).minCoeff(&i_min); - Eigen::Vector2f p_out; - p_out.resize(2); - p_out.coeffRef(0) = corners.coeff(0,i_min) - 1; - p_out.coeffRef(1) = p.coeff(1); - - // Now count intersections + int wn = 0; // the winding number counter + // loop through all edges of the polygon for (int i = 0, j = n_corners-1 ; i < n_corners ; j=i++) { - - // Check first if the point is on the segment - // We count the border as inside the polygon - c1c2p = ccw3p(corners.col(i), corners.col(j), p); - if (c1c2p == 0) + if (ccw3p(corners.col(j), corners.col(i), p_test) == 0 && on_segment(corners.col(j), corners.col(i), p_test)) { - // Here we know that p is co-linear with the two corners - float x_down, x_up, y_down, y_up; - x_down = fminf(corners.coeff(0,i), corners.coeff(0,j)); - x_up = fmaxf(corners.coeff(0,i), corners.coeff(0,j)); - y_down = fminf(corners.coeff(1,i), corners.coeff(1,j)); - y_up = fmaxf(corners.coeff(1,i), corners.coeff(1,j)); - if (x_down <= p.coeff(0) && p.coeff(0) <= x_up && y_down <= p.coeff(1) && p.coeff(1) <= y_up) - return 1; + // point is on the edge, we consider it inside + return 1; } - - // Now check intersection with standard method - c1c2p0 = ccw3p(corners.col(i), corners.col(j), p_out); - if (c1c2p == c1c2p0) // no intersection - continue; - - pp0c1 = ccw3p(p, p_out, corners.col(i)); - pp0c2 = ccw3p(p, p_out, corners.col(j)); - if (pp0c1 == pp0c2) // no intersection - continue; - - // at this point we are sure there is an intersection - - // the second condition takes care of horizontal edges and intersection on vertex - float c_max = fmaxf(corners.coeff(1,i), corners.coeff(1,j)); - if (p.coeff(1) + libroom_eps < c_max) - { - is_inside = !is_inside; + if (corners.coeff(1,j) <= p_test.coeff(1)) { // start y <= P.y + if (corners.coeff(1,i) > p_test.coeff(1)) { // an upward crossing + if (is_left(corners.col(j), corners.col(i), p_test) > 0) // P left of edge + ++wn; // have a valid up intersect + } + } else { // start y > P.y (no test needed) + if (corners.coeff(1, i) <= p_test.coeff(1)) { // a downward crossing + if (is_left(corners.col(j), corners.col(i), p_test) < 0) // P right of edge + --wn; // have a valid down intersect + } } - } - // for a odd number of intersections, the point is in the polygon - if (is_inside) - return 0; // point strictly inside + if (wn == 0) + return -1; else - return -1; // point is outside + return 0; } - float area_2d_polygon(const Eigen::Matrix &corners) { /* diff --git a/pyroomacoustics/libroom_src/room.cpp b/pyroomacoustics/libroom_src/room.cpp index 3f9431af..540ccc85 100644 --- a/pyroomacoustics/libroom_src/room.cpp +++ b/pyroomacoustics/libroom_src/room.cpp @@ -729,7 +729,7 @@ std::tuple < Vectorf, int, float > Room::next_wall_hit( Wall & w = scattered_ray ? walls[obstructing_walls[i]] : walls[i]; // To store the result of this iteration - Vectorf temp_hit; + Vectorf temp_hit = Vectorf::Zero(); // As a side effect, temp_hit gets a value (VectorXf) here int ret = w.intersection(start, end, temp_hit); @@ -948,6 +948,7 @@ void Room::simul_ray( auto p_hit = (1 - sqrt(1 - mic_radius_sq / std::max(mic_radius_sq, r_sq))); energy = transmitted / (r_sq * p_hit); // energy = transmitted / (travel_dist_at_mic - sqrtf(fmaxf(0.f, travel_dist_at_mic * travel_dist_at_mic - mic_radius_sq))); + microphones[k].log_histogram(travel_dist_at_mic, energy, start); } } diff --git a/pyroomacoustics/tests/tests_libroom/test_is_inside_2d_polygon.py b/pyroomacoustics/tests/tests_libroom/test_is_inside_2d_polygon.py index 8841eb15..3a6a8a63 100644 --- a/pyroomacoustics/tests/tests_libroom/test_is_inside_2d_polygon.py +++ b/pyroomacoustics/tests/tests_libroom/test_is_inside_2d_polygon.py @@ -30,18 +30,8 @@ polygons = [ np.array( [ # this one is clockwise - [ - 0, - 4, - 4, - 0, - ], - [ - 0, - 0, - 4, - 4, - ], + [0, 4, 4, 0], + [0, 0, 4, 4], ] ), np.array( diff --git a/pyroomacoustics/tests/tests_libroom/test_ray_energy.py b/pyroomacoustics/tests/tests_libroom/test_ray_energy.py index 107147cc..7189c92f 100644 --- a/pyroomacoustics/tests/tests_libroom/test_ray_energy.py +++ b/pyroomacoustics/tests/tests_libroom/test_ray_energy.py @@ -83,10 +83,10 @@ def test_square_room(self): print("Creating the python room (polyhedral)") walls_corners = [ - np.array([[0, 2, 2, 0], [0, 0, 0, 0], [0, 0, 2, 2]]), # left - np.array([[0, 0, 2, 2], [2, 2, 2, 2], [0, 2, 2, 0]]), # right - np.array([[0, 0, 0, 0], [0, 2, 2, 0], [0, 0, 2, 2]]), # front` - np.array([[2, 2, 2, 2], [0, 0, 2, 2], [0, 2, 2, 0]]), # back + np.array([[0, 2, 2, 0], [0, 0, 0, 0], [0, 0, 2, 2]]), # front + np.array([[0, 0, 2, 2], [2, 2, 2, 2], [0, 2, 2, 0]]), # back + np.array([[0, 0, 0, 0], [0, 2, 2, 0], [0, 0, 2, 2]]), # left` + np.array([[2, 2, 2, 2], [0, 0, 2, 2], [0, 2, 2, 0]]), # right np.array( [ [0, 2, 2, 0], @@ -157,8 +157,6 @@ def test_square_room(self): histogram_rt_poly = np.array(h_poly[0])[: len(histogram_gt)] histogram_rt_cube = np.array(h_cube[0])[: len(histogram_gt)] - print(abs(histogram_rt_poly - histogram_gt).max()) - print(abs(histogram_rt_cube - histogram_gt).max()) self.assertTrue(np.allclose(histogram_rt_poly, histogram_gt, atol=1e-5)) self.assertTrue(np.allclose(histogram_rt_cube, histogram_gt, atol=1e-5))