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