There are $n$ points in a 2-D plane and each is given by its $x$ and $y$ coordinates. They are stored in an array in an ascending order with respect to $x$.
All points are connected together by line segments ($n\choose 2$ lines). Within the complete graph created by these line segments, how can we find a line which has the most number of intersections?
In the given example, $AD$ (marked green) is the line with the most intersections (4) in a set 6 points.
You can use the Bentley–Ottmann algorithm for this. Given a set of $n$ line segments with $k$ intersections, the algorithm can identify all intersections in $O((n+k)\log n)$ time and $O(n)$ space. In cases where $k=o(\frac{n^2}{\log n})$ (that is, cases in which $k$ has an appropriate upper bound) this offers time savings versus a naive $O(n^2)$ algorithm that simply compares all segments.
More generally, this problem can be approached using any of a number of sweep line algorithms.
The trick, then, is to increment a segment's value in a scoring hash table each time it is involved in an intersection.
A couple of hours of fiddling with CGAL didn't reveal an obvious way to do this and other implementations produced incorrect answers due to floating-point issues at the lines' end points. Nonetheless, this is the most computationally-efficient way to approach the problem.
I've copied my stab at building a CGAL implementation below with the appropriate spot to edit the code noted:
//Compile with: g++ -g 23222-line-with-most-intersections.cpp -lCGAL -lgmp -lmpfr
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Arr_segment_traits_2.h>
#include <CGAL/Surface_sweep_2.h>
#include <CGAL/Surface_sweep_2_algorithms.h>
//#include <CGAL/Sweep_line_2.h>
#include <CGAL/Surface_sweep_2/Default_visitor.h>
#include <CGAL/Surface_sweep_2/Surface_sweep_2_utils.h>
#include <list>
#include <vector>
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef Kernel::Point_2 Point_2;
typedef CGAL::Arr_segment_traits_2<Kernel> Traits_2;
typedef Traits_2::Curve_2 Segment_2;
namespace CGAL {
namespace Surface_sweep_2 {
template <typename GeometryTraits_2, typename OutputIterator,
typename Allocator_ = CGAL_ALLOCATOR(int)>
class IntersectionCounter :
public Default_visitor<IntersectionCounter<GeometryTraits_2,
OutputIterator,
Allocator_>,
GeometryTraits_2, Allocator_>
{
public:
typedef GeometryTraits_2 Geometry_traits_2;
typedef OutputIterator Output_iterator;
typedef Allocator_ Allocator;
private:
typedef Geometry_traits_2 Gt2;
typedef IntersectionCounter<Gt2, Output_iterator, Allocator>
Self;
typedef Default_visitor<Self, Gt2, Allocator> Base;
public:
typedef typename Base::Event Event;
typedef typename Base::Subcurve Subcurve;
typedef typename Subcurve::Status_line_iterator Status_line_iterator;
typedef typename Gt2::X_monotone_curve_2 X_monotone_curve_2;
typedef typename Gt2::Point_2 Point_2;
typedef typename Base::Surface_sweep_2 Surface_sweep_2;
protected:
Output_iterator m_out; // The output points.
public:
IntersectionCounter(Output_iterator out) :
m_out(out)
{}
template <typename CurveIterator>
void sweep(CurveIterator begin, CurveIterator end)
{
std::vector<X_monotone_curve_2> curves_vec;
std::vector<Point_2> points_vec;
curves_vec.reserve(std::distance(begin,end));
make_x_monotone(begin, end,
std::back_inserter(curves_vec),
std::back_inserter(points_vec),
this->traits());
//Original curves get converted into x-monotone curves here, but, since they
//are segments, their ordering and data appears to be unaltered
std::cout<<"x-monotone curves\n";
for(auto &x: curves_vec)
std::cout<<x<<" "<<(&x)<<std::endl;
std::cout<<"x-monotone points\n";
for(auto &x: points_vec)
std::cout<<x<<std::endl;
//Perform the sweep
Surface_sweep_2* sl = this->surface_sweep();
sl->sweep(curves_vec.begin(), curves_vec.end(),
points_vec.begin(), points_vec.end());
}
bool after_handle_event(Event* event,
Status_line_iterator /* iter */,
bool /* flag */)
{
//TODO: Magic should happen here
if ((
event->is_intersection() ||
event->is_weak_intersection()) && event->is_closed())
{
*m_out = event->point();
++m_out;
}
return true;
}
Output_iterator output_iterator() { return m_out; }
};
} // namespace Surface_sweep_2
namespace Ss2 = Surface_sweep_2;
template <typename CurveInputIterator, typename OutputIterator, typename Traits>
OutputIterator CountIntersections(
CurveInputIterator curves_begin,
CurveInputIterator curves_end,
OutputIterator points,
Traits &tr
){
// Define the surface-sweep types:
typedef Ss2::IntersectionCounter<Traits, OutputIterator> Visitor;
typedef Ss2::Surface_sweep_2<Visitor> Surface_sweep;
// Perform the sweep and obtain the intersection points.
Visitor visitor(points);
Surface_sweep surface_sweep(&tr, &visitor);
visitor.sweep(curves_begin, curves_end);
return visitor.output_iterator();
}
template <typename CurveInputIterator, typename OutputIterator>
OutputIterator CountIntersections(
CurveInputIterator curves_begin,
CurveInputIterator curves_end,
OutputIterator points
){
typedef typename std::iterator_traits<CurveInputIterator>::value_type Curve;
typename Default_arr_traits<Curve>::Traits traits;
return CountIntersections(curves_begin, curves_end, points, traits);
}
} // namespace CGAL
int main(){
//Points as extracted from https://scicomp.stackexchange.com/q/23222/17088
const std::vector<Point_2> pts = {
Point_2( 57,931),
Point_2(447,699),
Point_2(899,748),
Point_2(863,137),
Point_2(530, 67),
Point_2(142,282)
};
//Points are fully connected
std::vector<Segment_2> segments;
for(int i=0; i<pts.size();i++)
for(int j=i+1;j<pts.size();j++){
segments.emplace_back(pts[i],pts[j]);
std::cout<<pts[i]<<"\n"<<pts[j]<<"\n\n";
}
// Compute all intersection points.
std::list<Point_2> ipts;
CGAL::CountIntersections(segments.begin(), segments.end(), std::back_inserter(ipts));
for(const auto &x: segments)
std::cout<<(&x)<<std::endl;
// Print the result.
std::cout << "Found " << ipts.size() << " intersection points: " << std::endl;
std::copy(ipts.begin(), ipts.end(),
std::ostream_iterator<Point_2>(std::cout, "\n"));
return 0;
}
The answer to this depends on the hypotheses you consider for the points. If one or more points are collinear or if more than two segments meet at a point then things may become more complicated.
So I'll assume the simple configuration where no three points are collinear and no three segments meet. Then the number of points a segment intersects is $L \times R$ where $L$ denotes the number of points on the "left" and $R$ on the "right" of the given segment. Since $R+L+2$ should give the total number of points, which is fixed, the maximal value of $L\times R$ is attained when $R$ and $L$ are as close as possible. A double loop on the vertices should give you the right answer.
(I'm sure other more optimal algorithms could exist)
My heuristic method is that the line with the most intersections will be a line joining points on the convex hull. You could use a variant of the rotating calipers algorithm and test the lines between antipodal points on the hull. The number of intersections depends on the points "between" the antipodal pair. The exact number depends on how many points are on each side of the line joining the antipodal pair. The number of points on one side multiplied by the number of points on the other determines the maximum number of intersections possible. The exact number may be less because multiple lines may intersect in the same place. Thus, actual evaluations will be needed but this basic heuristic can narrow down the number of lines to check.
In the given example, line AD are basically antipodal points. Since points B and C are on one side and points E and F are on the other, then there are 2 x 2 = 4 lines between them and thus 4 intersections maximum (which agrees with the diagram). If B was on the other side of AD, then there would be 3 and 1 points on either side, respectively and BF and BE wouldn't intersect AD, but BC would leaving 3 intersections. The intersections Q and R are close to each other, and even if the points are not collinear Q and R could be coincident resulting in fewer intersections than the maximum given by this heuristic. If 3 or more points are collinear, then there will also be fewer intersections than the maximum.
You simply take each $n\choose 2$ line segment and then for each other line segment see if they intersect. This is tallied for every single line segment. The line segment with the most intersection is your segment.
^ The above part is easy... It's like this is somebody's idea of a joke. Anyways...
Now the method to see if there is an intersection has occurred is the more technical part. I'm not well versed in Computer Graphics, so I wouldn't know conclusively right away the best method.
