Algorithme de la médiane roulante en C

Question

Je travaille actuellement sur un algorithme permettant d'implémenter un filtre médian glissant (analogue à un filtre moyen glissant) en C. D'après mes recherches dans la littérature, il semble y avoir deux façons raisonnablement efficaces de le faire. La première consiste à trier la fenêtre initiale de valeurs, puis à effectuer une recherche binaire pour insérer la nouvelle valeur et supprimer la valeur existante à chaque itération.

La seconde (de Hardle et Steiger, 1995, JRSS-C, Algorithme 296) construit une structure de tas à double extrémité, avec un tas maximum à une extrémité, un tas minimum à l'autre, et la médiane au milieu. Cela donne un algorithme en temps linéaire au lieu d'un algorithme qui est O(n log n).

Voici mon problème : la mise en œuvre de la première est faisable, mais je dois l'exécuter sur des millions de séries chronologiques, donc l'efficacité compte beaucoup. La seconde s'avère très difficile à mettre en œuvre. J'ai trouvé du code dans le fichier Trunmed.c du code du paquet stats de R, mais il est plutôt indéchiffrable.

Quelqu'un connaît-il une implémentation C bien écrite de l'algorithme de la médiane mobile en temps linéaire ?

Edit : Lien vers le code de Trunmed.c http://google.com/codesearch/p?hl=en&sa=N&cd=1&ct=rc#mYw3h_Lb_e0/R-2.2.0/src/library/stats/src/Trunmed.c

Answer 1

Je sais que ça fait deux jours, mais je poste juste au cas où ça aiderait.

Ce problème a été posé une fois lors d'une compétition TopCoder. Vous pouvez trouver une discussion sur plusieurs solutions possibles dans la section Match Editorial (faites défiler vers le bas jusqu'à FloatingMedian). Pour moi, l'approche multiset C++ semble la plus facile à coder.

Answer 2

J'ai regardé chez R. src/library/stats/src/Trunmed.c plusieurs fois, car je voulais quelque chose de similaire dans une classe C++ autonome ou une sous-routine C. Notez que ce sont en fait deux implémentations en une, voir src/library/stats/man/runmed.Rd (la source du fichier d'aide) qui dit

\details{
  Apart from the end values, the result \code{y = runmed(x, k)} simply has
  \code{y[j] = median(x[(j-k2):(j+k2)])} (k = 2*k2+1), computed very
  efficiently.

  The two algorithms are internally entirely different:
  \describe{
    \item{"Turlach"}{is the Härdle-Steiger
      algorithm (see Ref.) as implemented by Berwin Turlach.
      A tree algorithm is used, ensuring performance \eqn{O(n \log
        k)}{O(n * log(k))} where \code{n <- length(x)} which is
      asymptotically optimal.}
    \item{"Stuetzle"}{is the (older) Stuetzle-Friedman implementation
      which makes use of median \emph{updating} when one observation
      enters and one leaves the smoothing window.  While this performs as
      \eqn{O(n \times k)}{O(n * k)} which is slower asymptotically, it is
      considerably faster for small \eqn{k} or \eqn{n}.}
  }
}

Ce serait bien de le voir réutilisé d'une manière plus autonome. Vous êtes volontaire ? Je peux aider avec certains des bits R.

Edit 1 : Outre le lien vers l'ancienne version de Trunmed.c ci-dessus, voici les copies SVN actuelles de

• Srunmed.c (pour la version Stuetzle)
• Trunmed.c (pour la version Turlach)
• runmed.R pour la fonction R qui appelle ces

Edit 2 : Ryan Tibshirani a un peu de code C et Fortran sur binning médian rapide ce qui peut constituer un point de départ approprié pour une approche fenêtrée.

Answer 3

Je n'ai pas pu trouver une implémentation moderne d'une structure de données c++ avec une statistique d'ordre et j'ai donc fini par implémenter les deux idées dans le lien top coders suggéré par MAK ( Match Editorial : faites défiler vers le bas jusqu'à FloatingMedian).

Deux multisets

La première idée consiste à partitionner les données en deux structures de données (tas, multisets, etc.) avec O(ln N) par insertion/suppression, ce qui ne permet pas de modifier dynamiquement le quantile sans un coût important. C'est-à-dire que nous pouvons avoir une médiane glissante, ou un 75% glissant mais pas les deux en même temps.

Arbre des segments

La deuxième idée utilise un arbre de segments qui est O(ln N) pour les insertions/suppressions/requêtes mais qui est plus flexible. Le meilleur de tous, c'est que le "N" est la taille de votre plage de données. Ainsi, si votre médiane glissante a une fenêtre d'un million d'éléments, mais que vos données varient de 1 à 65536, seules 16 opérations sont nécessaires par mouvement de la fenêtre glissante de 1 million !

Le code c++ est similaire à celui que Denis a posté plus haut ("Voici un algorithme simple pour les données quantifiées")

Arbres statistiques d'ordre GNU

Juste avant d'abandonner, j'ai découvert que stdlibc++ contient des arbres statistiques d'ordre ! !!

Ceux-ci ont deux opérations critiques :

iter = tree.find_by_order(value)
order = tree.order_of_key(value)

Voir manuel libstdc++ policy_based_data_structures_test (recherche de "split and join").

J'ai emballé l'arbre pour l'utiliser dans un en-tête pratique pour les compilateurs supportant les typedefs partiels de style c++0x/c++11 :

#if !defined(GNU_ORDER_STATISTIC_SET_H)
#define GNU_ORDER_STATISTIC_SET_H
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

// A red-black tree table storing ints and their order
// statistics. Note that since the tree uses
// tree_order_statistics_node_update as its update policy, then it
// includes its methods by_order and order_of_key.
template <typename T>
using t_order_statistic_set = __gnu_pbds::tree<
                                  T,
                                  __gnu_pbds::null_type,
                                  std::less<T>,
                                  __gnu_pbds::rb_tree_tag,
                                  // This policy updates nodes'  metadata for order statistics.
                                  __gnu_pbds::tree_order_statistics_node_update>;

#endif //GNU_ORDER_STATISTIC_SET_H

Answer 4

J'ai fait un Mise en œuvre C aquí . Quelques détails supplémentaires figurent dans cette question : Médiane roulante en C - Mise en œuvre de Turlach .

Exemple d'utilisation :

int main(int argc, char* argv[])
{
   int i, v;
   Mediator* m = MediatorNew(15);

   for (i=0; i<30; i++) {
      v = rand() & 127;
      printf("Inserting %3d \n", v);
      MediatorInsert(m, v);
      v = MediatorMedian(m);
      printf("Median = %3d.\n\n", v);
      ShowTree(m);
   }
}

Answer 5

Voici un algorithme simple pour les données quantifiées (des mois plus tard) :

""" median1.py: moving median 1d for quantized, e.g. 8-bit data

Method: cache the median, so that wider windows are faster.
    The code is simple -- no heaps, no trees.

Keywords: median filter, moving median, running median, numpy, scipy

See Perreault + Hebert, Median Filtering in Constant Time, 2007,
    http://nomis80.org/ctmf.html: nice 6-page paper and C code,
    mainly for 2d images

Example:
    y = medians( x, window=window, nlevel=nlevel )
    uses:
    med = Median1( nlevel, window, counts=np.bincount( x[0:window] ))
    med.addsub( +, - )  -- see the picture in Perreault
    m = med.median()  -- using cached m, summ

How it works:
    picture nlevel=8, window=3 -- 3 1s in an array of 8 counters:
        counts: . 1 . . 1 . 1 .
        sums:   0 1 1 1 2 2 3 3
                        ^ sums[3] < 2 <= sums[4] <=> median 4
        addsub( 0, 1 )  m, summ stay the same
        addsub( 5, 1 )  slide right
        addsub( 5, 6 )  slide left

Updating `counts` in an `addsub` is trivial, updating `sums` is not.
But we can cache the previous median `m` and the sum to m `summ`.
The less often the median changes, the faster;
so fewer levels or *wider* windows are faster.
(Like any cache, run time varies a lot, depending on the input.)

See also:
    scipy.signal.medfilt -- runtime roughly ~ window size
    http://stackoverflow.com/questions/1309263/rolling-median-algorithm-in-c

"""

from __future__ import division
import numpy as np  # bincount, pad0

__date__ = "2009-10-27 oct"
__author_email__ = "denis-bz-py at t-online dot de"

#...............................................................................
class Median1:
    """ moving median 1d for quantized, e.g. 8-bit data """

    def __init__( s, nlevel, window, counts ):
        s.nlevel = nlevel  # >= len(counts)
        s.window = window  # == sum(counts)
        s.half = (window // 2) + 1  # odd or even
        s.setcounts( counts )

    def median( s ):
        """ step up or down until sum cnt to m-1 < half <= sum to m """
        if s.summ - s.cnt[s.m] < s.half <= s.summ:
            return s.m
        j, sumj = s.m, s.summ
        if sumj <= s.half:
            while j < s.nlevel - 1:
                j += 1
                sumj += s.cnt[j]
                # print "j sumj:", j, sumj
                if sumj - s.cnt[j] < s.half <= sumj:  break
        else:
            while j > 0:
                sumj -= s.cnt[j]
                j -= 1
                # print "j sumj:", j, sumj
                if sumj - s.cnt[j] < s.half <= sumj:  break
        s.m, s.summ = j, sumj
        return s.m

    def addsub( s, add, sub ):
        s.cnt[add] += 1
        s.cnt[sub] -= 1
        assert s.cnt[sub] >= 0, (add, sub)
        if add <= s.m:
            s.summ += 1
        if sub <= s.m:
            s.summ -= 1

    def setcounts( s, counts ):
        assert len(counts) <= s.nlevel, (len(counts), s.nlevel)
        if len(counts) < s.nlevel:
            counts = pad0__( counts, s.nlevel )  # numpy array / list
        sumcounts = sum(counts)
        assert sumcounts == s.window, (sumcounts, s.window)
        s.cnt = counts
        s.slowmedian()

    def slowmedian( s ):
        j, sumj = -1, 0
        while sumj < s.half:
            j += 1
            sumj += s.cnt[j]
        s.m, s.summ = j, sumj

    def __str__( s ):
        return ("median %d: " % s.m) + \
            "".join([ (" ." if c == 0 else "%2d" % c) for c in s.cnt ])

#...............................................................................
def medianfilter( x, window, nlevel=256 ):
    """ moving medians, y[j] = median( x[j:j+window] )
        -> a shorter list, len(y) = len(x) - window + 1
    """
    assert len(x) >= window, (len(x), window)
    # np.clip( x, 0, nlevel-1, out=x )
        # cf http://scipy.org/Cookbook/Rebinning
    cnt = np.bincount( x[0:window] )
    med = Median1( nlevel=nlevel, window=window, counts=cnt )
    y = (len(x) - window + 1) * [0]
    y[0] = med.median()
    for j in xrange( len(x) - window ):
        med.addsub( x[j+window], x[j] )
        y[j+1] = med.median()
    return y  # list
    # return np.array( y )

def pad0__( x, tolen ):
    """ pad x with 0 s, numpy array or list """
    n = tolen - len(x)
    if n > 0:
        try:
            x = np.r_[ x, np.zeros( n, dtype=x[0].dtype )]
        except NameError:
            x += n * [0]
    return x

#...............................................................................
if __name__ == "__main__":
    Len = 10000
    window = 3
    nlevel = 256
    period = 100

    np.set_printoptions( 2, threshold=100, edgeitems=10 )
    # print medians( np.arange(3), 3 )

    sinwave = (np.sin( 2 * np.pi * np.arange(Len) / period )
        + 1) * (nlevel-1) / 2
    x = np.asarray( sinwave, int )
    print "x:", x
    for window in ( 3, 31, 63, 127, 255 ):
        if window > Len:  continue
        print "medianfilter: Len=%d window=%d nlevel=%d:" % (Len, window, nlevel)
            y = medianfilter( x, window=window, nlevel=nlevel )
        print np.array( y )

# end median1.py