386 lines
11 KiB
Python
Executable File
386 lines
11 KiB
Python
Executable File
"""
|
|
Library functions for graham algorithm implementation (find the convex hull of a given list of points).
|
|
"""
|
|
|
|
from copy import deepcopy
|
|
|
|
import numpy as np
|
|
|
|
|
|
def clean_ROI(image):
|
|
"""
|
|
Remove instruments borders from an observation.
|
|
"""
|
|
H, J = [], []
|
|
|
|
shape = np.array(image.shape)
|
|
row, col = np.indices(shape)
|
|
|
|
for i in range(0, shape[0]):
|
|
r = row[i, :][image[i, :] > 0.0]
|
|
c = col[i, :][image[i, :] > 0.0]
|
|
if len(r) > 1 and len(c) > 1:
|
|
H.append((r[0], c[0]))
|
|
H.append((r[-1], c[-1]))
|
|
H = np.array(H)
|
|
for j in range(0, shape[1]):
|
|
r = row[:, j][image[:, j] > 0.0]
|
|
c = col[:, j][image[:, j] > 0.0]
|
|
if len(r) > 1 and len(c) > 1:
|
|
J.append((r[0], c[0]))
|
|
J.append((r[-1], c[-1]))
|
|
J = np.array(J)
|
|
xmin = np.min([H[:, 1].min(), J[:, 1].min()])
|
|
xmax = np.max([H[:, 1].max(), J[:, 1].max()]) + 1
|
|
ymin = np.min([H[:, 0].min(), J[:, 0].min()])
|
|
ymax = np.max([H[:, 0].max(), J[:, 0].max()]) + 1
|
|
return np.array([xmin, xmax, ymin, ymax])
|
|
|
|
|
|
# Define angle and vectors operations
|
|
def vector(A, B):
|
|
"""
|
|
Define a vector from the inputed points A, B.
|
|
"""
|
|
return (B[0] - A[0], B[1] - A[1])
|
|
|
|
|
|
def determinant(u, v):
|
|
"""
|
|
Compute the determinant of 2 vectors.
|
|
--------
|
|
Inputs:
|
|
u, v : vectors (list of length 2)
|
|
Vectors for which the determinant must be computed.
|
|
----------
|
|
Returns:
|
|
determinant : real
|
|
Determinant of the input vectors.
|
|
"""
|
|
return u[0] * v[1] - u[1] * v[0]
|
|
|
|
|
|
def pseudo_angle(A, B, C):
|
|
"""
|
|
Define a pseudo-angle from the inputed points A, B, C.
|
|
"""
|
|
u = vector(A, B)
|
|
v = vector(A, C)
|
|
return determinant(u, v)
|
|
|
|
|
|
def distance(A, B):
|
|
"""
|
|
Compute the euclidian distance betwwen 2 points of the plan.
|
|
----------
|
|
Inputs:
|
|
A, B : points (list of length 2)
|
|
Points between which the distance must be computed.
|
|
----------
|
|
Returns:
|
|
distance : real
|
|
Euclidian distance between A, B.
|
|
"""
|
|
x, y = vector(A, B)
|
|
return np.sqrt(x**2 + y**2)
|
|
|
|
|
|
# Define lexicographic and composition order
|
|
def lexico(A, B):
|
|
"""
|
|
Define the lexicographic order between points A, B.
|
|
A <= B if:
|
|
- y_a < y_b
|
|
or - y_a = y_b and x_a < x_b
|
|
----------
|
|
Inputs:
|
|
A, B : points (list of length 2)
|
|
Points compared for the lexicographic order
|
|
----------
|
|
Returns:
|
|
lexico : boolean
|
|
True if A<=B, False otherwise.
|
|
"""
|
|
return (A[1] < B[1]) or (A[1] == B[1] and A[0] <= B[0])
|
|
|
|
|
|
def min_lexico(s):
|
|
"""
|
|
Return the minimum of a list of points for the lexicographic order.
|
|
----------
|
|
Inputs:
|
|
s : list of points
|
|
List of points in which we look for the minimum for the lexicographic
|
|
order
|
|
----------
|
|
Returns:
|
|
m : point (list of length 2)
|
|
Minimum for the lexicographic order in the input list s.
|
|
"""
|
|
m = s[0]
|
|
for x in s:
|
|
if lexico(x, m):
|
|
m = x
|
|
return m
|
|
|
|
|
|
def comp(Omega, A, B):
|
|
"""
|
|
Test the order relation between points such that A <= B :
|
|
This is True if :
|
|
- (Omega-A,Omega-B) is a right-hand basis of the plan.
|
|
or - (Omega-A,Omega-B) are linearly dependent.
|
|
----------
|
|
Inputs:
|
|
Omega, A, B : points (list of length 2)
|
|
Points compared for the composition order
|
|
----------
|
|
Returns:
|
|
comp : boolean
|
|
True if A <= B, False otherwise.
|
|
"""
|
|
t = pseudo_angle(Omega, A, B)
|
|
if t > 0:
|
|
return True
|
|
elif t < 0:
|
|
return False
|
|
else:
|
|
return distance(Omega, A) <= distance(Omega, B)
|
|
|
|
|
|
# Implement quicksort
|
|
def partition(s, left, right, order):
|
|
"""
|
|
Take a random element of a list 's' between indexes 'left', 'right' and place it
|
|
at its right spot using relation order 'order'. Return the index at which
|
|
it was placed.
|
|
----------
|
|
Inputs:
|
|
s : list
|
|
List of elements to be ordered.
|
|
left, right : int
|
|
Index of the first and last elements to be considered.
|
|
order : func: A, B -> bool
|
|
Relation order between 2 elements A, B that returns True if A<=B,
|
|
False otherwise.
|
|
----------
|
|
Returns:
|
|
index : int
|
|
Index at which have been placed the element chosen by the function.
|
|
"""
|
|
i = left - 1
|
|
for j in range(left, right):
|
|
if order(s[j], s[right]):
|
|
i = i + 1
|
|
temp = deepcopy(s[i])
|
|
s[i] = deepcopy(s[j])
|
|
s[j] = deepcopy(temp)
|
|
temp = deepcopy(s[i + 1])
|
|
s[i + 1] = deepcopy(s[right])
|
|
s[right] = deepcopy(temp)
|
|
return i + 1
|
|
|
|
|
|
def sort_aux(s, left, right, order):
|
|
"""
|
|
Sort a list 's' between indexes 'left', 'right' using relation order 'order' by
|
|
dividing it in 2 sub-lists and sorting these.
|
|
"""
|
|
if left <= right:
|
|
# Call partition function that gives an index on which the list will be divided
|
|
q = partition(s, left, right, order)
|
|
sort_aux(s, left, q - 1, order)
|
|
sort_aux(s, q + 1, right, order)
|
|
|
|
|
|
def quicksort(s, order):
|
|
"""
|
|
Implement the quicksort algorithm on the list 's' using the relation order
|
|
'order'.
|
|
"""
|
|
# Call auxiliary sort on whole list.
|
|
sort_aux(s, 0, len(s) - 1, order)
|
|
|
|
|
|
def sort_angles_distances(Omega, s):
|
|
"""
|
|
Sort the list of points 's' for the composition order given reference point
|
|
Omega.
|
|
"""
|
|
|
|
def order(A, B):
|
|
return comp(Omega, A, B)
|
|
|
|
quicksort(s, order)
|
|
|
|
|
|
# Define fuction for stacks (use here python lists with stack operations).
|
|
def empty_stack():
|
|
return []
|
|
|
|
|
|
def stack(S, A):
|
|
S.append(A)
|
|
|
|
|
|
def unstack(S):
|
|
S.pop()
|
|
|
|
|
|
def stack_top(S):
|
|
return S[-1]
|
|
|
|
|
|
def stack_sub_top(S):
|
|
return S[-2]
|
|
|
|
|
|
# Alignement handling
|
|
def del_aux(s, k, Omega):
|
|
"""
|
|
Returne the index of the first element in 's' after the index 'k' that is
|
|
not linearly dependent of (Omega-s[k]), first element not aligned to Omega
|
|
and s[k].
|
|
----------
|
|
Inputs:
|
|
s : list
|
|
List of points
|
|
k : int
|
|
Index from which to look for aligned points in s.
|
|
Omega : point (list of length 2)
|
|
Reference point in 's' to test alignement.
|
|
----------
|
|
Returns:
|
|
index : int
|
|
Index of the first element not aligned to Omega and s[k].
|
|
"""
|
|
p = s[k]
|
|
j = k + 1
|
|
while (j < len(s)) and (pseudo_angle(Omega, p, s[j]) == 0):
|
|
j = j + 1
|
|
return j
|
|
|
|
|
|
def del_aligned(s, Omega):
|
|
"""
|
|
Deprive a list of points 's' from all intermediary points that are aligned
|
|
using Omega as a reference point.
|
|
----------
|
|
Inputs:
|
|
s : list
|
|
List of points
|
|
Omega : point (list of length 2)
|
|
Point of reference to test alignement.
|
|
-----------
|
|
Returns:
|
|
s1 : list
|
|
List of points where no points are aligned using reference point Omega.
|
|
"""
|
|
s1 = [s[0]]
|
|
n = len(s)
|
|
k = 1
|
|
while k < n:
|
|
# Get the index of the last point aligned with (Omega-s[k])
|
|
k = del_aux(s, k, Omega)
|
|
s1.append(s[k - 1])
|
|
return s1
|
|
|
|
|
|
# Implement Graham algorithm
|
|
def convex_hull(H):
|
|
"""
|
|
Implement Graham algorithm to find the convex hull of a given list of
|
|
points, handling aligned points.
|
|
----------
|
|
Inputs:
|
|
A : list
|
|
List of points for which the the convex hull must be returned.
|
|
----------
|
|
S : list
|
|
List of points defining the convex hull of the input list A.
|
|
"""
|
|
S = empty_stack()
|
|
Omega = min_lexico(H)
|
|
sort_angles_distances(Omega, H)
|
|
A = del_aligned(H, Omega)
|
|
if len(A) < 3:
|
|
return H
|
|
else:
|
|
stack(S, A[0])
|
|
stack(S, A[1])
|
|
stack(S, A[2])
|
|
for i in range(3, len(A)):
|
|
while pseudo_angle(stack_sub_top(S), stack_top(S), A[i]) <= 0:
|
|
unstack(S)
|
|
stack(S, A[i])
|
|
return S
|
|
|
|
|
|
def image_hull(image, step=5, null_val=0.0, inside=True):
|
|
"""
|
|
Compute the convex hull of a 2D image and return the 4 relevant coordinates
|
|
of the maximum included rectangle (ie. crop image to maximum rectangle).
|
|
----------
|
|
Inputs:
|
|
image : numpy.ndarray
|
|
2D array of floats corresponding to an image in intensity that need to
|
|
be cropped down.
|
|
step : int, optional
|
|
For images with straight edges, not all lines and columns need to be
|
|
browsed in order to have a good convex hull. Step value determine
|
|
how many row/columns can be jumped. With step=2 every other line will
|
|
be browsed.
|
|
Defaults to 5.
|
|
null_val : float, optional
|
|
Pixel value determining the threshold for what is considered 'outside'
|
|
the image. All border pixels below this value will be taken out.
|
|
Defaults to 0.
|
|
inside : boolean, optional
|
|
If True, the cropped image will be the maximum rectangle included
|
|
inside the image. If False, the cropped image will be the minimum
|
|
rectangle in which the whole image is included.
|
|
Defaults to True.
|
|
----------
|
|
Returns:
|
|
vertex : numpy.ndarray
|
|
Array containing the vertex of the maximum included rectangle.
|
|
"""
|
|
H = []
|
|
shape = np.array(image.shape)
|
|
row, col = np.indices(shape)
|
|
for i in range(0, int(min(shape) / 2), step):
|
|
r1, r2 = row[i, :][image[i, :] > null_val], row[-i, :][image[-i, :] > null_val]
|
|
c1, c2 = col[i, :][image[i, :] > null_val], col[-i, :][image[-i, :] > null_val]
|
|
if r1.shape[0] > 1:
|
|
H.append((r1[0], c1[0]))
|
|
H.append((r1[-1], c1[-1]))
|
|
if r2.shape[0] > 1:
|
|
H.append((r2[0], c2[0]))
|
|
H.append((r2[-1], c2[-1]))
|
|
S = np.array(convex_hull(H))
|
|
|
|
x_min, y_min = S[:, 0] < S[:, 0].mean(), S[:, 1] < S[:, 1].mean()
|
|
x_max, y_max = S[:, 0] > S[:, 0].mean(), S[:, 1] > S[:, 1].mean()
|
|
# Get the 4 extrema
|
|
# S0 = S[x_min*y_min][np.argmin(S[x_min*y_min][:, 0])]
|
|
# S1 = S[x_min*y_max][np.argmax(S[x_min*y_max][:, 1])]
|
|
# S2 = S[x_max*y_min][np.argmin(S[x_max*y_min][:, 1])]
|
|
# S3 = S[x_max*y_max][np.argmax(S[x_max*y_max][:, 0])]
|
|
S0 = S[x_min * y_min][np.abs(0 - S[x_min * y_min].sum(axis=1)).min() == np.abs(0 - S[x_min * y_min].sum(axis=1))][0]
|
|
S1 = S[x_min * y_max][np.abs(shape[1] - S[x_min * y_max].sum(axis=1)).min() == np.abs(shape[1] - S[x_min * y_max].sum(axis=1))][0]
|
|
S2 = S[x_max * y_min][np.abs(shape[0] - S[x_max * y_min].sum(axis=1)).min() == np.abs(shape[0] - S[x_max * y_min].sum(axis=1))][0]
|
|
S3 = S[x_max * y_max][np.abs(shape.sum() - S[x_max * y_max].sum(axis=1)).min() == np.abs(shape.sum() - S[x_max * y_max].sum(axis=1))][0]
|
|
# Get the vertex of the biggest included rectangle
|
|
if inside:
|
|
f0 = np.max([S0[0], S1[0]])
|
|
f1 = np.min([S2[0], S3[0]])
|
|
f2 = np.max([S0[1], S2[1]])
|
|
f3 = np.min([S1[1], S3[1]])
|
|
else:
|
|
f0 = np.min([S0[0], S1[0]])
|
|
f1 = np.max([S2[0], S3[0]])
|
|
f2 = np.min([S0[1], S2[1]])
|
|
f3 = np.max([S1[1], S3[1]])
|
|
|
|
return np.array([f0, f1, f2, f3]).astype(int)
|