3
\$\begingroup\$

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

\$\endgroup\$
  • \$\begingroup\$ Mini-review (I haven't had the time I expected): it would be good to #undef VAL once you've finished with it. Not a problem, just a good habit to be in. (I'm surprised you didn't continue to use it for the final three statements in that function...) \$\endgroup\$ – Toby Speight Feb 28 '18 at 17:31
  • \$\begingroup\$ @TobySpeight Good remark. I forgot about the last three statements. I mainly used this macro for readability to be below the 80 chars limit. \$\endgroup\$ – nowox Feb 28 '18 at 22:29
2
\$\begingroup\$

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, ...);
\$\endgroup\$
  • \$\begingroup\$ Right, I renamed median into medfilt at last, but forgot to do it properly. I also tried to build it with full MISRA-C3 compatibility and I discovered some issues. I should have used a typedef Float32 float du to the lack of native float32_t. #pragma once is not yet a standard, but it is supported by all the compilers I am using (Keil, IAR, gcc, clang, armgcc) and again MISRA complains about it. Protecting median_init against a size of 0 is a good advise. NaN is good advise too. I should add an init value. Good code review, thanks very much :) \$\endgroup\$ – nowox Feb 28 '18 at 21:55
  • \$\begingroup\$ I don't know how to deal with Allocation responsibility in pure C. I thought about this question and came with 3 different implementations: 1. Fully static array defined with a #define KERNEL_SIZE (N), unfortunately not possible to build a modular library with that. 2. Dynamic allocation, proper medfilt_create, but requires medfilt_destroy at the end, and still issue with the allocation responsability. 3. The current implementation. This code will for sure be more elegant written in C++ and with the use of templates, however I would loose it for MCU embedded applications. \$\endgroup\$ – nowox Feb 28 '18 at 21:59
  • \$\begingroup\$ About O(n) complexity, I wanted to avoid iterating over the whole kernel to find the eldest. Your solution looks interesting and it would make me benchmark different approaches. I take your parameters order suggestion. \$\endgroup\$ – nowox Feb 28 '18 at 22:09
  • \$\begingroup\$ I updated the GitHub project with most of your suggestions. \$\endgroup\$ – nowox Feb 28 '18 at 22:28
  • 1
    \$\begingroup\$ @nowox If n is sizable, perhaps spend effort going for a log2(n) solution. More complex and more memory, yet do-able. \$\endgroup\$ – chux Feb 28 '18 at 22:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.