Return to Question

/**
 * Moving Median Filter.
 * @author CHEVALLIER Yves <[email protected]>
 */
#include <stdint.h>
#include <stdlib.h>
#include "medfilt.h"

static void swap(MedfiltNode **a, MedfiltNode **b)
 {
    // Swap two node references in the sorted table.
    MedfiltNode *temp = *a;
    *a = *b;
    *b = temp;

    // Preserve index. Used to retrive the node position in the sorted table.
    size_t index = (*a)->index;
    (*a)->index = (*b)->index;
    (*b)->index = index;
}

void median(float input, MedfiltData *data, float *median, float *min, float *max)
 {
    // New value replaces the oldest
    MedfiltNode *n = data->kernel;
    MedfiltNode *node = data->oldest;
    node->value = input;
    data->oldest = node->parent;

    // Sort the kernel
    #define VAL(x) (n[x].sorted->value)
    for (size_t i = node->index; i < data->length - 1 && VAL(i) > VAL(i + 1); i++) {
        swap(&n[i].sorted, &n[i + 1].sorted);
    }

    for (size_t i = node->index; i > 0 && VAL(i) < VAL(i - 1); i++) {
        swap(&n[i].sorted, &n[i - 1].sorted);
    }

    // Get kernel information from sorted entries
    *min = n[0].sorted->value;
    *max = n[data->length - 1].sorted->value;
    *median = n[data->length / 2].sorted->value;
}

void median_init(MedfiltData *data, MedfiltNode *nodes, size_t length)
 {
    data->kernel = nodes;
    data->length = length;

    // Linked list initialization
    data->oldest = &data->kernel[length - 1];
    for (size_t i = 0; i < length; data->oldest = &data->kernel[i], i++) {
        data->kernel[i] = (MedfiltNode) {
            .value = 0,
            .parent = data->oldest,
            .index = i,
            .sorted = &data->kernel[i]
        };
    }
}

GitHub project

#include <stdint.h>
#include <stdlib.h>
#include "medfilt.h"

static void swap(MedfiltNode **a, MedfiltNode **b) {
    // Swap two node references in the sorted table.
    MedfiltNode *temp = *a;
    *a = *b;
    *b = temp;

    // Preserve index. Used to retrive the node position in the sorted table.
    size_t index = (*a)->index;
    (*a)->index = (*b)->index;
    (*b)->index = index;
}

void median(float input, MedfiltData *data, float *median, float *min, float *max) {
    // New value replaces the oldest
    MedfiltNode *n = data->kernel;
    MedfiltNode *node = data->oldest;
    node->value = input;
    data->oldest = node->parent;

    // Sort the kernel
    #define VAL(x) (n[x].sorted->value)
    for (size_t i = node->index; i < data->length - 1 && VAL(i) > VAL(i + 1); i++) {
        swap(&n[i].sorted, &n[i + 1].sorted);
    }

    for (size_t i = node->index; i > 0 && VAL(i) < VAL(i - 1); i++) {
        swap(&n[i].sorted, &n[i - 1].sorted);
    }

    // Get kernel information from sorted entries
    *min = n[0].sorted->value;
    *max = n[data->length - 1].sorted->value;
    *median = n[data->length / 2].sorted->value;
}

void median_init(MedfiltData *data, MedfiltNode *nodes, size_t length) {
    data->kernel = nodes;
    data->length = length;

    // Linked list initialization
    data->oldest = &data->kernel[length - 1];
    for (size_t i = 0; i < length; data->oldest = &data->kernel[i], i++) {
        data->kernel[i] = (MedfiltNode) {
            .value = 0,
            .parent = data->oldest,
            .index = i,
            .sorted = &data->kernel[i]
        };
    }
}

Test program

#include <stdio.h>
#include <stdlib.h>
#include "medfilt.h"

static void kernel2str(MedfiltData *data) {
    printf("[ ");
    for (size_t i = 0; i < data->length; i++) {
        printf("%d ", (int)data->kernel[i].sorted->value);
    }
    printf("]");
}

static void display(MedfiltData *data, size_t i, int in, int min, int mid, int max) {
    printf("u[%d] = %d -> min=%d mid=%d max=%d kernel=", (int)i, in, min, mid, max);
    kernel2str(data);
    printf("\n");
}


#define KERNEL  (5)

int main(void) {
    const int in[] = { 3, 1, 6, 5, 7, 8, 9, 0, 1 };

    MedfiltData data;
    MedfiltNode nodes[KERNEL];
    medfilt_init(&data, nodes, KERNEL);

    for (size_t i = 0; i < sizeof(in) / sizeof(in[0]); i++) {
        float min, mid, max;
        medfilt(in[i], &data, &mid, &min, &max);
        display(&data, i, in[i], (int)min, (int)mid, (int)max);
    }
    return EXIT_SUCCESS;
}

GitHub project

GitHub project

https://github.com/nowox/moving-median

I am using a linked list to keep track of the oldest value in the kernel and an array of nodes sorted which is sorted.

The sorted pointed is stored along the MedfiltNode because it saves one array. It is a bit ugly though.

Regarding the sorting algorithm. It is a half-bubble-sort implementation where the value is either bubbled toward the top of the bottom of the sorted list. Average complexity is N/2. Depending on the kernel size it might worth using a binary-insertion-sort algorithm instead, but I would like to avoid to use any recursive algorithm.