Circle intersection with lines

README.md

Circle intersection with (many) (oriented) lines with optimal runtime complexity

This code sample demonstrates how to compute the intersection with a circle of radius 1, and many oriented lines. Technically, oriented lines are equivalent to half-planes.

The function circle_line_intersection is doing all the work. The lines are specified as the 3 parameters of a 2d line equation, ax + by + c = 0. The function returns a list of pairs of points : the extremities of each arcs in counterclock wise order.

The algorithmic worst case complexity is O(n log(n)) with for n input segments, which I believe to be optimal. The algorithm behind the scene relies on Laguerre Voronoi diagrams on circles.

circle-lines-intersection.py

import numpy
from scipy.spatial import ConvexHull
import matplotlib.pyplot as plot
import matplotlib.patches as patches

	
	
# --- Sampling routines -------------------------------------------------------

'''
Pick N points uniformly from the unit disc
This sampling algorithm does not use rejection sampling.
'''
def disc_uniform_pick(N):
	angle = (2 * numpy.pi) * numpy.random.random(N)
	out = numpy.stack([numpy.cos(angle), numpy.sin(angle)], axis = 1)
	out *= numpy.sqrt(numpy.random.random(N))[:,None]
	return out



def enumerate_ring_consecutive_pairs(A):
	for item in zip(A[:-1], A[1:]):
		yield item
	yield (A[-1], A[0])
		
		
		
# --- Display -----------------------------------------------------------------

def arc_patch(theta1, theta2, ax=None, resolution=50, **kwargs):
	ax = plot.gca()  
	theta = numpy.linspace(numpy.radians(theta1), numpy.radians(theta2), resolution)
	points = numpy.vstack((numpy.cos(theta), numpy.sin(theta)))
	return patches.Polygon(points.T, closed=True, **kwargs)
	
	
	
def display(XY, theta, V, arc_list):
	fig = plot.figure()
	plot.axis('equal')
	plot.axis('off')	
	
	# Display input segments
	for XY_i, theta_i in zip(XY, theta):
		if theta_i != 0:
			alpha = numpy.arctan2(XY_i[1], XY_i[0])
			plot.gca().add_patch(arc_patch(numpy.degrees(alpha - .5 * theta_i), numpy.degrees(alpha + .5 * theta_i), zorder = 0, ec = '#f08080', fc = '#f0c0c0'))
		
			cos_theta_i, sin_theta_i = numpy.cos(.5 * theta_i), numpy.sin(.5 * theta_i)
			plot.plot([ cos_theta_i * XY_i[0] - sin_theta_i * XY_i[1],  cos_theta_i * XY_i[0] + sin_theta_i * XY_i[1]], 
			          [ sin_theta_i * XY_i[0] + cos_theta_i * XY_i[1], -sin_theta_i * XY_i[0] + cos_theta_i * XY_i[1]],
			          c = '#f08080')	
	# Display arcs
	for A, B in arc_list:
		plot.gca().add_patch(patches.Arc([0, 0], 2, 2, theta1 = numpy.degrees(numpy.arctan2(A[1], A[0])), theta2 = numpy.degrees(numpy.arctan2(B[1], B[0])), ls = '--', ec = '#8080f0', lw = 2))	

	arc_list = numpy.reshape(numpy.roll(numpy.reshape(arc_list, (2 * arc_list.shape[0], 2)), -1, axis = 0), (arc_list.shape[0], 2, 2))
		
	for A, B in arc_list:
		plot.gca().add_patch(patches.Arc([0, 0], 2, 2, theta1 = numpy.degrees(numpy.arctan2(A[1], A[0])), theta2 = numpy.degrees(numpy.arctan2(B[1], B[0])), ec = 'k', lw = 2, fc = '#f08080'))
	
	# Show all
	plot.show()
	
	
	
# --- Circular Laguerre-Voronoi diagram ---------------------------------------

def circle_line_intersection(line_eqns):
	half_sqrt_2 = .5 * 2 ** .5
	rot_pi_4 = half_sqrt_2 * numpy.array([[1, 1], [-1, 1]])
	
	# Unit equilateral triangle coordinates
	sqrt_3 = 3 ** .5
	equilateral_tri_coords = numpy.array([
		[  0.,  1.],
		[ -.5 * sqrt_3, -.5],
		[  .5 * sqrt_3, -.5]])	

	# Normalize the input
	XY_norm = numpy.sum(line_eqns[:,:2] ** 2, axis = 1) ** .5
	XY, Z = line_eqns[:,:2] / XY_norm[:,None], -line_eqns[:,2] / XY_norm
	
	# Break down circular sectors with angle > pi into pi / 2 sectors as
	# it allow to keep the algorithm simple yet optimal
	XYp, Zp = [], []
	for XY_i, Z_i in zip(XY, Z):
		if Z_i > half_sqrt_2:
			XYp.append(XY_i)
			Zp.append(Z_i)
		else:
			cos_theta_i, sin_theta_i = Z_i, (1. - Z_i ** 2) ** .5
			M = numpy.array([[cos_theta_i, sin_theta_i], [-sin_theta_i, cos_theta_i]])
			
			if Z_i <= half_sqrt_2 and Z_i > 0:
				M = numpy.dot(rot_pi_4, M.T)
				XYp.extend([numpy.dot(XY_i, M), numpy.dot(XY_i, M.T)])
				Zp.extend([half_sqrt_2, half_sqrt_2])
			elif Z_i <= 0 and Z_i > -half_sqrt_2:
				M = numpy.dot(rot_pi_4, M.T)				
				XYp.extend([numpy.dot(XY_i, M), numpy.dot(XY_i, M.T), XY_i])
				Zp.extend([half_sqrt_2, half_sqrt_2, half_sqrt_2])
			else:
				Ma = numpy.dot(rot_pi_4.T, M)
				XYp.extend([numpy.dot(XY_i, Ma), numpy.dot(XY_i, Ma.T)])
				Mb = numpy.dot(rot_pi_4, M)				
				XYp.extend([numpy.dot(-XY_i, Mb), numpy.dot(-XY_i, Mb.T)])			
				Zp.extend([half_sqrt_2, half_sqrt_2, half_sqrt_2, half_sqrt_2])
		
	XYp, Zp = numpy.array(XYp), numpy.array(Zp)
	
	# Augment the data with equilateral triangle coordinates, to ensure that
	# we have a convex hull even with an empty input
	XYp = numpy.concatenate([XYp, equilateral_tri_coords], axis = 0)
	Zp = numpy.concatenate([Zp, numpy.ones(3)], axis = 0)
	
	# Generate the Laguerre-Delaunay diagram
	delaunay_hull = ConvexHull(XYp / Zp[:,None]).vertices
	
	# Compute Laguerre-Voronoi vertices
	voronoi_vertices = numpy.empty((len(delaunay_hull), 2))
	for i, (j, k) in enumerate(enumerate_ring_consecutive_pairs(delaunay_hull)):
		voronoi_vertices[i] = numpy.linalg.solve(XYp[[j, k]], Zp[[j, k]])

	# Compute arcs
	voronoi_vertices_norm = numpy.sqrt(numpy.sum(voronoi_vertices ** 2, axis = 1))
	V = voronoi_vertices / voronoi_vertices_norm[:,None]

	arc_list = []
	for i, j in enumerate(delaunay_hull):
		cos_theta_j, sin_theta_j = Zp[j], (1 - Zp[j] ** 2) ** .5
		A = [cos_theta_j * XYp[j, 0] + sin_theta_j * XYp[j, 1], -sin_theta_j * XYp[j, 0] + cos_theta_j * XYp[j, 1]]
		B = [cos_theta_j * XYp[j, 0] - sin_theta_j * XYp[j, 1],  sin_theta_j * XYp[j, 0] + cos_theta_j * XYp[j, 1]]
				
		if voronoi_vertices_norm[i-1] >= 1:
			if voronoi_vertices_norm[i] >= 1:
				if numpy.cross(V[i] - cos_theta_j * XYp[j], XYp[j]) <= -sin_theta_j and numpy.cross(V[i-1] - cos_theta_j * XYp[j], XYp[j]) >= sin_theta_j:
					arc_list.append([(A, False), (B, False)])			
			else:
				arc_list.append([(A, False), (V[i], True)])
		else:
			if voronoi_vertices_norm[i] >= 1:		
				arc_list.append([(V[i-1], True), (B, False)])
			else:
				arc_list.append([(V[i-1], True), (V[i], True)])	

		if len(arc_list) >= 2:
			if arc_list[-2][1][1] and arc_list[-1][0][1]: 
				arc_list[-2][1] = arc_list[-1][1]
				arc_list.pop()

	if arc_list[-1][1][1] and arc_list[0][0][1]:
		arc_list[0][0] = arc_list[-1][0]
		arc_list.pop()

	arc_list = numpy.array([(A, B) for (A, a), (B, b) in arc_list if a == False and b == False and numpy.fabs(numpy.dot(A, B) - 1) > 1e-7])
	
	# Job done
	display(XY, 2 * numpy.arccos(Z), voronoi_vertices, arc_list)
	


# --- Main entry point --------------------------------------------------------

def get_random_dataset(sample_count, radius):
	XY = disc_uniform_pick(sample_count) 
	XY_norm = numpy.sqrt(numpy.sum(XY ** 2, axis = 1))
	XY /= XY_norm[:,None]
	theta = 2 * numpy.arccos(XY_norm) * radius
	return XY, theta
	
	
	
XY, theta = get_random_dataset(4, numpy.pi / 2)
circle_line_intersection(numpy.concatenate([XY, -numpy.cos(.5 * theta)[:,None]], axis = 1))