diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..ee7ce76 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "src/stb"] + path = src/stb + url = https://github.com/nothings/stb diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..fdddb29 --- /dev/null +++ b/LICENSE @@ -0,0 +1,24 @@ +This is free and unencumbered software released into the public domain. + +Anyone is free to copy, modify, publish, use, compile, sell, or +distribute this software, either in source code form or as a compiled +binary, for any purpose, commercial or non-commercial, and by any +means. + +In jurisdictions that recognize copyright laws, the author or authors +of this software dedicate any and all copyright interest in the +software to the public domain. We make this dedication for the benefit +of the public at large and to the detriment of our heirs and +successors. We intend this dedication to be an overt act of +relinquishment in perpetuity of all present and future rights to this +software under copyright law. + +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 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. + +For more information, please refer to diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..df4175b --- /dev/null +++ b/Makefile @@ -0,0 +1,12 @@ +PNAME = stbibalance +CFLAGS = -std=c99 -O2 -Wall -Wextra -Wno-unused-but-set-variable -Wno-unused-parameter -Werror +LDLIBS = -lm -s +SRCS = src/balance.c + +all: $(PNAME) + +$(PNAME): $(SRCS) + $(CC) $(CFLAGS) $^ $(LDLIBS) -o $@ + +clean: + rm -f $(PNAME) diff --git a/README.md b/README.md new file mode 100644 index 0000000..89a0027 --- /dev/null +++ b/README.md @@ -0,0 +1,51 @@ +# stbibalance + +Balance filter based of Gauss blur. +Used: +* Gauss blur: [iir_gauss_blur](https://github.com/arkanis/iir_gauss_blur). +* STB: [stb](https://github.com/nothings/stb). + +Balance is a thin filter (unlike others that use a thick difference between +the original image and its blurred version) that uses the difference between +a single and double radius blurred version of the original image +(single and double `sigma`). +This difference is used not directly, but by applying an `overlay` four times, +with the colors inverted. + +This filter was first applied in [STEX](https://github.com/ImageProcessing-ElectronicPublications/scantailor-experimental) (2024) +in a single-component version (Y, the values of the color components were aligned +in accordance with the brightness values before and after the filter). + +Here this filter is implemented in a full-color version. + +## Usage + +`./stbibalance [-h] [-s sigma] [-m mixed] input-file output-file` + +`-s sigma` The sigma of the gauss normal distribution (number >= 0.5). + Larger values result in a stronger blur. + +`-m mixed` The mixed coefficient. + +`-h` display this help and exit. + +You can use either sigma to specify the strengh of the blur. + +The performance is independent of the blur strengh (sigma). This tool is an +implementation of the paper "Recursive implementaion of the Gaussian filter" +by Ian T. Young and Lucas J. van Vliet. + +stb_image and stb_image_write by Sean Barrett and others is used to read and +write images. + +## Installation + +- Clone the repo or download the source `git clone --recurse-submodules https://github.com/ImageProcessing-ElectronicPublications/stbibalance` +- Execute `make` +- Done. Either use the `stbibalance` executable directly or copy it somewhere in your PATH. + +## Links + +* STB: [stb](https://github.com/nothings/stb). +* Gauss blur: [iir_gauss_blur](https://github.com/arkanis/iir_gauss_blur). +* Examples: [stbibalance-examples](https://github.com/ImageProcessing-ElectronicPublications/stbibalance-examples). diff --git a/src/balance.c b/src/balance.c new file mode 100644 index 0000000..6afd551 --- /dev/null +++ b/src/balance.c @@ -0,0 +1,226 @@ +#define _GNU_SOURCE +#include +#include + +#define STB_IMAGE_IMPLEMENTATION +#define STBI_FAILURE_USERMSG +#include "stb/stb_image.h" +#define STB_IMAGE_WRITE_IMPLEMENTATION +#include "stb/stb_image_write.h" +#define IIR_GAUSS_BLUR_IMPLEMENTATION +#include "iir_gauss_blur.h" + +void usage(char* progname) +{ + fprintf(stderr, + "%s %s %s\n", + "Usage:", progname, "[-h] [-s sigma] [-m mixed] input-file output.png\n" + "Balance filter an image and save it as PNG.\n" + ); +} + +void help(float sigma, float mix) +{ + fprintf(stderr, + "%s %f %s %f %s\n", + " -s sigma The sigma of the gauss normal distribution (number >= 0.5, default =", sigma, ").\n" + " Larger values result in a stronger blur.\n" + " -m mixed The mixed coefficient (number, default =", mix, ").\n" + " -h display this help and exit.\n" + "\n" + "You can use either sigma to specify the strengh of the blur.\n" + "\n" + "The performance is independent of the blur strengh (sigma). This tool is an\n" + "implementation of the paper \"Recursive implementaion of the Gaussian filter\"\n" + "by Ian T. Young and Lucas J. van Vliet.\n" + "\n" + "stb_image and stb_image_write by Sean Barrett and others is used to read and\n" + "write images.\n" + ); +} + +uint8_t* image_copy(unsigned int width, unsigned int height, unsigned char components, unsigned char* image) +{ + size_t image_size = height * width * components; + uint8_t* dest = (unsigned char*)malloc(image_size * sizeof(unsigned char)); + if (dest != NULL) + { + for (size_t i = 0; i < image_size; i++) + { + dest[i] = image[i]; + } + } + return dest; +} + +void image_overlay_blur(unsigned int width, unsigned int height, unsigned char components, unsigned char* blur1r, unsigned char* blur2r) +{ + size_t image_size = height * width * components; + if ((blur1r != NULL) && (blur2r != NULL)) + { + for (size_t i = 0; i < image_size; i++) + { + float b1r = blur1r[i]; + float b2r = blur2r[i]; + + /* overlay 1 */ + float base = 255.0f - b1r; + float overlay = b2r; + float retval1 = base; + if (base > 127.5f) + { + retval1 = 255.0f - retval1; + overlay = 255.0f - overlay; + } + retval1 *= overlay; + retval1 += retval1; + retval1 /= 255.0f; + if (base > 127.5f) + { + retval1 = 255.0f - retval1; + } + + /* overlay 2 */ + base = 255.0f - b2r; + overlay = b1r; + float retval2 = base; + if (base > 127.5f) + { + retval2 = 255.0f - retval2; + overlay = 255.0f - overlay; + } + retval2 *= overlay; + retval2 += retval2; + retval2 /= 255.0f; + if (base > 127.5f) + { + retval2 = 255.0f - retval2; + } + + blur1r[i] = (retval1 < 0.0f) ? 0 : ((retval1 < 255.0f) ? (uint8_t)(retval1 + 0.5f) : 255); + blur2r[i] = (retval2 < 0.0f) ? 0 : ((retval2 < 255.0f) ? (uint8_t)(retval2 + 0.5f) : 255); + } + } +} + +void image_overlay(unsigned int width, unsigned int height, unsigned char components, unsigned char* image, unsigned char* blur) +{ + size_t image_size = height * width * components; + if ((image != NULL) && (blur != NULL)) + { + for (size_t i = 0; i < image_size; i++) + { + float im = image[i]; + float bl = blur[i]; + + /* overlay*/ + float base = im; + float overlay = bl; + float retval = base; + if (base > 127.5f) + { + retval = 255.0f - retval; + overlay = 255.0f - overlay; + } + retval *= overlay; + retval += retval; + retval /= 255.0f; + if (base > 127.5f) + { + retval = 255.0f - retval; + } + + blur[i] = (retval < 0.0f) ? 0 : ((retval < 255.0f) ? (uint8_t)(retval + 0.5f) : 255); + } + } +} + +void image_mix(unsigned int width, unsigned int height, unsigned char components, float coef, unsigned char* image1, unsigned char* image2) +{ + size_t image_size = height * width * components; + if ((image1 != NULL) && (image2 != NULL)) + { + for (size_t i = 0; i < image_size; i++) + { + float im1 = image1[i]; + float im2 = image2[i]; + + /* overlay*/ + float retval = coef * im1 + (1.0 - coef) * im2; + + image2[i] = (retval < 0.0f) ? 0 : ((retval < 255.0f) ? (uint8_t)(retval + 0.5f) : 255); + } + } +} + +int main(int argc, char** argv) +{ + float sigma = 10.0f; + float mix = 1.0f; + + int opt; + while ( (opt = getopt(argc, argv, "s:m:h")) != -1 ) + { + switch(opt) + { + case 's': + sigma = strtof(optarg, NULL); + break; + case 'm': + mix = strtof(optarg, NULL); + break; + case 'h': + usage(argv[0]); + help(sigma, mix); + return 0; + default: + usage(argv[0]); + return 1; + } + } + + // Need at least two filenames after the last option + if (argc < optind + 2) + { + usage(argv[0]); + return 1; + } + + int width = 0, height = 0, components = 1; + uint8_t* image = stbi_load(argv[optind], &width, &height, &components, 0); + if (image == NULL) + { + fprintf(stderr, "Failed to load %s: %s.\n", argv[optind], stbi_failure_reason()); + return 2; + } + + uint8_t* blur1r = image_copy(width, height, components, image); + if (blur1r == NULL) + { + fprintf(stderr, "ERROR: not use memmory\n"); + return 3; + } + + uint8_t* blur2r = image_copy(width, height, components, image); + if (blur2r == NULL) + { + fprintf(stderr, "ERROR: not use memmory\n"); + return 3; + } + + iir_gauss_blur(width, height, components, blur1r, sigma); + iir_gauss_blur(width, height, components, blur2r, (sigma + sigma)); + + image_overlay_blur(width, height, components, blur1r, blur2r); + image_overlay(width, height, components, image, blur1r); + image_overlay(width, height, components, blur1r, blur2r); + image_mix(width, height, components, mix, blur2r, image); + + if ( stbi_write_png(argv[optind+1], width, height, components, image, 0) == 0 ) + { + fprintf(stderr, "Failed to save %s.\n", argv[optind+1]); + return 4; + } + + return 0; +} diff --git a/src/iir_gauss_blur.h b/src/iir_gauss_blur.h new file mode 100644 index 0000000..595a1ad --- /dev/null +++ b/src/iir_gauss_blur.h @@ -0,0 +1,214 @@ +/** + +IIR Gauss Filter v1.0 +By Stephan Soller +Based on the paper "Recursive implementation of the Gaussian filter" by Ian T. Young and Lucas J. van Vliet. +Licensed under the MIT license + +QUICK START + + #include ... + #include ... + #define IIR_GAUSS_BLUR_IMPLEMENTATION + #include "iir_gauss_blur.h" + ... + int width = 0, height = 0, components = 1; + uint8_t* image = stbi_load("foo.png", &width, &height, &components, 0); + float sigma = 10; + iir_gauss_blur(width, height, components, image, sigma); + stbi_write_png("foo.blurred.png", width, height, components, image, 0); + +This example uses stb_image.h to load the image, then blurs it and writes the result using stb_image_write.h. +`sigma` controls the strength of the blur. Higher values give you a blurrier image. + +DOCUMENTATION + +This is a single header file library. You'll have to define IIR_GAUSS_BLUR_IMPLEMENTATION before including this file to +get the implementation. Otherwise just the header will be included. + +The library only has a single function: iir_gauss_blur(width, height, components, image, sigma). + +- `width` and `height` are the dimensions of the image in pixels. +- `components` is the number of bytes per pixel. 1 for a grayscale image, 3 for RGB and 4 for RGBA. + The function can handle an arbitrary number of channels, so 2 or 7 will work as well. +- `image` is a pointer to the image data with `width * height` pixels, each pixel having `components` bytes (interleaved + 8-bit components). There is no padding between the scanlines of the image. + This is the format used by stb_image.h and stb_image_write.h and easy to work with. +- `sigma` is the strength of the blur. It's a number > 0.5 and most people seem to just eyeball it. + Start with e.g. a sigma of 5 and go up or down until you have the blurriness you want. + There are more informed ways to choose this parameter, see CHOOSING SIGMA below. + +The function mallocs an internal float buffer with the same dimensions as the image. If that turns out to be a +bottleneck fell free to move that out of the function. The source code is quite short and straight forward (even if the +math isn't). + +The function is an implementation of the paper "Recursive implementation of the Gaussian filter" by Ian T. Young and +Lucas J. van Vliet. It has nothing to do with recursive function calls, instead it's a special way to construct a +filter. Other (convolution based) gauss filters apply a kernel for each pixel and the kernel grows as sigma gets larger. +Meaning their performance degrades the more blurry you want your image to be. + +Instead The algorithm in the paper gets it done in just a few passes: A horizontal forward and backward pass and a +vertical forward and backward pass. The work done is independent of the blur radius and so you can have ridiculously +large blur radii without any performance impact. + +CHOOSING SIGMA + +There seem to be several rules of thumb out there to get a sigma for a given "blur radius". Usually this is something +like `radius = 2 * sigma`. So if you want to have a blur radius of 10px you can use `sigma = (1.0 / 2.0) * radius` to +get the sigma for it (5.0). I'm not sure what that "radius" is supposed to mean though. + +For my own projects I came up with two different kinds of blur radii and how to get a sigma for them: Given a big white +area on a black background, how far will the white "bleed out" into the surrounding black? How large is the distance +until the white (255) gets blurred down to something barely visible (smaller than 16) or even to nothing (smaller than +1)? There are to estimates to get the sigma for those radii: + + sigma = (1.0 / 1.42) * radius16; + sigma = (1.0 / 3.66) * radius1; + +Personally I use `radius16` to calculate the sigma when blurring normal images. Think: I want to blur a pixel across a +circle with the radius x so it's impact is barely visible at the edges. + +When I need to calculate padding I use `radius1`: When I have a black border of 100px around the image I can use a +`radius1` of 100 and be reasonable sure that I still got black at the edges. So given a `radius1` blur strength I can +use it as a padding width as well. + +I created those estimates by applying different sigmas (1 to 100) to a test image and measuring the effects with GIMP. +So take it with a grain of salt (or many). They're reasonable estimates but by no means exact. I tried to solve the +normal distribution to calculate the perfect sigma but gave up after a lot of confusion. If you know an exact solution +let me know. :) + +VERSION HISTORY + +v1.0 2018-08-30 Initial release + +**/ +#ifndef IIR_GAUSS_BLUR_HEADER +#define IIR_GAUSS_BLUR_HEADER +#ifdef __cplusplus + extern "C" { +#endif + +void iir_gauss_blur(unsigned int width, unsigned int height, unsigned char components, unsigned char* image, float sigma); + +#ifdef __cplusplus + } +#endif +#endif // IIR_GAUSS_BLUR_HEADER + +#ifdef IIR_GAUSS_BLUR_IMPLEMENTATION +#include +#include + +void iir_gauss_blur(unsigned int width, unsigned int height, unsigned char components, unsigned char* image, float sigma) { + // Create IDX macro but push any previous definition (and restore it later) so we don't overwrite a macro the user has possibly defined before us + #pragma push_macro("IDX") + #define IDX(x, y, n) ((y)*width*components + (x)*components + n) + + // Allocate buffers + float* buffer = (float*)malloc(width * height * components * sizeof(buffer[0])); + + // Calculate filter parameters for a specified sigma + // Use Equation 11b to determine q, do nothing if sigma is to small (should have no effect) or negative (doesn't make sense) + float q; + if (sigma >= 2.5) + q = 0.98711 * sigma - 0.96330; + else if (sigma >= 0.5) + q = 3.97156 - 4.14554 * sqrtf(1.0 - 0.26891 * sigma); + else + return; + + // Use equation 8c to determine b0, b1, b2 and b3 + float b0 = 1.57825 + 2.44413*q + 1.4281*q*q + 0.422205*q*q*q; + float b1 = 2.44413*q + 2.85619*q*q + 1.26661*q*q*q; + float b2 = -( 1.4281*q*q + 1.26661*q*q*q ); + float b3 = 0.422205*q*q*q; + // Use equation 10 to determine B + float B = 1.0 - (b1 + b2 + b3) / b0; + + // Horizontal forward pass (from paper: Implement the forward filter with equation 9a) + // The data is loaded from the byte image but stored in the float buffer + for(unsigned int y = 0; y < height; y++) { + float prev1[components], prev2[components], prev3[components]; + for(unsigned char n = 0; n < components; n++) { + prev1[n] = image[IDX(0, y, n)]; + prev2[n] = prev1[n]; + prev3[n] = prev2[n]; + } + + for(unsigned int x = 0; x < width; x++) { + for(unsigned char n = 0; n < components; n++) { + float val = B * image[IDX(x, y, n)] + (b1 * prev1[n] + b2 * prev2[n] + b3 * prev3[n]) / b0; + buffer[IDX(x, y, n)] = val; + prev3[n] = prev2[n]; + prev2[n] = prev1[n]; + prev1[n] = val; + } + } + } + + // Horizontal backward pass (from paper: Implement the backward filter with equation 9b) + for(unsigned int y = height-1; y < height; y--) { + float prev1[components], prev2[components], prev3[components]; + for(unsigned char n = 0; n < components; n++) { + prev1[n] = buffer[IDX(width-1, y, n)]; + prev2[n] = prev1[n]; + prev3[n] = prev2[n]; + } + + for(unsigned int x = width-1; x < width; x--) { + for(unsigned char n = 0; n < components; n++) { + float val = B * buffer[IDX(x, y, n)] + (b1 * prev1[n] + b2 * prev2[n] + b3 * prev3[n]) / b0; + buffer[IDX(x, y, n)] = val; + prev3[n] = prev2[n]; + prev2[n] = prev1[n]; + prev1[n] = val; + } + } + } + + // Vertical forward pass (from paper: Implement the forward filter with equation 9a) + for(unsigned int x = 0; x < width; x++) { + float prev1[components], prev2[components], prev3[components]; + for(unsigned char n = 0; n < components; n++) { + prev1[n] = buffer[IDX(x, 0, n)]; + prev2[n] = prev1[n]; + prev3[n] = prev2[n]; + } + + for(unsigned int y = 0; y < height; y++) { + for(unsigned char n = 0; n < components; n++) { + float val = B * buffer[IDX(x, y, n)] + (b1 * prev1[n] + b2 * prev2[n] + b3 * prev3[n]) / b0; + buffer[IDX(x, y, n)] = val; + prev3[n] = prev2[n]; + prev2[n] = prev1[n]; + prev1[n] = val; + } + } + } + + // Vertical backward pass (from paper: Implement the backward filter with equation 9b) + // Also write the result back into the byte image + for(unsigned int x = width-1; x < width; x--) { + float prev1[components], prev2[components], prev3[components]; + for(unsigned char n = 0; n < components; n++) { + prev1[n] = buffer[IDX(x, height-1, n)]; + prev2[n] = prev1[n]; + prev3[n] = prev2[n]; + } + + for(unsigned int y = height-1; y < height; y--) { + for(unsigned char n = 0; n < components; n++) { + float val = B * buffer[IDX(x, y, n)] + (b1 * prev1[n] + b2 * prev2[n] + b3 * prev3[n]) / b0; + image[IDX(x, y, n)] = val; + prev3[n] = prev2[n]; + prev2[n] = prev1[n]; + prev1[n] = val; + } + } + } + + // Free temporary buffers and restore any potential IDX macro + free(buffer); + #pragma pop_macro("IDX") +} +#endif // IIR_GAUSS_BLUR_IMPLEMENTATION diff --git a/src/stb b/src/stb new file mode 160000 index 0000000..5c20573 --- /dev/null +++ b/src/stb @@ -0,0 +1 @@ +Subproject commit 5c205738c191bcb0abc65c4febfa9bd25ff35234