"""
Classes and functions that facilitate handling geometry
"""

import numpy as np
import planar as pl
import matplotlib.pyplot as plt
import matplotlib.patches as mp
from matplotlib.collections import PatchCollection
from collections import OrderedDict

import workspace as ws
import algebra as alg
import tools as tls



def nearly_equal(a, b, tolerance):
	return np.abs(a - b) < tolerance

def dist(p1, p2):
	"""
	Distance between two points
	If one of the points is given as an array of points with the 
	correct shape, this function actually returns the result for the 
	array
	This shape is (2, npoints)
	"""
	return np.hypot(p2[1] - p1[1], p2[0] - p1[0])

def dist_cloud(x, y):
	""" 
	Mutual distances between a cloud of points 
	"""
	dists = np.hypot(alg.compl_mat(x), alg.compl_mat(y))
	pairs, dists = alg.one_side(dists)
	return pairs, dists

def pts_equal(p1, p2, tolerance):
	"""
	Check if p1 and p2 are the same point
	"""
	return nearly_equal(dist(p1, p2), 0, tolerance)

def get_midpoint(a, b):
	"""
	Get the point lying in the middle of points a and b
	"""
	mpx = 0.5 * (a[0] + b[0])
	mpy = 0.5 * (a[1] + b[1])
	return (mpx, mpy)

def points_in_circle(points, circle):
	"""
	Check which of a array of points is inside a circle
	The array of point needs to have the shape (2, npoints)
	IMPORTANT NOTE: Add a 1.e-10 % of the radius to make up for numerical 
	errors
	"""
	return dist(points, circle.c) <= circle.r * (1. + 1.e-10)

def get_slope(p1, p2):
	"""
	Get the slope of a line passing between p1 and p2
	"""
	try:
		slope = (p2[1] - p1[1]) / (p2[0] - p1[0])
	except ZeroDivisionError:
		slope = np.inf
	return slope

def get_orthoslope(slope):
	"""
	Get the slope of a line orthogonal to the line with a given slope 'slope'
	"""
	try:
		oslope = -1. / slope
	except ZeroDivisionError:
		oslope = np.inf
	return oslope

def angnegtopos(a):
	""" Given a negative angle 'a', return 2pi + a """
	twopi = 2. * np.pi
	n = abs(int(a / twopi)) + 1
	if a < 0:
		return n * twopi + a
	return a

def angle_segment(ref, p):
	"""
	Get angle that p has with respect to a reference point ref over the x-axis
	"""
	dx = p[0] - ref[0]
	dy = p[1] - ref[1]
	angle = np.angle(complex(dx, dy))
	# Turn it positive if it is negative
	angle = angnegtopos(angle)
	return angle

def angle_threepts(a, b, c):
	""" Get the angle at point b made by the segments ab and bc.
	Note: the order of a, b and c determines the value of the angle """
	# angle = angpostoneg(angle_segment(b, c)) - angpostoneg(angle_segment(b, a))
	angle = angle_segment(b, c) - angle_segment(b, a)
	# Turn it positive if it is negative
	return angnegtopos(angle)

def pts_angle(ref, p, numtolerance=1.e-9):
	"""
	Get angle that p has with respect to a reference point ref over the x-axis
	"""
	dx = p[0] - ref[0]
	dy = p[1] - ref[1]
	angle = np.angle(complex(dx, dy))
	if (angle < 0) & (np.abs(angle) > numtolerance):
		angle += 2. * np.pi
	return angle

def sort_by_angle(pts, ref=None):
	"""
	Sort a set of points pts by the angle they have with respect to ref
	"""
	# Reference point. If it is None, calculate it as a centroid of the points
	if ref is None:
		ref = pts.mean(axis=0).tolist()

	angles = []
	for p in pts:
		angles.append(pts_angle(ref, p))
	sorder = np.argsort(angles)

	return pts[sorder]

def intersection_two_lines(sl1, sl2):
	"""
	Calculate the intersection point between two straight lines in the 
	point-slope form
	"""
	ps1, ps2 = sl1.ps, sl2.ps
	p1, s1 = ps1
	p2, s2 = ps2
	r1 = sl1.equation
	r2 = sl2.equation
	# In case they are parallel
	if s1 == s2:
		return np.nan, np.nan
	# In case one of them has an infinite slope
	if np.abs(s1) == np.inf:
		x_int = p1[0]
		y_int = r2(x_int)
	elif np.abs(s2) == np.inf:
		x_int = p2[0]
		y_int = r1(x_int)
	else:
		x_int = (p2[1] - p1[1] + s1 * p1[0] - s2 * p2[0]) / (s1 - s2)
		y_int = r1(x_int)
	return x_int, y_int

def intersections_many_lines(lines_):
	""" Calculate all the intersections among one or various ensembles 
	of lines """

	# Get the vectors of number of lines, slopes and yintercpts from the ensembles
	nl = sum([lines.nl for lines in lines_])
	mvec = np.array([lines.m for lines in lines_]).flatten()
	nvec = np.array([lines.n for lines in lines_]).flatten()
	bxvec = np.array([lines.bxvec for lines in lines_]).flatten()
	byvec = np.array([lines.byvec for lines in lines_]).flatten()

	# Create the matrices to solve the problem
	m = alg.diagv(mvec)
	x = - alg.compl_mat(nvec) / alg.compl_mat(mvec)

	ivlines = np.where(np.abs(mvec) == np.inf)[0]
	nvlines = ivlines.size

	ihlines = np.where(np.abs(mvec) == 0.)[0]
	nhlines = ihlines.size

	# Modify values of x
	# Rows
	bx_rep_rows = np.repeat(bxvec[ivlines], nl).reshape((nvlines, nl))
	x[ivlines] = bx_rep_rows[:]

	# Columns
	bx_rep_cols = np.tile(bxvec[ivlines], nl).reshape((nl, nvlines))
	x[:, ivlines] = bx_rep_cols[:]

	# Infinities to nans
	x[np.where(np.abs(x) == np.inf)] = np.nan
	x = alg.make_sym_rem_nans(x)

	# Wherever x = nan, place any dummy number (it won't affect the 
	# theoretical result). This will avoid y being all nan
	xx = x.copy()
	xx[np.where(np.isnan(x))] = 0.

	# Compute y
	y = np.dot(m, xx) + alg.vectomat(nvec)

	# Modify values of y
	by_rep_rows = np.repeat(byvec[ihlines], nl).reshape((nhlines, nl))
	y[ihlines] = by_rep_rows[:]

	# Columns
	by_rep_cols = np.tile(byvec[ihlines], nl).reshape((nl, nhlines))
	y[:, ihlines] = by_rep_cols[:]

	# Infinities to nans
	y[np.where(np.abs(y) == np.inf)] = np.nan
	y = alg.make_sym_rem_nans(y)

	return x, y

def intersection_segments(a, b):
	""" 
	Intersection point between two segments a and b 
	"""

	x_int, y_int = intersection_two_lines(a.line, b.line)

	# Now check if this point belongs to the segments

	if a.contains_point((x_int, y_int), include_ends=False) \
		& b.contains_point((x_int, y_int), include_ends=False):
		return x_int, y_int
	else:
		return np.nan, np.nan

def crossing_segments(a, b):
	""" 
	Determine if two segments a and b are crossing (their intersection 
	point is not on any of their limits
	"""
	x_int, y_int = intersection_segments(a, b)
	if np.isnan(np.array([x_int, y_int])).any():
		return False, (x_int, y_int)

	a0x, a0y = a.a
	a1x, a1y = a.b
	b0x, b0y = b.a
	b1x, b1y = b.b
	x = np.array([a0x, a1x, b0x, b1x])
	y = np.array([a0y, a1y, b0y, b1y])
	_, dists = dist_cloud(x - x_int, y - y_int)
	if (dists <= 1.e-9).any():
		return False, (x_int, y_int)
	return True, (x_int, y_int)

def overlapping_segments(a, b):
	"""
	Determine if two segments overlap. For this: 
	1) check that their lines are the same
	2) check that a point of one of them (any) falls onto the other 
	segment (eg. a point in a contained by b or viceversa)
	"""

	# Get their value of equation(x=0.)
	y0a = a.line.equation(0)
	y0b = b.line.equation(0)

	tolerance = 1.e-9
	same_y0 = nearly_equal(y0a, y0b, tolerance=tolerance)
	same_slope = nearly_equal(a.slope, b.slope, tolerance=tolerance)
	same_lines = same_y0 and same_slope

	# If the lines are not the same, bye bye
	if not same_lines:
		return False
	
	# Check if any point from a falls on b or viceversa
	for p in a.points:
		if b.contains_point(p, include_ends=False):
			return True
	for p in b.points:
		if a.contains_point(p, include_ends=False):
			return True

	# If nothing happened so far, the segments don't overlap
	return False

def intersection_polygons(a, b):
	""" Return the polygon, if it exists, resulting from the 
	intersection between two polygons a and b """
	
	# List the points that are inside one another
	inter_points = []
	for pa in a.verts:
		if b.contains_point(pa):
			inter_points.append(pa)
	for pb in b.verts:
		if a.contains_point(pb):
			inter_points.append(pb)

	# Find the intersections between the segments

	for sa in a.segments.values():
		for sb in b.segments.values():
			xi, yi = intersection_segments(sa, sb)
			if (not np.isnan(xi)) and (not np.isnan(yi)):
				inter_points.append(np.array([xi, yi]))

	# Remove repeated points
	avoid_these_points = []
	for i, p1 in enumerate(inter_points):
		for j, p2 in enumerate(inter_points[i+1:]):
			j += i + 1
			# print(i, j)
			if pts_equal(p1, p2, 1.e-9):
				avoid_these_points.append(j)

	# Remove redundant points
	inter_points = np.array([p for (i, p) in enumerate(inter_points) if not i in avoid_these_points])

	# Finally, create the polygon/segment/point and return it if it has any points
	if len(inter_points) > 1:
		return Polygon(sort_by_angle(inter_points))
	else:
		return None

###############################################################################


class Point():
	"""
	Point.
	NOTE: The majority of the code here is written to treat points as just 
	pairs of values, not to use this class. But this class is defined here 
	to be used in some circumstances and also to use it in a future version 
	of the code.
	"""

	def __init__(self, x, y):
		self.x = x
		self.y = y
		self.array = np.array([x, y])

	def draw(self, ax, colour='k', marker='o', markersize=1, alpha=1.):
		ax.plot(
				self.x, 
				self.y, 
				color=colour, 
				marker=marker, 
				markersize=markersize, 
				alpha=alpha, 
			)

class StraightLine():
	"""
	Straight Line
	"""

	def __init__(self, point, slope):
		self.point = point
		self.slope = slope
		self.ps = (point, slope)
		self.set_equation()
		self.vertical = np.abs(self.slope) == np.inf
		self.horizontal = nearly_equal(np.abs(self.slope), 0., tolerance=1.e-9)

	def set_equation(self):
		"""
		Equation for the straight line
		"""
		p, s = self.point, self.slope
		self.y_ord = p[1] - s * p[0]
		y = self.y_ord
		self.equation = lambda x: y + s * x

	def draw(self, ax, colour='k', linewidth=1, alpha=1.):
		"""
		Draw the line on ax
		"""
		xleft, xright = ax.get_xlim()
		yleft, yright = self.equation(xleft), self.equation(xright)
		if (self.slope == np.inf) or (self.slope == -np.inf) or (np.isnan(self.slope)):
			xleft = xright = self.ps[0][0]
			yleft, yright = ax.get_ylim()
		ax.plot(
				[xleft, xright], 
				[yleft, yright], 
				color=colour, 
				lw=linewidth, 
				alpha=alpha, 
			)

	def contains_point(self, p):
		"""
		Check if a point p is on this line
		"""
		xp, yp = [tls.list_to_number(x) for x in p]
		if self.vertical:
			return nearly_equal(xp, self.point[0], tolerance=1.e-9)
		return nearly_equal(yp, self.equation(xp), tolerance=1.e-9)


class Segment():
	"""
	Segment
	"""

	def __init__(self, points):

		# Validation
		if len(points) != 2:
			msg = 'WARNING: Defined a segment with %i points:\n%s\nSegment not created. Returned None'%(len(points), str(points))
			try:
				ws.log(msg)
			except AttributeError:
				print(msg)
			return None

		# Clean the points. Give them a nice and clean format
		points = (
			tuple([tls.list_to_number(x) for x in points[0]]), 
			tuple([tls.list_to_number(x) for x in points[1]])
		)
		self.points_list = list(points)
		# Add attributes
		self.update(points)

	def __str__(self):
		""" Print the segment with its values """
		return 'Segment. a: (%f, %f); b: (%f, %f)'%tuple([x for x in self.points.flatten()])

	def update(self, points):
		""" Update the values of the segment """
		self.points = np.array(list(points))
		self.a = points[0]
		self.b = points[1]
		# Get the straight line onto which the segment is
		self.slope = get_slope(*self.points.tolist())
		self.line = StraightLine(self.a, self.slope)
		self.vertical = self.line.vertical
		self.horizontal = self.line.horizontal
		# Length
		self.length = dist(self.a, self.b)

	def change_point(self, index, point):
		""" Change one point of the segment; hence change it """
		# Clean the point
		point = tuple([tls.list_to_number(x) for x in point])
		# Add it to the list
		self.points_list[index] = point
		# Update the segment
		self.update(self.points_list)

	def contains_point(self, p, include_ends=True):
		""" Check if a point is on the segment """

		# Format point correctly
		p = tuple([tls.list_to_number(x) for x in p])

		# Now check if the line contains the point
		# If not, return False
		if not self.line.contains_point(p):
			return False

		# Now check that the point falls in the rectangle where seg is the 
		# diagonal
		xp, yp = p
		a, b = self.a, self.b

		# Note: if the segment is vertical, include_ends = True
		include_ends = include_ends or self.vertical or self.horizontal

		if include_ends:
			inx = (a[0] <= xp <= b[0]) | (a[0] >= xp >= b[0])
			iny = (a[1] <= yp <= b[1]) | (a[1] >= yp >= b[1])
		else:
			inx = (a[0] < xp < b[0]) | (a[0] > xp > b[0])
			iny = (a[1] < yp < b[1]) | (a[1] > yp > b[1])

		if inx and iny:
			return True
		else:
			return False

	def draw(self, ax, colour='k', linewidth=1, alpha=1.):
		""" Draw the segment on ax """
		ax.plot(
				(self.a[0], self.b[0]), 
				(self.a[1], self.b[1]), 
				c=colour, 
				lw=linewidth, 
				alpha=alpha, 
			)


class Circle():
	"""
	Circle
	"""

	def __init__(self, c, r):
		self.c = c
		self.r = r

	def draw(self, ax, colour='k', linewidth=1, alpha=1.):
		""" 
		Draw the circle 
		"""
		plt_circle = plt.Circle(self.c, self.r, color=colour, 
			linewidth=linewidth, alpha=alpha, fill=False)
		ax.add_artist(plt_circle)

		# If ax limits are smaller than the circle, expand them
		self.xmin = self.c[0] - self.r
		self.xmax = self.c[0] + self.r
		self.ymin = self.c[1] - self.r
		self.ymax = self.c[1] + self.r
		ax_xmin, ax_xmax = ax.get_xlim()
		ax_ymin, ax_ymax = ax.get_ylim()
		margin = 0.05 * max([self.xmax - self.xmin, self.ymax - self.ymin])
		# x-axis
		if self.xmin < ax_xmin:
			ax.set_xlim(left=self.xmin - margin)
		if self.xmax > ax_xmax:
			ax.set_xlim(right=self.xmax + margin)
		# y-axis
		if self.ymin < ax_ymin:
			ax.set_ylim(bottom=self.ymin - margin)
		if self.ymax > ax_ymax:
			ax.set_ylim(top=self.ymax + margin)


class Polygon(object):
	"""
	Class for any polygon that is defined with a list of verts
	*args and **kwargs account for any arguments and keyword arguments that may enter in instances of childs of this class
	"""

	def __new__(cls, verts, *args, **kwargs):

		lv = len(verts)

		# Validation
		if lv == 0:
			msg = 'ERROR: Point list or array is empty.'
			try:
				ws.log(msg)
			except AttributeError:
				print(msg)
			else:
				ws.terminate()

		elif lv == 1:
			msg = 'WARNING: Defined polygon with only 1 point:\n%s. Returing point'%str(verts)
			try:
				ws.log(msg)
			except AttributeError:
				print(msg)
			return Point(*verts[0])

		elif lv == 2:
			msg = 'WARNING: Defined polygon with only 2 points:\n%s. Returing Segment'%str(verts)
			try:
				ws.log(msg)
			except AttributeError:
				print(msg)
			return Segment(verts)

		# elif lv == 3:
			# This is a triangle
			# return super(Triangle, cls).__new__(cls)
			# return super(Triangle, cls).__new__(cls)
			# I didn't get this to work...
			pass

		elif lv >= 3:
			# Now, this is a proper polygon. So, proceed.
			return super(Polygon, cls).__new__(cls)
		
	def __init__(self, verts, *args):
		# print('2: Creating polygon')

		# Attributes
		self.verts = verts
		# If it's not an array, make it be
		try:
			verts.shape
		except AttributeError:
			self.verts = np.array(verts)
		self.nverts = len(verts)
		self.x, self.y = self.verts.T
		# Compute attributes
		self.set_segments()
		self.set_angles()
		self.get_area()
		# self.get_centroid()
		self.get_circumcircle()
		# planar.Polygon. I need it for planar.Polygon.contains_point
		self.plpol = pl.Polygon(self.verts.tolist())

	def set_segments(self):
		"""
		Set the segments of this polygon. Also their lengths and slopes
		"""

		segments = OrderedDict()
		sides = []
		slopes = []

		# Iterate over self.nverts - 2 points
		for i in range(self.nverts - 1):
			seg = Segment(self.verts[i:i+2])
			sides.append(seg.length)
			slopes.append(seg.slope)
			segments[(i, i+1)] = seg

		# Last segment
		seg = Segment((self.verts[i+1], self.verts[0]))
		sides.append(seg.length)
		slopes.append(seg.slope)
		segments[(i+1, 0)] = seg

		# Add as attribute
		self.segments = segments
		self.sides = np.array(sides)
		self.slopes = np.array(slopes)

	def set_angles(self):
		"""
		Get the angles for each one of the vertices
		Do it COUNTER-CLOCKWISE, which is how the points are oriented
		"""
		self.angles = np.zeros(self.nverts)
		for i in np.arange(0, self.nverts - 1, 1):
			c, b, a = self.verts[i - 1], self.verts[i], self.verts[i + 1]
			self.angles[i] = angle_threepts(a, b, c)
		# Last vertex
		c, b, a = self.verts[-2], self.verts[-1], self.verts[0]
		self.angles[-1] = angle_threepts(a, b, c)

		# Extra:
		# Tell if the polygon is convex
		self.is_convex = (self.angles <= np.pi).all()

		# Angles to degrees
		self.angles_deg = (180. / np.pi) * self.angles

	def get_area(self):
		"""
		Calculate the area of this polygon
		"""
		x, y = self.x, self.y
		self.signed_area = 0.5 * (np.dot(x, np.roll(y, 1)) \
			- np.dot(y, np.roll(x, 1)))
		# NOTE: This only works for polygons with no crossing segments, right?
		# WARNING: True, this is correct only for polygons with non-intersecting segments
		self.area = np.abs(self.signed_area)

	def get_centroid(self):
		"""
		Get the centroid of the polygon
		ERROR: SOMETHING IS WRONG WITH THIS
		"""
		oosixa = 1. / (6. * self.signed_area)
		x, y = self.x, self.y
		a = x[:-1] * y[1:] - x[1:] * y[:-1]
		cx = oosixa * ((x[:-1] + x[1:]) * a).sum()
		cy = oosixa * ((y[:-1] + y[1:]) * a).sum()
		self.centroid = (cx, cy)

	def get_circumcircle(self):
		""" 
		Create the circumcircle with its properties 
		"""

		# First and foremost, we have to try with the two most distant 
		# points. They will give us THE SMALLEST POSSIBLE CIRCLE
		# If this circle does not contain all the points (it does not
		# meet the condition), then we use the algorithm of the 
		# triangles (which always works)

		# Find the two most distant points
		xv, yv = self.verts.T
		pairs, distances = dist_cloud(xv, yv)
		self.most_distant_points = self.verts[pairs[np.where(distances == distances.max())][0]]
		# Get the centre
		a, b = self.most_distant_points
		circumcentre = get_midpoint(a, b)
		circumradius = dist(circumcentre, a)
		circumcircle = Circle(circumcentre, circumradius)

		# Provisionally save the previous circumcircle for comparison
		self.circumcircle_prov = circumcircle

		# Check if it meets the condition
			# self.circumcircle = circumcircle
		if not (points_in_circle(self.verts.T, circumcircle)).all():
		# else:
			# It didn't meet the condition, so go with the triangles

			# First off, get all the possible triangles
			# And get their circumcircles
			self.triangles = {}
			possible_circumcircles = []
			vertices = range(self.nverts)
			# Loop over all possible triangles
			for i in vertices:
				for j in vertices[i+1:]:
					for k in vertices[j+1:]:
						tverts = (
								self.verts[i], 
								self.verts[j], 
								self.verts[k]
							)

						# Create triangle
						# print('tverts:', tverts)
						triangle = Triangle(tverts)

						# Check if the circumcircle contains all points
						if (points_in_circle(
								self.verts.T, triangle.circumcircle)
							).all():
							pc = triangle.circumcircle

							if len(possible_circumcircles) == 0:
								# If this circumcircle is the smallest one so far, keep it
								circumcircle = pc
								circumcentre = pc.c
								circumradius = pc.r
							elif pc.r < circumcircle.r:
								# If we have more possibilities, choose it if it's the smallest one
								circumcircle = pc
								circumcentre = pc.c
								circumradius = pc.r

							# Even if this was positive, add it to the pc's 
							# that meet the condition
							possible_circumcircles.append(pc)
						# Save triangle
						self.triangles[(i, j, k)] = triangle
			# Save the possible circumcircles
			self.possible_circumcircles = possible_circumcircles

		# Save attributes
		self.circumcircle = circumcircle
		self.circumcentre = circumcentre
		self.circumradius = circumradius
		self.circumdiameter = 2. * circumradius

		# print('Polygon.')
		# print('Centroid:', self.centroid)
		# print('Circumcentre:', self.circumcentre)
		
	def contains_point(self, p, include_contour=True):
		"""
		Check if the polygon contains a point p
		"""

		# First, check if its planar.Polygon has it
		if self.plpol.contains_point(p):
			return True

		# If it doesn't, check if it's on the contour
		if include_contour:
			for seg in self.segments.values():
				if seg.contains_point(p, include_ends=True):
					return True

		# If we're here, the point is not in or on this polygon
		return False

	def draw(self, ax, colour='k', alpha=1., linestyle='-', linewidth=1., 
		markers=None, facecolor='none', dashes=(1, 1), zorder=0):
		"""
		Draw the polygon on ax
		"""
		# Add the last point
		ptstoplot = self.verts.tolist()
		ptstoplot.append(ptstoplot[0])
		ptstoplot = np.array(ptstoplot)
		if facecolor == 'none':
			try:
				if linestyle == '--':
					ax.plot(ptstoplot[:, 0], ptstoplot[:, 1], c=colour, alpha=alpha,
						ls=linestyle, dashes=dashes, lw=linewidth, marker=markers, zorder=zorder)
				else:
					ax.plot(ptstoplot[:, 0], ptstoplot[:, 1], c=colour, alpha=alpha,
						ls=linestyle, lw=linewidth, marker=markers, zorder=zorder)
			except IndexError:
				pass
		else:
			patches = []
			filled_pol = mp.Polygon(self.verts)
			patches.append(filled_pol)
			collection = PatchCollection(patches, facecolors=facecolor)
			ax.add_collection(collection)
			collection.set_alpha(alpha)
			collection.set_edgecolor(colour)
			collection.set_linewidth(linewidth)
			collection.set_linestyle(linestyle)


class Triangle(Polygon):
	"""
	Class for a Triangle that contains all the important parts
	I need
	"""
	# def __init__(self, verts, *args, indices=(0, 1, 2), fas=None):
	def __init__(self, verts, indices=(0, 1, 2), fas=None):
		# print('1: Creating triangle')
		self.indices = indices
		self.verts = np.array([verts[i] for i in indices])
		Polygon.__init__(self, self.verts)
		# Get relevant attributes
		# print('Here we have a triangle with %i vertices'%self.nverts)
		self.stablish_pairs()
		self.stablish_lines()
		# self.get_sides()
		self.get_incentre()
		self.get_medians()
		self.get_altitudes()
		self.get_circumcircle()
		self.get_orthodistances()
		# Check if it is inside a fascicle
		# Depending on that, the pseudo_incentre will be calculated or not
		if fas is not None:
			self.fas = fas
			self.fibers = fas.fibers
			self.get_pseudo_incentre()

	def stablish_pairs(self):
		self.pairs = []
		self.pairs.append((self.indices[1], self.indices[2]))
		self.pairs.append((self.indices[2], self.indices[0]))
		self.pairs.append((self.indices[0], self.indices[1]))

	def stablish_lines(self):
		# NOTE: Lines are actually 'segments', to be more precise
		# Lines
		# The index should be the same as the index of the opposite 
		# vertex
		self.lines = []
		self.lines.append([self.verts[1], self.verts[2]])
		self.lines.append([self.verts[2], self.verts[0]])
		self.lines.append([self.verts[0], self.verts[1]])

	# def get_sides(self):
	# 	# Line lengths or sides and slopes
	# 	sides = []
	# 	slopes = []
	# 	for line in self.lines:
	# 		sides.append(dist(line[0], line[1]))
	# 		# Slopes for the lines
	# 		slopes.append(get_slope(line[0], line[1]))
	# 	self.sides = np.array(sides)
	# 	self.slopes = np.array(slopes)

	def get_incentre(self):
		# Position of the incentre
		self.incentre = np.array([self.sides[i] \
		 * np.array(self.verts[i]) for i, p in \
		  enumerate(self.verts)]).sum(axis=0) / self.sides.sum()		
		# print(self.sides)
		a, b, c = self.sides
		s = 0.5 * self.sides.sum()
		self.inradius = np.sqrt((s - a) * (s - b) * (s - c) / s)		

	def get_pseudo_incentre(self):
		"""
		If the incentre is inside a cell, move it outwards
		"""
		# First, the pseudo_incentre is a copy of incentre
		self.pseudo_incentre = list(self.incentre)
		# Check if incentre is inside any of the triangle's cells
		for i in self.indices:
			cell = self.fibers[i]
			circ = cell.circumph
			c, r = circ
			# circ = sg.Circle(tup2sympoint(c), r)
			if inside_circle(self.incentre, circ):
				# Move pseudo_incentre
				# Place it on the cell's boundary, on the line 
				# perpendicular to the triangle's wall opposite to its
				# center. So first, find the equation of that line

				# Get the indices of the cells pairing on the side 
				# I want
				othercells = [j for j in self.indices if j != i]
				slope = get_slope(self.fibers[othercells[0]].center, 
					self.fibers[othercells[1]].center)
				orthoslope = get_orthoslope(slope)
				ortholine = StraightLine(c, orthoslope)
				# ortholine = sg.Line(cell.center, slope=orthoslope)
				inters = intersection_sl_circ(ortholine, circ)
				# Of the two points from inters, I need the closest 
				# one to the incentre
				self.pseudo_incentre = find_closest(inters, self.incentre)

	def get_medians(self):
		""" Calculate the medians of the triangle """
		a, b, c = self.sides
		self.medians = 0.5 * np.array([
			np.sqrt(2. * (b*b + c*c) - a*a),
			np.sqrt(2. * (a*a + c*c) - b*b), 
			np.sqrt(2. * (b*b + a*a) - c*c)
			])
		xa, xb, xc = self.verts

	def get_altitudes(self):
		""" Get altitudes from the sides """
		a, b, c = self.sides
		s = 0.5 * self.sides.sum()
		self.altitudes = \
			2. * np.sqrt(s * (s - a) * (s - b) * (s - c)) / self.sides

	def get_circumcircle(self):
		""" Get and save the circumcircle of this triangle """
		midpoints = []
		side_bisectors = []
		for slope, segment in zip(self.slopes, self.segments.values()):
			a, b = segment.a, segment.b
			mp = get_midpoint(a, b)
			midpoints.append(mp)
			side_bisectors.append(StraightLine(mp, get_orthoslope(slope)))
		
		# Find the circumcentre, which is the intersection between two
		# of these lines
		circumcentre = intersection_two_lines(
			side_bisectors[0], 
			side_bisectors[1]
			)
		
		# Circumradius
		circumradius = dist(circumcentre, self.verts[0])

		# Get circumcircle
		self.circumcircle = Circle(circumcentre, circumradius)

		# Save attributes
		self.circumcentre = circumcentre
		self.circumradius = circumradius
		self.midpoints = midpoints
		self.side_bisectors = side_bisectors

	def get_orthodistances(self):
		""" Get distances from the sides to the orthocenter """
		self.orthodistances = np.empty(3)
		for i in range(3):
			# sign = np.sign(np.cos(self.angles[i]))
			sign = -np.sign(np.cos(self.angles[i]))
			# vh = sign * np.sqrt(4 * self.circumradius ** 2 - self.sides[i] ** 2)
			vh = np.sqrt(4 * self.circumradius ** 2 - self.sides[i] ** 2)
			self.orthodistances[i] = sign * (self.altitudes[i] - vh)