# How can I find a line segment with the most intersections along with the coordinates of the intersection points?

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.

• @ChristianClason now? I'm sorry if this sounding rude, I'm not good in grammar. – Anuraj Kathait Feb 26 '16 at 19:26
• Yes, much better (and don't worry)! – Christian Clason Feb 26 '16 at 19:29
• The line you search for is an existent line which connects two points or it could be a different line? Can you provide more details regarding the constraints on the line you search for? – Aurelian Tutuianu Mar 2 '16 at 11:44

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)

• I don't understand your method. can you provide more details? a double loop on all vertices means you iterate over all pairs of vertices. I can't see how this will do it. if you check for each vertice if it is left or right of the line defined by a pair then you have another loop. consider the line defined by the pair B and E. A and C are on opposite sites of this line but the line segments do not intersect. – miracle173 May 30 '18 at 6:02
• You are right... With what I said you can only count the lines which intersect... If you speak about segments then not all pairs of points on opposite sides will give intersecting segments... – Beni Bogosel May 30 '18 at 14:21

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.

• put an arbitrary set of points in the interior of a triangle.you get three line segment if you join the three points of the convex hull. none of them intersects another line segment – miracle173 May 30 '18 at 5:42

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.

• so if n=10000 you have to check 2.5*10^15 pairs of line segments and if n=10^6 you have to check 2.5*10^23 pairs of line segments. That will take a while. – miracle173 May 30 '18 at 6:16