Tiny (yet fast) median stacking/blending filter, with minimal dependency (stb_image and/or libjpeg only!)

median.c

/*
 * MIT License
 *
 * Copyright (c) 2021 Davidson Francis <[email protected]>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

/*
 * median.c supports both stb_image and libjpeg(-turbo) libraries
 *
 * Compile with -DUSE_STB_IMAGE to use stb_image and with -DUSE_LIBJPEG
 * to use libjpeg(-turbo). In the absence of both, stb_image will be
 * used.
 *
 * stb_image instructions (default build):
 * ---------------------------------------
 *
 * Build with: gcc median.c -o median -O3 -fopenmp -lm
 *
 * libjpeg(-turbo) instructions:
 * -----------------------------
 *
 * Build with: gcc median.c -o median -O3 -fopenmp -ljpeg -DUSE_LIBJPEG
 *
 * Common for both libraries:
 * --------------------------
 *
 * If you _do not_ want OpenMP and/or want to reduce dependencies
 * just remove -fopenmp from your command-line options.
 *
 * You may also want to reduce binary size by passing the flags:
 * -fdata-sections -ffunction-sections -Wl,--gc-sections
 *
 * with this, file size reduces from 123k to 79k or a 35.7% reduction
 * for dynamic linking.
 *
 * Dependencies:
 * =============
 *
 * stb_image
 * ---------
 *
 * This file just requires a c99-compatible C compiler and 2 small header-only
 * libraries: stb_image.h and stb_image_write.h. Both can be obtained from:
 *   https://github.com/nothings/stb
 *
 * or more specifically:
 *   https://raw.githubusercontent.com/nothings/stb/master/stb_image.h
 *   https://raw.githubusercontent.com/nothings/stb/master/stb_image_write.h
 *   (versions used: v2.26 and v1.15, respectively, although I do not think
 *    that will have a huge break API in future versions).
 *
 * stb_image supports other file formats besides JPEG, but for simplicity,
 * I'm using only JPEG. To support other formats as well (such as PNG, BMP...),
 * comment the line '#define STBI_ONLY_JPEG' and/or define other macros. More
 * info about this inside the stb_image header.
 *
 * or libjpeg(-turbo):
 * -------------------
 *
 * median.c requires stb_image or libjpeg, not both!. Please check the specific
 * procedures of your Linux distribution for installing libjpeg.
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>

#if !defined(USE_STB_IMAGE) && !defined(USE_LIBJPEG)
#	define USE_STB_IMAGE
#endif

#if defined(USE_STB_IMAGE)
#	define FREE_IMAGE  stbi_image_free
#	define WRITE_IMAGE stbi_write_jpg
#	define READ_IMAGE  stbi_load
#	define STB_IMAGE_IMPLEMENTATION
#	define STBI_ONLY_JPEG
#	include "stb_image/stb_image.h"
#	define STB_IMAGE_WRITE_IMPLEMENTATION
#	include "stb_image/stb_image_write.h"
#else
#	define FREE_IMAGE  free
#	define WRITE_IMAGE write_jpeg_turbo
#	define READ_IMAGE  load_jpeg_turbo
#	include <jpeglib.h>
#	include <inttypes.h>
#endif

/* Pixel utilities. */
#define M_r(i,j) \
	( ((((i)*WIDTH)+(j))*CHANNELS)+0 )
#define M_g(i,j) \
	( ((((i)*WIDTH)+(j))*CHANNELS)+1 )
#define M_b(i,j) \
	( ((((i)*WIDTH)+(j))*CHANNELS)+2 )

uint8_t  *out;
uint8_t **imgs;
int img_count;
int WIDTH, HEIGHT, CHANNELS;

/**
 * @brief Exhibits a formatted message on stderr and aborts
 * the program.
 *
 * @param fmt Formatted message, just like printf.
 */
static inline void panic(const char *fmt, ...)
{
	va_list args;

	va_start(args, fmt);
	vfprintf(stderr, fmt, args);
	va_end(args);

	/* Free resources. */
	if (img_count)
	{
		if (imgs)
		{
			for (int i = 0; i < img_count; i++)
				if (imgs[i])
					FREE_IMAGE(imgs[i]);
			free(imgs);

			if (out)
				free(out);
		}
	}

	exit(EXIT_FAILURE);
}

#if defined(USE_LIBJPEG)

/**
 * @brief For a given @p filename, reads a compressed
 * JPEG file and outputs a decompressed one.
 *
 * @param filename File to be read.
 * @param width Returned image width.
 * @param height Returned image height.
 * @param channels Returned image channels.
 * @param unused Required channels (unused).
 *
 * @return Returns a decompressed JPEG image or NULL
 * otherwise.
 */
uint8_t* load_jpeg_turbo(const char *filename,
	int *width, int *height, int *channels, int unused)
{
	struct jpeg_decompress_struct cinfo;
	struct jpeg_error_mgr jerr;
	int row_stride;
	uint8_t *buf;
	uint8_t *out;
	FILE *fp;

	((void)unused);

	/* Allocate and initialize a JPEG decompression object. */
	cinfo.err = jpeg_std_error(&jerr);
	jpeg_create_decompress(&cinfo);

	/* Set input file. */
	fp = fopen(filename, "rb");
	if (fp == NULL)
		return (NULL);
	jpeg_stdio_src(&cinfo, fp);

	/* Read JPEG header. */
	jpeg_read_header(&cinfo, 1);

	/* Start decompressor. */
	jpeg_start_decompress(&cinfo);

	/* Set our initial parameters. */
	*width = cinfo.output_width;
	*height = cinfo.output_height;
	*channels = cinfo.output_components;

	/* Allocate our buffer. */
	out = malloc(cinfo.output_width * cinfo.output_height *
		cinfo.output_components);
	if (!out)
	{
		fclose(fp);
		return (NULL);
	}

	row_stride = cinfo.output_width * cinfo.output_components;

	while (cinfo.output_scanline < cinfo.output_height)
	{
		/*
		 * One scanline per time.
		 */
		buf = ((out + (row_stride * cinfo.output_scanline)));
		jpeg_read_scanlines(&cinfo, &buf, 1);
	}

	jpeg_finish_decompress(&cinfo);
	jpeg_destroy_decompress(&cinfo);
	fclose(fp);

	return (out);
}

/**
 * @brief For a given image, writes into @p filename, with
 * the given @p width, @p height, @p channels and @p quality.
 *
 * @param filename File to be written.
 * @param width Image width.
 * @param height Image height.
 * @param channels Image channels.
 * @param buf Image buffer to be written.
 * @param quality Image quality.
 *
 * @return Returns 1 if success, 0 otherwise.
 */
int write_jpeg_turbo(const char *filename, int width, int height,
	int channels, uint8_t *buf, int quality)
{
	struct jpeg_compress_struct cinfo;
	struct jpeg_error_mgr jerr;
	uint32_t stride;
	uint8_t *row;

	cinfo.err = jpeg_std_error(&jerr);
	jpeg_create_compress(&cinfo);

    FILE *fp = fopen(filename, "wb");
    if (fp == NULL)
        return (0);

    /* Set output file. */
    jpeg_stdio_dest(&cinfo, fp);

    /* Configs. */
    cinfo.image_width      = width;
    cinfo.image_height     = height;
    cinfo.input_components = channels;
    cinfo.in_color_space   = JCS_RGB;
    jpeg_set_defaults(&cinfo);
    jpeg_set_quality(&cinfo, quality, 1);

    /* Start compression. */
    jpeg_start_compress(&cinfo, 1);

    /* Write each row to the file. */
    row = buf;
    stride = width * channels;

    for (int i = 0; i < height; i++)
    {
        jpeg_write_scanlines(&cinfo, &row, 1);
        row += stride;
    }

    /* Finish. */
    jpeg_finish_compress(&cinfo);
    jpeg_destroy_compress(&cinfo);
    fclose(fp);

    return (1);
}
#endif

/**
 * @brief Quick select algorithm: for a given array @p arr
 * of size @p n, finds the median element without ordering
 * everything.
 *
 * In contrast to what I had said before (previous revisions
 * of the code), there are multiple ways to handle median for
 * an even number of elements. The 'classic', consists of
 * getting the average of the two central values (which
 * implies a great overhead for quickselect and that made me
 * abandon the idea at first). However, in [1] it is said
 * that in addition to the center mean, two other forms are
 * the 'lower median' and 'upper median', which correspond
 * to the values before and after the center and that both
 * are equivalent.
 *
 * In tests performed with insertion sort (center average)
 * and quickselect, they demonstrated that the resulting
 * image, although not exactly identical, has a large noise
 * removal. So I see no reason not to adopt quickselect.
 *
 * This algorithm was obtained from [2] and is licensed as
 * a public domain.
 *
 * References:
 * [1]: Thomas H. Cormen, Charles E. Leiserson,
 * Ronald L. Rivest, Clifford Stein, Introduction to algorithms,
 * MIT press.
 *
 * [2]: Fast median search: an ANSI C implementation
 * Nicolas Devillard - ndevilla AT free DOT fr, July 1998
 * http://ndevilla.free.fr/median/median/index.html
 *
 * @param arry RGB list to be sorted.
 * @param n List size.
 *
 * @return Returns the median value.
 */
#define ELEM_SWAP(a,b) \
	{ register uint8_t t=(a);(a)=(b);(b)=t; }

static inline uint8_t quick_select(uint8_t *arr, int n)
{
	int low, high;
	int median;
	int middle, ll, hh;

	low = 0; high = n-1; median = (low + high) / 2;
	for (;;)
	{
		/* One element only */
		if (high <= low)
			return arr[median];

		/* Two elements only */
		if (high == low + 1)
		{
			if (arr[low] > arr[high])
				ELEM_SWAP(arr[low], arr[high]);
			return arr[median];
		}

		/*
		 * Find median of low, middle and high items;
		 * swap into position low
		 */
		middle = (low + high) / 2;
		if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]);
		if (arr[low] > arr[high])    ELEM_SWAP(arr[low], arr[high]);
		if (arr[middle] > arr[low])  ELEM_SWAP(arr[middle], arr[low]);

		/*
		 * Swap low item (now in position middle) into
		 * position (low+1)
		 */
		ELEM_SWAP(arr[middle], arr[low+1]);

		/*
		 * Nibble from each end towards middle, swapping items
		 * when stuck
		 */
		ll = low + 1;
		hh = high;
		for (;;)
		{
			do ll++; while (arr[low] > arr[ll]);
			do hh--; while (arr[hh]	 > arr[low]);

			if (hh < ll)
				break;

			ELEM_SWAP(arr[ll], arr[hh]);
		}

		/*
		 * Swap middle item (in position low) back into correct
		 * position.
		 */
		ELEM_SWAP(arr[low], arr[hh]);

		/* Re-set active partition */
		if (hh <= median)
			low = ll;
		if (hh >= median)
			high = hh - 1;
	}
}
#undef ELEM_SWAP

/**
 * @brief For a given coordinate, obtains the median for the loaded
 * images.
 *
 * @param x X-Coordinate.
 * @param y Y-Coordinate.
 * @param r Pointer to the Red-component, this will be updated with
 *          the median value found.
 * @param g Pointer to the Green-component, this will be updated with
 *          the median value found.
 * @param b Pointer to the Blue-component, this will be updated with
 *          the median value found.
 */
static inline void get_pixel_median(int x, int y, uint8_t *r,
	uint8_t *g, uint8_t *b)
{
	/*
	 * VLAs... yeah, but I do not expect img_count to be *huge*,
	 * right? just 3 bytes per input, seems feasible... and
	 * dinamically allocate memory here just do not seems right.
	 */
	uint8_t mtx_r[img_count];
	uint8_t mtx_g[img_count];
	uint8_t mtx_b[img_count];

	for (int i = 0; i < img_count; i++)
	{
		mtx_r[i] = imgs[i][M_r(x,y)];
		mtx_g[i] = imgs[i][M_g(x,y)];
		mtx_b[i] = imgs[i][M_b(x,y)];
	}

	*r = quick_select(mtx_r, img_count);
	*g = quick_select(mtx_g, img_count);
	*b = quick_select(mtx_b, img_count);
}

/**
 * @brief Calculates the median =).
 */
static void do_image_median(void)
{
	uint8_t r, g, b;

	#pragma omp parallel for private(r,g,b)
	for (int i = 0; i < HEIGHT; i++)
	{
		for (int j = 0; j < WIDTH; j++)
		{
			get_pixel_median(i, j, &r, &g, &b);
			out[M_r(i,j)] = r;
			out[M_g(i,j)] = g;
			out[M_b(i,j)] = b;
		}
	}
}

/**
 * Main routine.
 */
int main(int argc, char **argv)
{
	int w, h, c;

	if (argc < 4)
	{
		panic("You must supply at least two input files and one output\n"
			"Usage: \n   %s <input file1> <input file2> ... <output.jpg>\n"
			"Last argument is always the output file!\n",
			argv[0]);
	}

	img_count = argc - 2;

	imgs = calloc(img_count, sizeof(uint8_t*));
	if (!imgs)
		panic("Cannot allocate image list\n");

	/* Load first image for width, height and channels reference. */
	imgs[0] = READ_IMAGE(argv[1], &WIDTH, &HEIGHT, &CHANNELS, 0);

	/* Load remaining images. */
	#pragma omp parallel for
	for (int i = 2; i < argc - 1; i++)
	{
		imgs[i - 1] = READ_IMAGE(argv[i], &w, &h, &c, 0);
		if (!imgs[i - 1])
			panic("Unable to read image (%s)\n",
				argv[i]);

		if (w != WIDTH || h != HEIGHT || c != CHANNELS)
		{
			panic("Width/Height/Channels (%d/%d/%d) for image (%s)\n"
				"differs from the reference image (%s) (%d/%d%d)\n",
				w,h,c,argv[i],argv[1],WIDTH,HEIGHT,CHANNELS);
		}
	}

	/* Output image. */
	out = malloc(WIDTH * HEIGHT * CHANNELS);
	if (!out)
		panic("Cannot allocate output file\n");

	do_image_median();

	WRITE_IMAGE(argv[argc-1], WIDTH, HEIGHT, CHANNELS, out, 92);

	/* Free everything =). */
	for (int i = 0; i < img_count; i++)
		FREE_IMAGE(imgs[i]);
	free(imgs);
	free(out);

	return (0);
}

Introduction

median.c is a small program that calculates the median/blending for multiple images and supports stb_image and libjpeg(-turbo) libraries.

Compile with -DUSE_STB_IMAGE to use stb_image and with -DUSE_LIBJPEG to use libjpeg(-turbo). In the absence of both,
stb_image will be used.

stb_image build instructions (default build):

gcc median.c -o median -O3 -fopenmp -lm

libjpeg(-turbo) build instructions:

gcc median.c -o median -O3 -fopenmp -ljpeg -DUSE_LIBJPEG

Common for both libraries:

If you do not want OpenMP and/or want to reduce dependencies just remove -fopenmp from your command-line options.

You may also want to reduce binary size by passing the flags: -fdata-sections -ffunction-sections -Wl,--gc-sections.
This reduces the file size from 123kB to 79kB or a 35.7% reduction for dynamic linking (w/ stb_image).

Dependencies:

stb_image

This file just requires a c99-compatible C compiler and 2 small header-only libraries: stb_image.h and stb_image_write.h.
Both can be obtained from: https://github.com/nothings/stb

or more specifically:
https://raw.githubusercontent.com/nothings/stb/master/stb_image.h
https://raw.githubusercontent.com/nothings/stb/master/stb_image_write.h
(versions used: v2.26 and v1.15, respectively, although I do not think that will have a huge break API in future versions).

stb_image supports other file formats besides JPEG, but for simplicity, I'm using only JPEG. To support other formats as well
(such as PNG, BMP...), comment the line #define STBI_ONLY_JPEG and/or define other macros. More info about this inside
the stb_image header.

or libjpeg(-turbo)

median.c requires stb_image or libjpeg, not both!. Please check the specific procedures of your Linux distribution for installing libjpeg.

Performance:

Eighteen JPEG images were tested with 4608x2592 (12MP). All builds were done with GCC v9.3 (w/ -O3) and run on a Core i5 7300HQ.
The following table summarizes the execution times and memory consumption for median.c, libvips (v8.10.6) and ImageMagick (v7.0.11-6).

(All tests were run 10 times, and the times represent the average runs. All input images were cached).

Programs	Memory Usage (MiB)	Time (seconds)	Speedup (with ImageMagick)
ImageMagick (v7.0.11-6)	2628 MiB	11.86s	1
libvips (v8.10.6)	173 MiB	3.73s	3.18x
median.c (w/ stb_image v2.26)	681 MiB	2,97s	3.98x
median.c (w/ libjpeg-turbo v2.0.5)	651 MiB	2.32s	5.11x

Overall, it can be noted that the memory consumption of median.c is significantly less (about ~4x) than ImageMagick, but
higher than libvips. As for execution times, median.c performs about 5.11 times better than ImageMagick and 1.60 times
better than libvips (using libjpeg-turbo).

It is worth noting that median.c is quite simple and perhaps should not be taken too seriously, but I am happy with the results
obtained =).