# lissage avec noyau circulaire
hist2 = np.sqrt(lissageHist(grads))
## hist2 = np.sqrt(lissageHist(grads))
## TO CHOOSE BETWEEN the two afterwards
## hist2 = np.sqrt(lissageHist2(grads))
x = (np.arange(np.size(hist2))*360./np.size(hist2))
## angle = np.ones(2)*(-999)
# detection des pics
mode = find_mode(np.abs(hist2), 0., 5., circ=True)
#--------------------
# lisser histogramme
#--------------------
def lissageHist(hist):
'''
Lisser histogram avec 4 operateurs
param/sortie: histogram des gradients
'''
# definition operators
Bx = np.array([1,2,1],float) * 1/4
Bx2 = np.array([1,0,2,0,1],float) * 1/4
Bx4 = np.array([1,0,0,0,2,0,0,0,1],float) * 1/4
Bx8 = np.array([1,0,0,0,0,0,0,0,2,0,0,0,0,0,0,0,1],float) * 1/4
taille = len(hist)
# convolution circulaire
temp = np.concatenate((hist,hist,hist)) # concat pour circularite
temp_init = temp.copy()
temp = signal.convolve(temp,Bx,'same')
temp = signal.convolve(temp,Bx2,'same')
temp2 = temp.copy()
# lissage supplementaire si pas assez de gradients presents
#if (abs(hist)>0).sum()>20:
temp = signal.convolve(temp,Bx4,'same')
temp3 = temp.copy()
#if (abs(hist)>0).sum()>30:
temp = signal.convolve(temp,Bx8,'same')
temp4 = temp.copy()
hist = temp[taille:taille*2] # prendre la partie centrale
return hist
def find_mode(tab00, smooth_fac, noise, circ=False):
tab0 = tab00.copy()
tab = np.asarray(tab0, dtype = 'float')
mm = np.max(tab)
if circ:
## duplicate and center
arg_max = np.argmax(tab)
sz = np.size(tab)
tab_o = tab.copy()
tab = np.hstack([tab, tab, tab])
# ;; TEST the input data
if min(tab) < 0:
logger.info('input data must be positive')
raise ValueError('input data must be positive')
minzero = 0
if min(tab) == 0:
minzero = 1
tab = tab + 1
sm = 0
# ;; marge in Percentage
if np.size(tab) < smooth_fac:
smooth_fac = 0
if smooth_fac > 0:
cnt = np.size(tab)
s_tab = tab.copy()
if cnt > smooth_fac:
if circ:
s_tab = utils.smooth(np.hstack([tab_o, tab, tab_o]), smooth_fac)
s_tab = s_tab[sz:np.size(tab)+sz]
else:
s_tab = utils.smooth(tab, smooth_fac)
else:
s_tab = tab.copy()
# ;; Gives all the modes without taking into account the possible
# ;; noise. Discrimation is done after
mode = find_mode_intern(np.copy(s_tab), smooth=0, circ=circ)
siz = np.size(mode.val)
final_val = np.zeros(siz) # replicate(0., siz)
final_loc = np.zeros(siz, np.int32) # replicate(0, siz)
final_beg_loc = np.zeros(siz, np.int32) # replicate(0, siz)
final_end_loc = np.zeros(siz, np.int32) # replicate(0, siz)
if noise > 0:
_val = mode.val[np.argsort(mode.loc)]
val_max = np.max(_val)
noise_abs = noise * val_max / 100.
h = 0
# while size(where(s_tab > 0), / dimensions) NE 0 do begin
while np.count_nonzero(s_tab) > 0:
mode = find_mode_intern(np.copy(s_tab), smooth=0, circ=circ)
siz = np.size(mode.val)
final_val_tmp = 0.
final_loc_tmp = 0
final_beg_loc_tmp = 0
final_end_loc_tmp = 0
s_mode_loc = np.argsort(mode.loc)
val_max = np.max(mode.val[s_mode_loc])
imax = np.argmax(mode.val[s_mode_loc])
loc_max = mode.loc[s_mode_loc][imax]
beg_loc_max = mode.beg_loc[s_mode_loc][imax]
end_loc_max = mode.end_loc[s_mode_loc][imax]
final_val_tmp = val_max
final_loc_tmp = loc_max
final_beg_loc_tmp = beg_loc_max
final_end_loc_tmp = end_loc_max
if siz >= 2:
# ;; Is this a real new mode or just tiny fluctuations?
# ;; Let's consider a noise of 10% relatively to the
# ;; highest maximum
# ;; Backward
if imax > 0:
val = mode.val[s_mode_loc][0:imax]
loc = mode.loc[s_mode_loc][0:imax]
beg_loc = mode.beg_loc[s_mode_loc][0:imax]
end_loc = mode.end_loc[s_mode_loc][0:imax]
ind = np.size(loc)
while ind >= 0:
if s_tab[final_beg_loc_tmp - 1] != 0:
dif = np.abs(s_tab[loc[ind - 1]] - s_tab[final_beg_loc_tmp])
if dif <= noise_abs:
# ;; the new mode is extended
final_beg_loc_tmp = int(np.floor(beg_loc[ind - 1]))
# ;;final_end_loc_tmp = floor(end_loc_max)
if ind >= 2:
ind = ind - 1
val = val[0:ind]
loc = loc[0:ind]
end_loc = end_loc[0:ind]
beg_loc = beg_loc[0:ind]
else:
break
else:
break
else:
break
if imax < siz - 1:
val = mode.val[s_mode_loc][imax + 1:siz]
loc = mode.loc[s_mode_loc][imax + 1:siz]
beg_loc = mode.beg_loc[s_mode_loc][imax + 1:siz]
end_loc = mode.end_loc[s_mode_loc][imax + 1:siz]
ind = 0
n = np.size(loc)
while n >= 1:
if s_tab[final_end_loc_tmp + 1] != 0:
dif = abs(s_tab[loc[0]] - s_tab[final_end_loc_tmp])
if dif <= noise_abs:
# ;; the new mode is extended
final_end_loc_tmp = end_loc[0]
if n >= 2:
val = val[1:n]
loc = loc[1:n]
end_loc = end_loc[1:n]
beg_loc = beg_loc[1:n]
n = np.size(loc)
else:
break
else:
break
else:
break
if np.size(final_end_loc_tmp) > 1:
final_end_loc_tmp = final_end_loc_tmp[0]
if np.size(final_beg_loc_tmp) > 1:
final_beg_loc_tmp = final_beg_loc_tmp[0]
final_beg_loc[h] = final_beg_loc_tmp
final_end_loc[h] = final_end_loc_tmp
final_val[h] = final_val_tmp
final_loc[h] = final_loc_tmp
h = h + 1
# ;; We now repeat the operation for all the remaining modes
# pt = indgen(final_end_loc_tmp-final_beg_loc_tmp+1)+final_beg_loc_tmp
pt = np.arange(final_end_loc_tmp - final_beg_loc_tmp + 1) + final_beg_loc_tmp
if pt[0] > 0:
if s_tab[pt[0] - 1] > 0:
pt = pt[1:np.size(pt)]
nn = np.size(pt)
if pt[nn - 1] < np.size(s_tab) - 1:
if s_tab[pt[nn - 1] + 1] > 0:
pt = pt[0:np.size(pt) - 1]
s_tab[pt] = 0.
final_beg_loc = final_beg_loc[0:h]
final_end_loc = final_end_loc[0:h]
final_val = final_val[0:h]
final_loc = final_loc[0:h]
else:
final_beg_loc = mode.beg_loc
final_end_loc = mode.end_loc
final_val = mode.val
final_loc = mode.loc
s_final_loc = np.argsort(final_loc)
# ;; Elements are given by increasing location
final_beg_loc = final_beg_loc[s_final_loc]
final_end_loc = final_end_loc[s_final_loc]
final_val = final_val[s_final_loc]
final_loc = final_loc[s_final_loc]
res = PropertyDict(val=final_val,
loc=final_loc,
beg_loc=final_beg_loc,
end_loc=final_end_loc)
if circ:
## Arrange outputs:
wres = np.where((res.loc >= sz) & (res.loc < 2*sz))[0]
if np.size(wres) >= 1:
final_loc = np.mod(res.loc[wres], sz)
final_val = res.val[wres]
final_beg_loc = np.mod(res.beg_loc[wres], sz)
final_end_loc = np.mod(res.end_loc[wres], sz)
final_val = res.val[wres]
wres_min = np.argmin(wres)
wres_max = np.argmax(wres)
final_beg_loc[0] = np.mod(res.end_loc[wres[wres_min]-1], sz)
final_end_loc[-1] = np.mod(res.beg_loc[wres[wres_max]+1], sz)
res = PropertyDict(val=final_val,
loc=final_loc,
beg_loc=final_beg_loc,
end_loc=final_end_loc)
if minzero == 1:
res.val = res.val - 1
# ;;; To a check
for a in range(0, np.size(res.val) - 1):
if res.end_loc[a] != res.beg_loc[a + 1]:
logger.error('Bug in the system - check the find_mode analysis')
#print 'Bug in the system - check the find_mode analysis'
raise AssertionError('Bug in the system - check the find_mode analysis')
return res
def smooth(v1, w):
"""
try to give the same result than the IDL smooth function.
Maybe some difference in the result for pair values
:param v1: array to smooth
:param w: width of the smoothing window
:return: smoothed array
"""
if (w < 1):
w = 1
else:
w = int(math.floor(w)) # we handle the case where the value is not an integer
wa = np.ones(w) / float(w)
return ndimage.convolve1d(v1, wa)
def find_mode_intern(hist, smooth=None, circ=False):
if smooth > 0:
hist = utils.smooth(hist, smooth)
hist_tmp = np.asarray(hist)
Opt = []
Oval = []
cond = np.size(np.argwhere(hist != 0)) > 0
while cond: # .ma.masked_equal(hist_tmp,0).count()>0:
m = max(hist_tmp)
i = np.argmax(hist_tmp)
# ;; The max of the mode is found, let's find its extension
pt_tmp = np.zeros(np.size(hist_tmp), np.int32) # replicate(0, n_elements(hist_tmp))
pt_tmp[0] = i
val_tmp = np.zeros(np.size(hist_tmp)) # replicate(0., n_elements(hist_tmp))
val_tmp[0] = m
# ;; Backward
h = 1
iloc = i
out = search_mode(hist_tmp, pt_tmp, val_tmp, -1, h, iloc)
h = out.h # in these case we need to get back the h
# ;; Reverse order
val_tmp = utils.invert_tab(out.val, h)
pt_tmp = utils.invert_tab(out.pt, h)
# ;; Forward
out = search_mode(hist_tmp, pt_tmp, val_tmp, 1, out.h, i)
pt = out.pt[0:out.h]
val = out.val[0:out.h]
Opt.append(pt)
Oval.append(val)
if pt[0] > 0:
if hist_tmp[pt[0] - 1] > 0:
pt = pt[1:np.size(pt)]
nn = np.size(pt)
if pt[nn - 1] < np.size(hist_tmp) - 1:
if hist_tmp[pt[nn - 1] + 1] > 0:
pt = pt[0:np.size(pt) - 1] # -2
for p in pt:
hist_tmp[p] = 0.
cond = np.size(np.argwhere(hist != 0)) > 0
value = np.zeros(len(Oval)) # replicate(0., Oval->n_elements())
location = np.zeros(len(Oval), np.int32) # replicate(0, Oval->n_elements())
beg_location = np.zeros(len(Oval), np.int32) # replicate(0, Oval->n_elements())
end_location = np.zeros(len(Oval), np.int32) # replicate(0, Oval->n_elements())
# ;;mean = replicate(0., Oval->n_elements())
# ;;std_dev = replicate(0., Oval->n_elements())
# ;;Pearson_md_sk = replicate(0., Oval->n_elements())
for a in range(0, len(Oval)):
value[a] = np.max(Oval[a])
wmax = np.argmax(Oval[a])
location[a] = Opt[a][wmax]
beg_location[a] = Opt[a][0]
end_location[a] = Opt[a][np.size(Opt[a]) - 1]
# ;; result is sorted by order of location
s_location = np.argsort(location)
value2 = value[s_location]
location2 = location[s_location]
beg_location2 = beg_location[s_location]
end_location2 = end_location[s_location]
out = PropertyDict(val=value2,
loc=location2,
beg_loc=beg_location2,
end_loc=end_location2)
return out
def search_mode(hist, pt, val, increm, h, iloc):
"""
This function look for all the values and points to associate to a
mode. Backward if increm = -1 and Forward if increm = +1
:param hist: ensemble des points non scanné pour ceux différent de zero
:param pt: indices de l'histogramme appartenant au mode
:param val: valeur de l'histogramme du mode
:param increm: sens de recherche
:param h: position dans le mode
:param iloc: position dans l'histogramme
:return:
"""
cond = ((iloc + increm) >= 0) and ((iloc + increm) <= (np.size(hist) - 1))
while cond:
if (hist[iloc + increm] > 0.) and (hist[iloc + increm] <= hist[iloc]):
pt[h] = iloc + increm
val[h] = hist[iloc + increm]
h = h + 1
iloc = iloc + increm
cond = ((iloc + increm) >= 0) and ((iloc + increm) <= np.size(hist) - 1)
else:
cond = 0
out = PropertyDict(pt=pt, val=val, h=h)
return out
def invert_tab(tab, i):
"""
;; invert the order of the i first values of the tab
:param tab:
:param i:
:return:
"""
# n = np.size(tab)
# tab_tmp = np.zeros(n)#replicate(0., n )
tab_tmp = np.copy(tab)
tab_tmp[:] = 0
for r in range(0, i):
tab_tmp[r] = tab[i - 1 - r]
return tab_tmp
class PropertyDict(dict):
"""
classe se comportant comme un dictionnaire et en plus comme une structure
>>> a = PropertyDict(A=12, B=32,C='toto')
>>> a['A']
12
>>> a.A
12
>>> a.C='titi'
>>> a['C']
'titi'
"""
def __init__(self, *args, **kw):
super(PropertyDict, self).__init__(*args, **kw)
self.__dict__ = self
def __getstate__(self):
"""
pour faire du peakle
"""
return self
def __setstate__(self, state):
"""
pour faire du peakle
"""
self.update(state)
self.__dict__ = self
if __name__ == '__main__':
import doctest
doctest.testmod()