| 894 | } |
| 895 | |
| 896 | IntersectionResult< absl::InlinedVector< Point3D, 2 > > |
| 897 | plane_circle_intersection( const Plane& plane, const Circle& circle ) |
| 898 | { |
| 899 | const auto& circle_plane = circle.plane(); |
| 900 | const auto planes_intersection = |
| 901 | plane_plane_intersection( plane, circle_plane ); |
| 902 | if( !planes_intersection.has_intersection() ) |
| 903 | { |
| 904 | // The planes are parallel or nonintersecting. |
| 905 | return { planes_intersection.type }; |
| 906 | } |
| 907 | // The planes intersect in a line. Locate one or two points that |
| 908 | // are on the circle and line. If the line is t*D+P, the circle |
| 909 | // center is C, and the circle radius is r, then |
| 910 | // r^2 = |t*D+P-C|^2 = |D|^2*t^2 + 2*Dot(D,P-C)*t + |P-C|^2 |
| 911 | // This is a quadratic equation of the form |
| 912 | // a2*t^2 + 2*a1*t + a0 = 0. |
| 913 | const auto& line = planes_intersection.result.value(); |
| 914 | const Vector3D diff{ circle_plane.origin(), line.origin() }; |
| 915 | const auto a2 = line.direction().dot( line.direction() ); |
| 916 | const auto a1 = diff.dot( line.direction() ); |
| 917 | const auto a0 = diff.dot( diff ) - circle.radius() * circle.radius(); |
| 918 | |
| 919 | const auto discr = a1 * a1 - a0 * a2; |
| 920 | if( discr < 0. ) |
| 921 | { |
| 922 | // No real roots, the circle does not intersect the plane. |
| 923 | return { INTERSECTION_TYPE::none }; |
| 924 | } |
| 925 | absl::InlinedVector< Point3D, 2 > result; |
| 926 | CorrectnessInfo< absl::InlinedVector< Point3D, 2 > >::Correctness |
| 927 | first_correctness; |
| 928 | CorrectnessInfo< absl::InlinedVector< Point3D, 2 > >::Correctness |
| 929 | second_correctness; |
| 930 | const auto compute_correctness = [&first_correctness, |
| 931 | &second_correctness, &plane, |
| 932 | &circle]( |
| 933 | const Point3D& intersection ) { |
| 934 | auto plane_output = point_plane_distance( intersection, plane ); |
| 935 | first_correctness.first = |
| 936 | std::get< 0 >( plane_output ) <= GLOBAL_EPSILON; |
| 937 | first_correctness.second.emplace_back( |
| 938 | std::move( std::get< 1 >( plane_output ) ) ); |
| 939 | auto circle_output = point_circle_distance( intersection, circle ); |
| 940 | second_correctness.first = |
| 941 | std::get< 0 >( circle_output ) <= GLOBAL_EPSILON; |
| 942 | second_correctness.second.emplace_back( |
| 943 | std::move( std::get< 1 >( circle_output ) ) ); |
| 944 | }; |
| 945 | if( discr == 0. ) |
| 946 | { |
| 947 | // The quadratic polynomial has 1 real-valued repeated root. |
| 948 | // The circle just touches the plane. |
| 949 | compute_correctness( result.emplace_back( |
| 950 | line.origin() - line.direction() * ( a1 / a2 ) ) ); |
| 951 | } |
| 952 | else |
| 953 | { |