Commit d667ba64 authored by Félix Hartmann's avatar Félix Hartmann
Browse files

Better handling of skeleton detection stepsize and curvature computation.

parent 1055c3dd
This diff is collapsed.
......@@ -642,7 +642,8 @@ def get_steady_state(hdf5file):
def save_detection_step(hdf5file, step):
"""
Fonction pour sauvegarder le pas d'espace pour la détection des squelettes.
Save the spatial step size (in pixels) used in the algorithm of
skeleton detection.
"""
if step is None:
goods = None
......@@ -654,7 +655,8 @@ def save_detection_step(hdf5file, step):
def get_detection_step(hdf5file):
"""
Fonction pour récupérer le pas d'espace pour la détection des squelettes.
Return the spatial step size (in pixels) used in the algorithm of
skeleton detection.
"""
step = None
with h5py.File(hdf5file, 'r') as f:
......
......@@ -1049,20 +1049,6 @@ def load_results(file_name):
# TRAITEMENT DES TIGES EXTRAITES #
###############################################################################
def moving_avg_historical(x, W):
# no longer used
W = int(W)
# S = array( [ mean( x[i-W/2:i+W/2] ) for i in arange(W/2, len(x)-W/2)] )
xT = ma.hstack([x[W / 2 - 1:0:-1], x, x[-1:-W / 2 - 1:-1]])
wind = ones(W)
# Comptibilitée version de numpy
try:
S = ma.convolve(wind / wind.sum(), xT, mode='valid')
except:
S = convolve(wind / wind.sum(), xT, mode='valid')
return S
def moving_average(x, window_width=11, window='hanning'):
"""Smooth the data using a window with requested size.
......@@ -1096,9 +1082,6 @@ def moving_average(x, window_width=11, window='hanning'):
return y[(window_width/2-1):-(window_width/2)] instead of just y.
"""
if x.ndim != 1:
raise ValueError("Smooth only accepts 1 dimension arrays.")
if x.size < window_width:
raise ValueError("Input vector needs to be bigger than window width.")
......@@ -1110,15 +1093,48 @@ def moving_average(x, window_width=11, window='hanning'):
raise ValueError("""Window is on of 'flat', 'hanning', 'hamming', """
"""'bartlett', 'blackman'""")
s = numpy.r_[x[window_width-1:0:-1], x, x[-2:-window_width-1:-1]]
if window == 'flat':
w = ones(window_width, 'd')
else:
w = eval('numpy.'+window+'(window_width)')
y = numpy.convolve(w/w.sum(), s, mode='valid')
# return y
return y[(window_width//2-1):-(window_width//2)]
if x.ndim == 1:
s = numpy.r_[x[window_width-1:0:-1], x, x[-2:-window_width-1:-1]]
y = convolve(w/w.sum(), s, mode='valid')
return y[(window_width//2-1):-(window_width//2)]
elif x.ndim == 2:
y = zeros_like(x) - 3000
for i in len(x): # loop on rows
si = numpy.r_[x[i,window_width-1:0:-1], x[i], x[i,-2:-window_width-1:-1]]
y[i] = convolve(
w/w.sum(), si, mode='valid')[(window_width//2-1):-(window_width//2)]
return y
else:
raise ValueError("Smooth only accepts 1 dimension arrays.")
def compute_moving_window_size(tiges_data, tige_id, hdf5file, image=None, mask=None):
"""Compute and return the size of the moving window for smoothing.
Parameters:
-----------
tiges_data: TigesManager instance
tige_id: int, index of the tige
hdf5file: path to hdf5 data file
image: int, for restricting the smoothing to a single image, default None
mask: a mask array for excluding images, default None
Not compatible with the 'image' parameter.
"""
diameters = tiges_data.diam[tige_id]
if image is not None:
diameters = diameters[image]
elif mask is not None:
diameters = diameters[mask]
# Retrieve the detection step from the h5 file
detection_step = h5store.get_detection_step(hdf5file)
# Size of the averaging window
W = int(4 * round(diameters.mean() / detection_step))
return W
def traite_tige2(xc, yc, R, pas, cutoff=5):
......@@ -1246,7 +1262,8 @@ def traite_tiges2(tiges, itige=None, pas=0.3, causette=True,
bad_value = -30000
Ntime = len(tiges.xc[0, :, 0])
Max_tige_lengh = len(tiges.xc[0, 0, :]) # Taille de la convolution possible
Max_tige_lengh = len(tiges.xc[0, 0, :]) + 1 # Taille de la convolution possible
# print("Max_tige_lengh:", Max_tige_lengh)
x = ma.masked_equal(zeros((Ntige, Ntime, Max_tige_lengh)) + bad_value, bad_value)
y = ma.masked_equal(zeros((Ntige, Ntime, Max_tige_lengh)) + bad_value, bad_value)
s = ma.masked_equal(zeros((Ntige, Ntime, Max_tige_lengh)) + bad_value, bad_value)
......@@ -1382,28 +1399,131 @@ def get_arc_length(tiges_data, tige, image, cutoff=5, pas=0.3, scale=None):
return s
def get_tige_curvature(angles, s, smoothing=False, window_width=30):
"""Return the curvature along a tige.
def get_tige_curvatures_for_one_image(
angles, s, angle_smoothing=False, window_width=30, return_angles=False):
"""Return the curvatures along an organ for a single image.
If 'smoothing' is True, 'angles' is smoothed with a moving window
of width 'window_width'.
Parameters:
-----------
angles: 1-D array of floats
Angles along the organ skeleton.
s: 1-D array of floats
Curvilinear abscissa along the organ skeleton.
angle_smoothing: Bool, default False
If 'angle_smoothing' is True, 'angles' is smoothed with a
moving window of width 'window_width'.
window_width: int, default 30
Width of the window used for smoothing angles,
if 'angle_smoothing' is True.
return_angles: Bool, default False
If True, return the angles as a 1-D masked array, smoothed if
'angle_smoothing' is true.
"""
if smoothing:
if angle_smoothing:
smooth_angles = moving_average(angles, window_width=window_width)
else:
smooth_angles = angles
# Computing the gradient with central differences
y1 = hstack((smooth_angles[0], smooth_angles[:-1]))
y2 = hstack((smooth_angles[1:], smooth_angles[-1]))
dx1 = hstack((0, diff(s)))
dx2 = hstack((diff(s), 0))
curvatures = (y2 - y1) / (dx2 + dx1)
if smoothing:
curvatures = gradient(smooth_angles, s)
if return_angles:
return curvatures, smooth_angles
else:
return curvatures
def get_tige_curvatures(angles, s, radii=None, angle_smoothing=False, window_width=30,
trim_left=0, trim_right=0, mask=None,
return_angles=False, return_s=False, return_radii=False):
"""Return the curvatures along an organ for all images.
Parameters:
-----------
angles: 2-D array of floats
Angle along the organ skeleton, for all images
s: 2-D array of floats
Curvilinear abscissa along the organ skeleton, for all images.
radii; 2-D array of floats, optional, default None.
radius along the organ skeleton, for all images.
angle_smoothing: Bool, default False
If 'smoothing' is True, 'angles' is smoothed with a moving
window of width 'window_width'.
window_width: int, default 30
Width of the window used for the smoothing.
trim_left: int, default 0
Number of points to trim on the left of the arrays.
trim_right: int, default 0
Number of points to trim on the right of the arrays.
mask: a mask array for excluding images, default None
return_angles: Bool, default False
If True, return the angles as a 2-D masked array, smoothed if
'angle_smoothing' is True, with the same shape as the
curvature array.
return_s: Bool, default False
If True, return the organ curvilinear abscissa as a 2-D
masked array, with the same shape as the curvature array.
return_radii: Bool, default False
If True and 'radii' given, return the organ radii as a 2-D
masked array, with the same shape as the curvature array.
"""
curvatures = zeros_like(angles) - 3000 # curvatures
if return_angles:
new_angles = zeros_like(angles) - 3000
if return_s:
new_s = zeros_like(angles) - 3000
if return_radii and radii is not None:
new_radii = zeros_like(angles) - 3000
for i in range(len(angles)): # loop on images
angles_i_bool = angles[i].mask == False
if not angles_i_bool.any(): # if no measurement for image i
continue
imax = flatnonzero(angles_i_bool).max() + 1
angles_i = angles[i, trim_left:imax-trim_right]
s_i = s[i, trim_left:imax-trim_right]
C_i, smooth_angles = get_tige_curvatures_for_one_image(
angles_i, s_i, angle_smoothing=True, window_width=window_width,
return_angles=True)
curvatures[i, :len(C_i)] = C_i
if return_angles:
new_angles[i, :len(smooth_angles)] = smooth_angles
if return_s:
new_s[i, :len(s_i)] = s_i
if return_radii and radii is not None:
R_i = radii[i, trim_left:imax-trim_right]
new_radii[i, :len(R_i)] = R_i
curvatures = ma.masked_equal(curvatures, -3000)
curvatures = ma.masked_invalid(curvatures)
if mask is None:
mask = ones(len(curvatures), dtype=bool)
curvatures = curvatures[mask]
arrays_to_return = [curvatures]
if return_angles:
new_angles = ma.masked_equal(new_angles, -3000)
new_angles = ma.masked_invalid(new_angles)
new_angles = new_angles[mask]
arrays_to_return.append(new_angles)
if return_s:
new_s = ma.masked_equal(new_s, -3000)
new_s = ma.masked_invalid(new_s)
new_s = new_s[mask]
arrays_to_return.append(new_s)
if return_radii:
new_radii = ma.masked_equal(new_radii, -3000)
new_radii = ma.masked_invalid(new_radii)
new_radii = new_radii[mask]
arrays_to_return.append(new_radii)
if len(arrays_to_return) == 1:
return curvatures
else:
return tuple(arrays_to_return)
def get_tiges_lines(xc, yc):
"""
Fonction qui retourne un liste pour chaque tiges une liste de
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment