Moving median filter (medfilt) in C

Question

I have implemented a Moving Median Filter and I would like a review.

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.

Header

/**
* Moving Median Filter.
*
* This algorithm is iterative. Each call will compute the next point.
* In the example below, the kernel has a size of 3. Notice that the
* values in the kernel are alway sorted. The left value is therefore
* the minimum in the kernel, the center value is the median and the
* right value is the maximum value.
*
* Input data:         0 1 2 5 4 8 3
* Kernel:      |0 0 0|            .   (min=0, med=0, max=0)
*                |0 0 0|          .   (min=0, med=0, max=0)
*                  |0 0 1|        .   (min=0, med=0, max=1)
*                    |0 1 2|      .   (min=0, med=1, max=2)
*                      |1 2 5|    .   (min=1, med=2, max=5)
*                        |2 4 5|  .   (min=2, med=4, max=5)
*                          |4 5 8|.   (min=4, med=5, max=8)
*                            |3 4 8|  (min=3, med=4, max=8)
*/
#pragma once
#include <stdlib.h>

typedef struct MedfiltNode {
    float value;
    size_t index; // Node index in the sorted table
    struct MedfiltNode* parent;
    struct MedfiltNode* sorted;
} MedfiltNode;

typedef struct MedfiltData {
    MedfiltNode *kernel; // Working filter memory
    MedfiltNode *oldest; // Reference to the oldest value in the kernel
    size_t length; // Number of nodes
} MedfiltData;

void median_init(MedfiltData *data, MedfiltNode *nodes, size_t length);
void median(float input, MedfiltData *data, float *median, float *min, float *max);

Source Code

#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

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

Answer 1

Presentation

Code formatted and presented well. Refreshing. Compiled properly - almost.

Initialization

With median_init(), the queue is pre-populated with the value of 0. This implies that the return values of min,median,max are all bogus until n values added. it would make sense to either

1. Allow an initialization value: median_init(...., float init_value)

2. Or better, functions on min,median,max as if it had less than n elements are in it until n elements are added.

Complexity

Instead of the structures involved, this task could be solved with less data usage: a length, eldest index and a single array as below and still maintain O(n) complexity. Adding is a simple as scanning for the eldest index in a sorted array, replacing it with the new value. Then bubble sort the new value into place.

typedef struct MedfiltDataAlt {
    float value;
    size_t queue_order;
} MedfiltDataAlt;

Allocation responsibility

This code obliges the calling code to provide the space needed. This fine. Yet I find it more useful for the median_...() functions to handle that. OP's choice on this, just wanted to provide an alternative POV.

Added functionality

A means to report the min,median,max without adding a new number would be useful.

Parameter order

I'd expect the state data first.

// void median(float input, MedfiltData *data, float *median, float *min, float *max)
   void median(MedfiltData *data, float input, float *median, float *min, float *max)

NAN

As float can take on a value of not-a-number. Code's use of > and < does not clearly push the NAN is a proscribed manner. When NANs are involved, > is not the opposite of <=. If code was amended to deal with a growing array during the first n calls of median(), consider a median(NAN, ...) to have the effect of shortening the queue and not saving the NAN.

Naming

The .c file (assume medfilt.c to match medfilt.h) name, the type names and function names would be more uniform if they all began with the same prefix and case rather than Medfilt, median and medfilt.

This is manifested in the minor code error in main.c

// medfilt_init(&data, nodes, KERNEL);
median_init(&data, nodes, KERNEL);

Portability

#pragma once is not defined in the C standard

Function at 0

Minor: Both median() and median_init() fail with a size of 0. I would expect code to tolerate this edge case and not cause undefined behavior.

display()

Minor: I'd expect code to use "%zu" to match size_t or a cast to a wide unsigned type to maintain same sign-ness. Also int min, int mid, int max --> float min, float mid, float max and use "%g" or "%e". A goal of display is diagnostic. It should avoid data truncation.

// printf("u[%d] =...", (int)i, ...);
printf("u[%zu] =...", i, ...);
printf("u[%lu] =...", (unsigned long)i, ...);