def outlierRemoval(time, flux):
	fluxDetrended = medianDetrend(flux, 3)
	out1 = plateau(fluxDetrended, 5 * np.std(fluxDetrended))
	out2 = plateau(-fluxDetrended, 5 * np.std(fluxDetrended))
	if out1 == [] and out2 == []:
		singleOutlierIndices = []
	else:
		outliers = np.append(out1, out2).reshape(-1,2)
		# Only want groups of one outlier, since > 1 may be transit points
		singleOutlierIndices = np.sort(outliers[(outliers[:,1] - outliers[:,0] == 1)][:,0])
	# Check periodicity of outliers, with PRECISION of 0.0205 days
	# 0.0205 days = 29.52 minutes = ~length of long cadence
	precision = 0.0205
	outlierTimes = time[singleOutlierIndices]
	diffs = [outlierTimes[i+1] - outlierTimes[i] for i in range(0, len(outlierTimes)-1)]
	diffs = [round(d, 5) for d in diffs]
	if len(singleOutlierIndices) >= 4:
		if len(set(diffs)) == len(diffs):
			possibleTimes = np.array([])
		else:
			period = max(set(diffs), key = diffs.count) # period = most common difference
			epoch = outlierTimes[diffs.index(period)]
			possibleTimes = np.arange(epoch, outlierTimes[-1] + 0.5*period, period)
		notOutliers = []
		for i in range(len(outlierTimes)):
			if np.any((abs(possibleTimes - outlierTimes[i]) < precision)):
				notOutliers.append(i)
		singleOutlierIndices = np.delete(singleOutlierIndices, notOutliers)
	elif len(singleOutlierIndices) == 3:
		if abs(diffs[0] - diffs[1]) < precision:
			singleOutlierIndices = []
	# Uncomment to see how the plotting algorithm worked for a lightcurve
	# ----------------------------- PLOTTING ----------------------------- #
	# plt.subplot(311)
	# plt.scatter(time, flux, marker = '.', s = 1, color = 'k', alpha = 1)
	# plt.scatter(time[singleOutlierIndices], flux[singleOutlierIndices], 
	#             s = 30, marker = 'o', facecolors = 'none', edgecolors = 'r')
	# plt.title('Original')
	# plt.subplot(312)
	# plt.scatter(time, fluxDetrended, marker = '.', s = 1, color = 'k', alpha = 1)
	# plt.scatter(time[singleOutlierIndices], fluxDetrended[singleOutlierIndices], 
	#             s = 30, marker = 'o', facecolors = 'none', edgecolors = 'r')
	# x1, x2, y1, y2 = plt.axis()
	# plt.hlines([-5*np.std(fluxDetrended), 5*np.std(fluxDetrended)], x1, x2, 
	#            color = 'b', linestyles = 'dashed')
	# plt.axis([x1, x2, y1, y2])
	# plt.title('Detrended')
	# plt.subplot(313)
	# plt.scatter(np.delete(time, singleOutlierIndices), np.delete(flux, singleOutlierIndices), 
	#             marker = '.', s = 1, color = 'k', alpha = 1)
	# plt.title('Outliers removed: ' + str(len(singleOutlierIndices)))
	# plt.show()
	# -------------------------------------------------------------------- #
	return np.delete(time, singleOutlierIndices), np.delete(flux, singleOutlierIndices)