Implement Otsu's algorithm in autothresh0ld

src/filter/CMakeLists.txt

@@ -18,6 +18,7 @@ endif (${Cairo_FOUND})
 add_subdirectory (3dflippo)
 add_subdirectory (aech0r)
 add_subdirectory (alpha0ps)
+add_subdirectory (autothresh0ld)
 add_subdirectory (balanc0r)
 add_subdirectory (baltan)
 add_subdirectory (bluescreen0r)

src/filter/autothresh0ld/CMakeLists.txt

@@ -0,0 +1,11 @@
+set (SOURCES autothresh0ld.c)
+set (TARGET autothresh0ld)
+
+if (MSVC)
+  set (SOURCES ${SOURCES} ${FREI0R_DEF})
+endif (MSVC)
+
+add_library (${TARGET} MODULE ${SOURCES})
+
+# No «lib» prefix (name.so instead of libname.so)
+set_target_properties (${TARGET} PROPERTIES PREFIX "")

src/filter/autothresh0ld/autothresh0ld.c

@@ -0,0 +1,192 @@
+/**
+ * (c) Copyright 2025 Cynthia <cynthia2048@proton.me>
+ *
+ * This is based on Otsu's algorithm to segment an image into foreground
+ * or background based on the shape of the histogram. Instead of doing a
+ * hard threshold, we use the threshold obtained from Otsu's algorithm
+ * as the base for a sigmoidal transfer with a high slope; this produces
+ * a more eye-soothing threshold effect as seen in ImageMagick.
+ *
+ * This has the added benefit that, whereas the sigmoidal filter's base
+ * requires manual tuning, here it is determined algorithimically and thus
+ * can adapt on a frame-to-frame basis.
+ *
+ * This file is a Frei0r plugin.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ */
+
+#include <stdint.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "frei0r.h"
+#include "frei0r/math.h"
+
+#include "variance.h"
+
+typedef struct {
+    unsigned int width, height;
+    float slope;
+
+    uint8_t lut[256][256];
+    uint8_t *lumaframe;
+} s0ft0tsu_t;
+
+static void gen_sigmoid_lut (uint8_t lut[256][256], float slope)
+{
+    float k = expf(slope) / 255.0;
+
+    for (int j = 0; j < 256; ++j)
+        for (int i = 0; i < 256; ++i)
+            lut[j][i] = CLAMP (255.0 / (1.0 + expf(-k * (i - j))), 0, 255);
+}
+
+int f0r_init()
+{
+    return 0;
+}
+
+void f0r_deinit() {}
+
+f0r_instance_t f0r_construct(unsigned int width, unsigned int height)
+{
+    s0ft0tsu_t* s = calloc(1, sizeof(s0ft0tsu_t));
+
+    s->width = width;
+    s->height = height;
+    s->slope = 4.0;
+    s->lumaframe = malloc(s->width * s->height);
+
+    gen_sigmoid_lut(s->lut, s->slope);
+
+    return s;
+}
+
+void f0r_get_plugin_info(f0r_plugin_info_t* info)
+{
+    info->name = "autothreshold";
+    info->author = "Cynthia";
+    info->explanation = "Automatically threshold moving pictures";
+    info->major_version = 0;
+    info->minor_version = 1;
+    info->frei0r_version = FREI0R_MAJOR_VERSION;
+    info->plugin_type = F0R_PLUGIN_TYPE_FILTER;
+    info->color_model = F0R_COLOR_MODEL_RGBA8888;
+    info->num_params = 1;
+}
+
+void f0r_get_param_info(f0r_param_info_t *info, int index)
+{
+    switch (index)
+    {
+    case 0:
+        info->name = "Slope";
+        info->explanation = "Slope of sigmoidal transfer";
+        info->type = F0R_PARAM_DOUBLE;
+    break;
+    }
+}
+
+void f0r_get_param_value(f0r_instance_t inst, f0r_param_t param, int index)
+{
+    s0ft0tsu_t* s = inst;
+
+    switch (index)
+    {
+    case 0:
+        *(double*)param = s->slope / 7.0;
+    break;
+    }
+}
+
+void f0r_set_param_value(f0r_instance_t inst, f0r_param_t param, int index)
+{
+    s0ft0tsu_t* s = inst;
+
+    switch (index)
+    {
+    case 0:
+        s->slope = *(double*)param * 7.0;
+
+        gen_sigmoid_lut(s->lut, s->slope);
+    break;
+    }
+}
+
+void f0r_update(f0r_instance_t inst, double time,
+                const uint32_t* inframe, uint32_t* outframe)
+{
+    s0ft0tsu_t *s = inst;
+
+    uint8_t *src = (uint8_t*)inframe;
+    uint8_t *dst = s->lumaframe;
+    uint8_t r, g, b, luma;
+    size_t len = s->width * s->height;
+
+    // Normalised histogram (L = 256)
+    float hist[256] = {0};
+
+    for (int i = 0; i < len; ++i)
+    {
+        r = *src++;
+        g = *src++;
+        b = *src++;
+        src++; // Ignore alpha
+
+        luma = CLAMP(0.299 * r + 0.587 * g + 0.114 * b, 0, 255);
+        *dst++ = luma;
+
+        ++hist[luma]; // Add to histogram
+    }
+
+    // Normalise histogram; becomes a probability distribution
+    for (int i = 0; i < 256; ++i)
+        hist[i] /= len;
+
+    uint8_t thresh = 0.0;
+    float max_var = 0.0; // Maximum inter-class variance.
+
+    for (int i = 1; i < 256; ++i)
+    {
+        float var = find_variance(hist, i);
+        if (var > max_var)
+        {
+            thresh = i;
+            max_var = var;
+        }
+    }
+
+    src = (uint8_t*)s->lumaframe; // Replenish
+    dst = (uint8_t*)outframe;
+
+    for (int i = 0; i < len; ++i)
+    {
+        // Select the appropriate transfer
+        uint8_t* lut = s->lut[thresh];
+
+        luma = lut[*src++];
+
+        *dst++ = luma;
+        *dst++ = luma;
+        *dst++ = luma;
+        *dst++ = 0xFF; // TODO Copy alpha
+    }
+}
+
+void f0r_destruct(f0r_instance_t s)
+{
+    free(s);
+}

src/filter/autothresh0ld/variance.h

@@ -0,0 +1,84 @@
+#ifdef __SSE2__
+#include <emmintrin.h>
+#endif
+
+static float find_variance(float *hist, int thresh)
+{
+    float w_0 = 0.0, mu_0 = 0.0, w_1 = 0.0, mu_1 = 0.0;
+
+    for (int i = 0; i < thresh; ++i)
+    {
+        w_0 += hist[i];
+        mu_0 += i * hist[i];
+    }
+    for (int i = thresh; i < 256; ++i)
+    {
+        w_1 += hist[i];
+        mu_1 += i * hist[i];
+    }
+
+    float mu_diff = (mu_0/w_0 - mu_1/w_1);
+    return w_0*w_1*mu_diff*mu_diff;
+}
+
+#ifdef __SSE2__
+static float _sse1_hadd_ps(__m128 v)
+{
+    // Based on a StackOverflow answer.
+    __m128 shuf = _mm_shuffle_ps(v, v, _MM_SHUFFLE(2, 3, 0, 1));
+    __m128 sums = _mm_add_ps(v, shuf);
+    shuf        = _mm_movehl_ps(shuf, sums);
+    sums        = _mm_add_ss(sums, shuf);
+    return _mm_cvtss_f32(sums);
+}
+
+static float find_variance_sse2(float *hist, int thresh)
+{
+    __m128 vw_0 = _mm_setzero_ps(), vmu_0 = _mm_setzero_ps(),
+           vw_1 = _mm_setzero_ps(), vmu_1 = _mm_setzero_ps();
+    __m128 vcnt = _mm_set_ps(0, 1, 2, 3);
+
+    int thresh_low = (thresh / 4 + 0) * 4,
+        thresh_up  = (thresh / 4 + 1) * 4;
+
+    for (int i = 0; i < thresh_low; i+=4)
+    {
+        __m128 vhist = _mm_castsi128_ps(_mm_loadu_si128((void*)&hist[i]));
+
+        vw_0  = _mm_add_ps(vhist, vw_0);
+        vmu_0 = _mm_add_ps(_mm_mul_ps(vcnt, vhist), vw_0);
+        vcnt  = _mm_add_ps(vcnt, _mm_set1_ps(4));
+    }
+
+    // This is skipped, and handled as an edge case by non-SSE code.
+    vcnt = _mm_add_ps(vcnt, _mm_set1_ps(4));
+
+    for (int i = thresh_up; i < 256; i+=4)
+    {
+        __m128 vhist = _mm_castsi128_ps(_mm_loadu_si128((void*)&hist[i]));
+
+        vw_1  = _mm_add_ps(vhist, vw_1);
+        vmu_1 = _mm_add_ps(_mm_mul_ps(vcnt, vhist), vw_1);
+        vcnt  = _mm_add_ps(vcnt, _mm_set1_ps(4));
+    }
+
+    float w_0 = _sse1_hadd_ps(vw_0), mu_0 = _sse1_hadd_ps(vmu_0),
+          w_1 = _sse1_hadd_ps(vw_1), mu_1 = _sse1_hadd_ps(vmu_1);
+
+    // This edge case is here, because thresh may not be a multiple of 4.
+    for (int i = thresh_low; i < thresh; ++i)
+    {
+        w_0 += hist[i];
+        mu_0 += i * hist[i];
+    }
+
+    for (int i = thresh; i < thresh_up; ++i)
+    {
+        w_1 += hist[i];
+        mu_1 += i * hist[i];
+    }
+
+    float mu_diff = (mu_0/w_0 - mu_1/w_1);
+    return w_0*w_1*mu_diff*mu_diff;
+}
+#endif