[gegl] operations: fattal02 tone mapping operator
- From: Martin Nordholts <martinn src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gegl] operations: fattal02 tone mapping operator
- Date: Wed, 11 Aug 2010 18:03:31 +0000 (UTC)
commit cc6ab994c6816320b054af156c0a55489e654253
Author: Danny Robson <danny blubinc net>
Date: Wed Jun 9 18:31:08 2010 +1000
operations: fattal02 tone mapping operator
Implements the paper 'Gradient Domain High Dynamic Range Compression' by
Fattal et al. It is a local tone mapping operator, whose effect has
proven popular online
The implementation is a reasonably straight forward port from pfstmo
based code. There has been no attempt to improve/clean-up the base
implementation in this iteration.
operations/common/fattal02.c | 1323 +++++++++++++++++++++++++++++
tests/compositions/fattal02.xml | 17 +
tests/compositions/reference/fattal02.png | Bin 0 -> 200169 bytes
3 files changed, 1340 insertions(+), 0 deletions(-)
---
diff --git a/operations/common/fattal02.c b/operations/common/fattal02.c
new file mode 100644
index 0000000..d02f5bf
--- /dev/null
+++ b/operations/common/fattal02.c
@@ -0,0 +1,1323 @@
+/* This file is an image processing operation for GEGL
+ *
+ * GEGL is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 3 of the License, or (at your option) any later version.
+ *
+ * GEGL 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
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GEGL; if not, see <http://www.gnu.org/licenses/>.
+ *
+ * TMO:
+ * Copyright 2010 Danny Robson <danny blubinc net>
+ * (pfstmo) 2003-2004 Grzegorz Krawczyk <krawczyk mpi-sb mpg de>
+ *
+ * PDE:
+ * 2003-2004 Grzegorz Krawczyk <krawczyk mpi-sb mpg de>
+ * Rafal Mantiuk <mantiuk mpi-sb mpg de>
+ * Some code from Numerical Recipes in C
+ */
+
+#include "config.h"
+#include <glib/gi18n-lib.h>
+
+
+#ifdef GEGL_CHANT_PROPERTIES
+
+gegl_chant_double (alpha, _("Alpha"),
+ 0.0, 2.0, 1.0,
+ _("Gradient threshold for detail enhancement"))
+gegl_chant_double (beta, _("Beta"),
+ 0.1, 2.0, 0.9,
+ _("Strength of local detail enhancement"))
+gegl_chant_double (saturation, _("Saturation"),
+ 0.0, 1.0, 0.8,
+ _("Global colour saturation factor"))
+gegl_chant_double (noise, _("Noise"),
+ 0.0, 1.0, 0.0,
+ _("Gradient threshold for lowering detail enhancement"))
+
+
+#else
+
+#define GEGL_CHANT_TYPE_FILTER
+#define GEGL_CHANT_C_FILE "fattal02.c"
+
+#include "gegl-chant.h"
+#include "gegl-debug.h"
+#include <stdlib.h>
+
+static const gchar *OUTPUT_FORMAT = "RGB float";
+static const gint MINIMUM_PYRAMID = 32;
+
+/* I: pixel buffer, luminance with stride of 1
+ * R: rectangle, describes the buffer extent
+ * X: width coordinate
+ * Y: height coordinate
+ */
+#define _P(I,R,X,Y) ((I)[(Y) * (R)->width + (X)])
+
+/* The width/height of the pyramid at a level */
+#define LEVEL_WIDTH(extent, level) ((extent)->width / (1 << (level)))
+#define LEVEL_HEIGHT(extent, level) ((extent)->height / (1 << (level)))
+#define LEVEL_EXTENT(extent, level) \
+ ((GeglRectangle){ \
+ 0, 0, \
+ LEVEL_WIDTH ((extent), (level)), \
+ LEVEL_HEIGHT ((extent), (level)) \
+ })
+#define LEVEL_SIZE(extent, level) (LEVEL_EXTENT((extent), (level)).width * \
+ LEVEL_EXTENT((extent), (level)).height)
+
+#define MODYF 0 /* 1 or 0 (1 is better) */
+#define MINS 16 /* minimum size 4 6 or 100 */
+
+/* #define MODYF_SQRT -1.0f *//* -1 or 0 */
+#define SMOOTH_IT 1 /* minimum 1 */
+#define BCG_STEPS 20
+#define V_CYCLE 2 /* number of v-cycles */
+
+/* precision */
+#define EPS 1.0e-12
+
+static void
+linbcg (guint rows,
+ guint cols,
+ gfloat b[],
+ gfloat x[],
+ gint itol,
+ gfloat tol,
+ gint itmax,
+ gint *iter,
+ gfloat *err);
+
+
+/**
+ * Set all elements of the array to a give value.
+ *
+ * @param array array to modify
+ * @param value all elements of the array will be set to this value
+ */
+static inline void
+fattal02_set_array (gfloat *array,
+ guint size,
+ gfloat value)
+{
+ guint i;
+ for (i = 0; i < size; ++i)
+ array[i] = value;
+}
+
+
+static inline void
+fattal02_add_array (gfloat *accum,
+ guint size,
+ const gfloat *input)
+{
+ guint i;
+ for (i = 0; i < size; ++i)
+ accum[i] += input[i];
+}
+
+
+static inline void
+fattal02_copy_array (const gfloat *input,
+ gsize size,
+ gfloat *output)
+{
+ bcopy (input, output, size * sizeof (input[0]));
+}
+
+
+/*
+ * Full Multigrid Algorithm for solving partial differential equations
+ */
+
+static void
+fattal02_restrict (const gfloat *input,
+ const GeglRectangle *extent_i,
+ gfloat *output,
+ const GeglRectangle *extent_o)
+{
+ const guint inRows = extent_i->height,
+ inCols = extent_i->width;
+
+ const guint outRows = extent_o->height,
+ outCols = extent_o->width;
+
+ const gfloat dx = (gfloat)inCols / (gfloat)outCols,
+ dy = (gfloat)inRows / (gfloat)outRows;
+
+ const gfloat filterSize = 0.5;
+
+ gfloat sx, sy;
+ guint x, y;
+
+ for (y = 0, sy = dy / 2 - 0.5; y < outRows; ++y, sy += dy)
+ {
+ for (x = 0, sx = dx / 2 - 0.5; x < outCols; ++x, sx += dx )
+ {
+ gfloat pixVal = 0;
+ gfloat w = 0;
+ gint ix, iy;
+
+ for (ix = MAX (0, ceilf (sx - dx * filterSize));
+ ix <= MIN (floorf (sx + dx * filterSize), inCols - 1);
+ ++ix)
+ {
+ for (iy = MAX (0, ceilf (sy - dx * filterSize));
+ iy <= MIN (floorf (sy + dx * filterSize), inRows - 1);
+ ++iy)
+ {
+ pixVal += input[ix + iy * inCols];
+ w += 1;
+ }
+ }
+
+ output[x + y * outCols] = pixVal / w;
+ }
+ }
+}
+
+
+static void
+fattal02_prolongate (const gfloat *input,
+ const GeglRectangle *extent_i,
+ gfloat *output,
+ const GeglRectangle *extent_o)
+{
+ gfloat dx = (gfloat)extent_i->width / (gfloat)extent_o->width,
+ dy = (gfloat)extent_i->height / (gfloat)extent_o->height;
+
+ const guint outRows = extent_o->height,
+ outCols = extent_o->width;
+
+ const gfloat inRows = extent_i->height,
+ inCols = extent_i->width;
+
+ const float filterSize = 1;
+
+ gfloat sx, sy;
+ guint x, y;
+
+ for (y = 0, sy = -dy / 2; y < outRows; ++y, sy += dy)
+ {
+ for (x = 0, sx = -dx / 2; x < outCols; ++x, sx += dx )
+ {
+ gfloat pixVal = 0;
+ gfloat weight = 0;
+ gfloat ix, iy;
+
+ for (ix = MAX (0, ceilf (sx - filterSize));
+ ix <= MIN (floorf (sx + filterSize), inCols - 1);
+ ++ix)
+ {
+ for (iy = MAX (0, ceilf (sy - filterSize));
+ iy <= MIN (floorf (sy + filterSize), inRows - 1);
+ ++iy)
+ {
+ const gfloat fx = fabs (sx - ix),
+ fy = fabs (sy - iy),
+ fval = (1 - fx) * (1 - fy);
+
+ pixVal += input[(guint)ix + (guint)iy * (guint)inCols] * fval;
+ weight += fval;
+ }
+ }
+
+ g_return_if_fail (weight != 0);
+
+ output [x + y * outCols] = pixVal / weight;
+ }
+ }
+}
+
+
+static void
+fattal02_exact_solution (gfloat *F,
+ const GeglRectangle *extent_f,
+ gfloat *U,
+ const GeglRectangle *extent_u)
+{
+ /* pfstmo suggests that successive over-relaxation should be used here,
+ * followed by scaling by the square of the inverse of the sqrt of the array
+ * length. However it was commented out due to 'incorrect results', and the
+ * array zeroing was used in its place.
+ */
+ fattal02_set_array (U, extent_u->width * extent_u->height, 0.0f);
+ return;
+}
+
+
+/* smooth u using f at level */
+static void
+fattal02_smooth (gfloat *U,
+ const GeglRectangle *extent_u,
+ gfloat *F,
+ const GeglRectangle *extent_f)
+{
+ gint iter;
+ gfloat err;
+
+ linbcg (extent_u->height,
+ extent_u->width,
+ F, U, 1, 0.001,
+ BCG_STEPS, &iter, &err);
+
+ /* pfstmo notes here that 'gauss relaxation is too slow'. */
+}
+
+
+static void
+fattal02_calculate_defect (gfloat *D,
+ const GeglRectangle *extent_d,
+ gfloat *U,
+ const GeglRectangle *extent_u,
+ gfloat *F,
+ const GeglRectangle *extent_f)
+{
+ guint sx = extent_f->width,
+ sy = extent_f->height;
+ guint x, y;
+
+ for (y = 0; y < sy; ++y)
+ {
+ for (x = 0; x < sx; ++x)
+ {
+ guint w = (x == 0 ? 0 : x - 1),
+ n = (y == 0 ? 0 : y - 1),
+ s = (y + 1 == sy ? y : y + 1),
+ e = (x + 1 == sx ? x : x + 1);
+
+ _P (D, extent_d, x, y) = _P (F, extent_f, x, y) - (
+ _P (U, extent_u, e, y) +
+ _P (U, extent_u, w, y) +
+ _P (U, extent_u, x, n) +
+ _P (U, extent_u, x, s) -
+ 4.0 * _P (U, extent_u, x, y)
+ );
+ }
+ }
+}
+
+
+static void
+fattal02_solve_pde_multigrid (gfloat *F,
+ const GeglRectangle *extent_f,
+ gfloat *U,
+ const GeglRectangle *extent_u)
+{
+ guint xmax = extent_f->width,
+ ymax = extent_f->height;
+
+ gint i, /* index for simple loops */
+ k, /* index for iterating through levels */
+ k2; /* index for iterating through levels in V-cycles */
+
+ gint levels;
+
+ gfloat **RHS, /* given function f restricted on levels */
+ **IU, /* approximate initial sollutions on levels */
+ **VF; /* target functions in cycles (approximate sollution error (uh - ~uh) ) */
+
+ /* 1. restrict f to coarse-grid (by the way count the number of levels)
+ * k=0: fine-grid = f
+ * k=levels: coarsest-grid
+ */
+ {
+ guint mins = MIN (xmax, ymax);
+ levels = 0;
+
+ while (mins >= MINS)
+ {
+ levels++;
+ mins = mins / 2 + MODYF;
+ }
+ }
+
+ RHS = g_new (gfloat*, levels + 1);
+ IU = g_new (gfloat*, levels + 1);
+ VF = g_new (gfloat*, levels + 1);
+
+ RHS[0] = F;
+ VF[0] = g_new (gfloat, xmax * ymax);
+ IU[0] = g_new (gfloat, xmax * ymax);
+ fattal02_copy_array (U, xmax * ymax, IU[0]);
+
+ for (k = 0; k < levels; ++k)
+ {
+ RHS[k + 1] = g_new (gfloat, LEVEL_SIZE (extent_f, k + 1));
+ IU[k + 1] = g_new (gfloat, LEVEL_SIZE (extent_f, k + 1));
+ VF[k + 1] = g_new (gfloat, LEVEL_SIZE (extent_f, k + 1));
+
+ /* restrict from level k to level k+1 (coarser-grid) */
+ fattal02_restrict (RHS[k ], &LEVEL_EXTENT (extent_f, k ),
+ RHS[k + 1], &LEVEL_EXTENT (extent_f, k + 1));
+ }
+
+ /* 2. find exact solution at the coarsest-grid (k=levels) */
+ fattal02_exact_solution (RHS[levels], &LEVEL_EXTENT (extent_f, levels),
+ IU[levels], &LEVEL_EXTENT (extent_f, levels));
+
+ /* 3. nested iterations */
+ for (k = levels - 1; k >= 0; --k)
+ {
+ guint cycle;
+
+ /* 4. interpolate sollution from last coarse-grid to finer-grid
+ * interpolate from level k+1 to level k (finer-grid)
+ */
+ fattal02_prolongate (IU[k + 1], &LEVEL_EXTENT (extent_f, k + 1),
+ IU[k ], &LEVEL_EXTENT (extent_f, k ));
+
+ /* 4.1. first target function is the equation target function
+ * (following target functions are the defect)
+ */
+ fattal02_copy_array (RHS[k], LEVEL_SIZE (extent_f, k), VF[k]);
+
+ /* 5. V-cycle (twice repeated) */
+ for (cycle = 0; cycle < V_CYCLE; ++cycle)
+ {
+ /* 6. downward stroke of V */
+ for (k2 = k; k2 < levels; ++k2)
+ {
+ gfloat *D;
+
+ /* 7. pre-smoothing of initial sollution using target function
+ * zero for initial guess at smoothing
+ * (except for level k when iu contains prolongated result)
+ */
+ if (k2 != k)
+ {
+ fattal02_set_array (IU[k2], LEVEL_SIZE (extent_f, k2), 0.0f);
+ }
+
+ for (i = 0; i < SMOOTH_IT; ++i)
+ {
+ fattal02_smooth (IU[k2], &LEVEL_EXTENT (extent_f, k2),
+ VF[k2], &LEVEL_EXTENT (extent_f, k2));
+ }
+
+ /* 8. calculate defect at level
+ * d[k2] = Lh * ~u[k2] - f[k2]
+ */
+ D = g_new (gfloat, LEVEL_SIZE (extent_f, k2));
+ fattal02_calculate_defect ( D, &LEVEL_EXTENT (extent_f, k2),
+ IU[k2], &LEVEL_EXTENT (extent_f, k2),
+ VF[k2], &LEVEL_EXTENT (extent_f, k2));
+
+ /* 9. restrict deffect as target function for next coarser-grid
+ * def -> f[k2+1]
+ */
+ fattal02_restrict ( D, &LEVEL_EXTENT (extent_f, k2 ),
+ VF[k2 + 1], &LEVEL_EXTENT (extent_f, k2 + 1));
+ g_free (D);
+ }
+
+ /* 10. solve on coarsest-grid (target function is the deffect)
+ * iu[levels] should contain sollution for
+ * the f[levels] - last deffect, iu will now be the correction
+ */
+ fattal02_exact_solution (VF[levels], &LEVEL_EXTENT (extent_f, levels),
+ IU[levels], &LEVEL_EXTENT (extent_f, levels));
+
+ /* 11. upward stroke of V */
+ for (k2 = levels - 1; k2 >= k; --k2)
+ {
+ /* 12. interpolate correction from last coarser-grid to finer-grid
+ * iu[k2+1] -> cor
+ */
+ gfloat *C = g_new (gfloat, LEVEL_SIZE (extent_f, k2));
+ fattal02_prolongate (IU[k2 + 1], &LEVEL_EXTENT (extent_f, k2 + 1),
+ C, &LEVEL_EXTENT (extent_f, k2 ));
+
+ /* 13. add interpolated correction to initial sollution at level k2 */
+ fattal02_add_array (IU[k2], LEVEL_SIZE (extent_f, k2), C);
+ g_free (C);
+
+ /* 14. post-smoothing of current sollution using target function */
+ for (i = 0; i < SMOOTH_IT; ++i)
+ fattal02_smooth (IU[k2], &LEVEL_EXTENT (extent_f, k2),
+ VF[k2], &LEVEL_EXTENT (extent_f, k2));
+ }
+
+ } /*--- end of V-cycle */
+
+ } /*--- end of nested iteration */
+
+ /* 15. final sollution
+ * IU[0] contains the final sollution
+ */
+
+ fattal02_copy_array (IU[0], extent_f->width * extent_f->height, U);
+
+ g_free (VF[0]);
+ g_free (IU[0]);
+
+ for (k = 1; k <= levels; ++k)
+ {
+ g_free (RHS[k]);
+ g_free ( IU[k]);
+ g_free ( VF[k]);
+ }
+
+ g_free (RHS);
+ g_free ( IU);
+ g_free ( VF);
+}
+
+
+static void
+asolve (gulong n,
+ gfloat b[],
+ gfloat x[],
+ gint itrnsp)
+{
+ guint i;
+
+ for (i = 0; i < n; ++i)
+ x[i] = -4 * b[i];
+}
+
+static void
+atimes (guint rows,
+ guint cols,
+ gfloat x[],
+ gfloat res[],
+ gint itrnsp)
+{
+ guint r, c;
+
+#define IDX(R,C) ((R) * cols + (C))
+
+ for (r = 1; r < rows - 1; ++r)
+ {
+ for (c = 1; c < cols - 1; ++c)
+ {
+ res[IDX (r,c)] = x[IDX (r-1,c)] + x[IDX (r+1,c)] +
+ x[IDX (r,c-1)] + x[IDX (r,c+1)] - 4*x[IDX (r,c)];
+ }
+ }
+
+ for (r = 1; r < rows - 1; ++r)
+ {
+ res[IDX (r, 0)] = x[IDX (r - 1, 0)] +
+ x[IDX (r + 1, 0)] +
+ x[IDX (r , 1)] -
+ 3 * x[IDX (r , 0)];
+
+ res[IDX (r, cols - 1)] = x[IDX (r - 1, cols - 1)] +
+ x[IDX (r + 1, cols - 1)] +
+ x[IDX (r , cols - 2)] -
+ 3 * x[IDX (r , cols - 1)];
+ }
+
+ for (c = 1; c < cols - 1; ++c)
+ {
+ res[IDX (0, c)] = x[IDX (1, c )] +
+ x[IDX (0, c - 1)] +
+ x[IDX (0, c + 1)] -
+ 3 * x[IDX (0, c )];
+
+ res[IDX (rows - 1, c)] = x[IDX (rows - 2, c )] +
+ x[IDX (rows - 1, c - 1)] +
+ x[IDX (rows - 1, c + 1)] -
+ 3 * x[IDX (rows - 1, c )];
+ }
+
+ res[IDX (0 , 0)] = x[IDX (1 , 0)] +
+ x[IDX (0 , 1)] -
+ 2 * x[IDX (0 , 0)];
+ res[IDX (rows - 1, 0)] = x[IDX (rows - 2, 0)] +
+ x[IDX (rows - 1, 1)] -
+ 2 * x[IDX (rows - 1, 0)];
+ res[IDX (0 , cols - 1)] = x[IDX (1 , cols - 1)] +
+ x[IDX (0 , cols - 2)] -
+ 2 * x[IDX (0 , cols - 1)];
+ res[IDX (rows - 1, cols - 1)] = x[IDX (rows - 2, cols - 1)] +
+ x[IDX (rows - 1, cols - 2)] -
+ 2 * x[IDX (rows - 1, cols - 1)];
+}
+
+static gfloat
+snrm (gulong n,
+ gfloat sx[],
+ gint itol)
+{
+ gulong i;
+
+ if (itol <= 3)
+ {
+ gfloat ans = 0.0;
+ for (i = 0; i < n; ++i)
+ ans += sx[i] * sx[i];
+ return sqrtf (ans);
+ }
+ else
+ {
+ gulong isamax = 0;
+ for (i = 0; i < n; ++i)
+ if (fabs (sx[i]) > fabs (sx[isamax]))
+ isamax = i;
+ return fabs (sx[isamax]);
+ }
+}
+
+
+/**
+ * Biconjugate Gradient Method
+ * from Numerical Recipes in C
+ */
+static void
+linbcg (guint rows,
+ guint cols,
+ gfloat b[],
+ gfloat x[],
+ gint itol,
+ gfloat tol,
+ gint itmax,
+ gint *iter,
+ gfloat *err)
+{
+ guint n = rows * cols;
+
+ gulong j;
+ gfloat ak,akden,bk,bkden,bknum,bnrm,dxnrm,xnrm,zm1nrm,znrm;
+ gfloat *p,*pp,*r,*rr,*z,*zz;
+
+ p = g_new (gfloat, n);
+ pp = g_new (gfloat, n);
+ r = g_new (gfloat, n);
+ rr = g_new (gfloat, n);
+ z = g_new (gfloat, n);
+ zz = g_new (gfloat, n);
+
+ *iter=0;
+ atimes (rows, cols, x, r, 0);
+ for (j = 0; j < n; ++j)
+ {
+ r[j] = b[j] - r[j];
+ rr[j] = r[j];
+ }
+
+ atimes (rows, cols, r, rr, 0); /* minimum residual */
+ znrm = 1.0;
+
+ if (itol == 1)
+ {
+ bnrm = snrm (n, b, itol);
+ }
+ else if (itol == 2)
+ {
+ asolve (n, b, z, 0);
+ bnrm = snrm (n, z, itol);
+ }
+ else if (itol == 3 || itol == 4)
+ {
+ asolve (n, b, z, 0);
+ bnrm = snrm (n, z, itol);
+ asolve (n, r, z, 0);
+ znrm = snrm (n, z, itol);
+ }
+ else
+ {
+ g_warning ("illegal itol in linbcg");
+ }
+
+ asolve (n, r, z, 0);
+
+ while (*iter <= itmax)
+ {
+ ++(*iter);
+
+ zm1nrm = znrm;
+ asolve (n, rr, zz, 1);
+ for (bknum = 0.0, j = 0; j < n; ++j)
+ {
+ bknum += z[j] * rr[j];
+ }
+
+ if (*iter == 1)
+ {
+ for (j = 0; j < n; ++j)
+ {
+ p[j] = z[j];
+ pp[j] = zz[j];
+ }
+ }
+ else
+ {
+ bk = bknum / bkden;
+
+ for (j = 0; j < n; ++j)
+ {
+ p[j] = bk * p[j] + z[j];
+ pp[j] = bk * pp[j] + zz[j];
+ }
+ }
+
+ bkden = bknum;
+ atimes (rows, cols, p, z, 0);
+
+ for (akden = 0.0, j = 0; j < n; ++j)
+ {
+ akden += z[j] * pp[j];
+ }
+
+ ak = bknum / akden;
+ atimes (rows, cols, pp, zz, 1);
+
+ for (j = 0; j < n; ++j)
+ {
+ x[j] += ak * p[j];
+ r[j] -= ak * z[j];
+ rr[j] -= ak * zz[j];
+ }
+
+ asolve (n, r, z, 0);
+
+ if (itol == 1 || itol == 2)
+ {
+ znrm = 1.0;
+ *err = snrm (n, r, itol) / bnrm;
+ }
+ else if (itol == 3 || itol == 4)
+ {
+ znrm = snrm (n, z, itol);
+
+ if (fabs (zm1nrm - znrm) > EPS * znrm)
+ {
+ dxnrm = fabs (ak) * snrm (n, p, itol);
+ *err = znrm / fabs (zm1nrm - znrm) * dxnrm;
+ }
+ else
+ {
+ *err = znrm / bnrm;
+ continue;
+ }
+
+ xnrm = snrm (n, x, itol);
+ if (*err <= 0.5 * xnrm)
+ {
+ *err /= xnrm;
+ }
+ else
+ {
+ *err=znrm/bnrm;
+ continue;
+ }
+ }
+
+ if (*err <= tol)
+ break;
+ }
+
+ g_free (p);
+ g_free (pp);
+ g_free (r);
+ g_free (rr);
+ g_free (z);
+ g_free (zz);
+}
+
+
+/* Downscale the input buffer by a factor of two. Extent describes the input
+ * buffer. Assumes a pixel stride of 1, as we're really only dealing with
+ * luminance. Output should be preallocated with a size that is half of the
+ * input.
+ */
+static void
+fattal02_downsample (const gfloat *input,
+ const GeglRectangle *extent,
+ gfloat *output)
+{
+ guint x, y, width, height;
+ g_return_if_fail (input);
+ g_return_if_fail (extent);
+ g_return_if_fail (output);
+
+ width = extent->width / 2,
+ height = extent->height / 2;
+
+ g_return_if_fail (width > 0);
+ g_return_if_fail (height > 0);
+
+ for (y = 0; y < height; ++y)
+ {
+ for (x = 0; x < width; ++x)
+ {
+ gfloat p = 0.0f;
+
+ /* Sum the 4 pixels from the input which directly contribute to the
+ * output, and divide by four.
+ */
+ p += input[(2 * y + 0) * extent->width + (2 * x + 0)];
+ p += input[(2 * y + 0) * extent->width + (2 * x + 1)];
+ p += input[(2 * y + 1) * extent->width + (2 * x + 0)];
+ p += input[(2 * y + 1) * extent->width + (2 * x + 1)];
+
+ output [y * width + x] = p / 4.0f;
+ }
+ }
+}
+
+
+/* Blur the input buffer with a one pixel radius. Output should be
+ * preallocated with the same size as the input buffer. This must perform
+ * correctly when input and ouput alias.
+ */
+static void
+fattal02_gaussian_blur (const gfloat *input,
+ const GeglRectangle *extent,
+ gfloat *output)
+{
+ const guint width = extent->width,
+ height = extent->height,
+ size = width * height;
+ guint x, y;
+ gfloat *temp;
+
+ g_return_if_fail (input);
+ g_return_if_fail (extent);
+ g_return_if_fail (output);
+ g_return_if_fail (size > 0);
+
+ temp = g_new (gfloat, size);
+
+ /* horizontal blur */
+ for (y = 0; y < height; ++y)
+ {
+ for (x = 1; x < width - 1; ++x)
+ {
+ gfloat p = 2 * input[x + y * width];
+ p += input[x - 1 + y * width];
+ p += input[x + 1 + y * width];
+
+ temp[x + y * extent->width] = p / 4.0f;
+ }
+
+ temp[0 + y * width] = (3 * input[0 + y * width] +
+ input[1 + y * width]) / 4.0f;
+ temp[width - 1 + y * width] = (3 * input[width - 1 + y * width] +
+ input[width - 2 + y * width]) / 4.0f;
+ }
+
+ /* vertical blur */
+ for (x = 0; x < width; ++x)
+ {
+ for (y = 1; y < height - 1; ++y)
+ {
+ gfloat p = 2 * temp[x + y * width];
+ p += temp[x + (y - 1) * width];
+ p += temp[x + (y + 1) * width];
+
+ output[x + y * width] = p / 4.0f;
+ }
+
+ output[x + 0 * width] = (3 * output[x + 0 * width] +
+ output[x + 1 * width]) / 4.0f;
+ output[x + (height - 1) * width] = (3 * output[x + (height - 1) * width] +
+ output[x + (height - 2) * width]) / 4.0f;
+ }
+
+ g_free (temp);
+}
+
+
+static void
+fattal02_create_gaussian_pyramids (const gfloat *zero,
+ const GeglRectangle *extent,
+ gfloat **pyramid,
+ gint levels)
+{
+ gint i;
+ gfloat *blur;
+ GeglRectangle level_extent = *extent;
+
+ /* Copy the first level of the pyramid into place */
+ pyramid[0] = g_new (gfloat, level_extent.width * level_extent.height);
+ for (i = 0; i < level_extent.width * level_extent.height; ++i)
+ {
+ pyramid[0][i] = zero[i];
+ }
+
+ /* Establish a temporary blur buffer. The allocated memory will be used for
+ * progressively smaller levels, and we don't free this until the end.
+ */
+ blur = g_new (gfloat, level_extent.width * level_extent.height);
+ fattal02_gaussian_blur (pyramid[0], &level_extent, blur);
+
+ for (i = 1; i < levels; ++i)
+ {
+ level_extent.width /= 2;
+ level_extent.height /= 2;
+
+ g_return_if_fail (level_extent.width >= MINIMUM_PYRAMID);
+ g_return_if_fail (level_extent.height >= MINIMUM_PYRAMID);
+
+ /* Downsample the blurred buffer into the pyramid */
+ pyramid[i] = g_new (gfloat, LEVEL_SIZE (extent, i));
+ fattal02_downsample (blur, &LEVEL_EXTENT (extent, i - 1), pyramid[i]);
+
+ /* Blur the current level over the blur buffer */
+ fattal02_gaussian_blur (pyramid[i], &level_extent, blur);
+ }
+
+ g_free (blur);
+}
+
+
+/********************************************************************/
+
+static gfloat
+fattal02_calculate_gradients (const gfloat *input, /* H */
+ const GeglRectangle *extent, /* */
+ gfloat *output, /* G */
+ gint k)
+{
+ guint width = extent->width,
+ height = extent->height;
+ gfloat divider = powf (2.0f, k + 1),
+ average = 0.0f;
+
+ guint x, y;
+
+ for (y = 0; y < height; ++y)
+ {
+ for (x = 0; x < width; ++x)
+ {
+ gfloat gx, gy;
+ gint w, n, e, s;
+
+ w = (x == 0 ? 0 : x - 1);
+ n = (y == 0 ? 0 : y - 1);
+ s = (y + 1 == height ? y : y + 1);
+ e = (x + 1 == width ? x : x + 1);
+
+ gx = (input[w + y * width] - input[e + y * width]) / divider;
+ gy = (input[x + s * width] - input[x + n * width]) / divider;
+
+ output[x + y * width] = sqrtf (gx * gx + gy * gy);
+ average += output[x + y * width];
+ }
+ }
+
+ return average / (width * height);
+}
+
+
+/********************************************************************/
+
+static void
+fattal02_upsample (const gfloat *input,
+ const GeglRectangle *extent,
+ gfloat *output)
+{
+ guint width_i = extent->width,
+ height_i = extent->height,
+ width_o = width_i * 2,
+ height_o = height_i * 2;
+ guint x_o, y_o;
+
+ for (y_o = 0; y_o < height_o; ++y_o)
+ {
+ for (x_o = 0; x_o < width_o; ++x_o)
+ {
+ guint x_i = x_o / 2,
+ y_i = y_o / 2;
+
+ x_i = (x_i < width_i) ? x_i : width_i - 1;
+ y_i = (y_i < height_i) ? y_i : height_i - 1;
+
+ output[x_o + y_o * width_o] = input[x_i + y_i * width_i];
+ }
+ }
+}
+
+
+static void
+fattal02_FI_matrix (gfloat *FI,
+ const GeglRectangle *extent,
+ gfloat **gradients,
+ const gfloat *averages,
+ const gint levels,
+ const gfloat alfa,
+ const gfloat beta,
+ const gfloat noise)
+{
+ GeglRectangle level_extent = *extent;
+ gint i;
+ gfloat **fi;
+
+ level_extent.width = LEVEL_WIDTH (extent, levels - 1);
+ level_extent.height = LEVEL_HEIGHT (extent, levels - 1);
+
+ fi = g_new (gfloat*, levels);
+
+ fi[levels - 1] = g_new (gfloat, level_extent.width * level_extent.height);
+
+ for (i = 0; i < level_extent.width * level_extent.height; ++i)
+ {
+ fi[levels - 1][i] = 1.0f;
+ }
+
+ for (i = levels - 1; i >= 0; --i)
+ {
+ gint x, y;
+
+ level_extent.width = LEVEL_WIDTH (extent, i);
+ level_extent.height = LEVEL_HEIGHT (extent, i);
+
+ for (y = 0; y < level_extent.height; ++y)
+ for (x = 0; x < level_extent.width; ++x)
+ {
+ gfloat grad = gradients[i][x + y * level_extent.width],
+ a = alfa * averages[i],
+ value = 1.0f;
+
+ if (grad > 1e-4f)
+ value = a / (grad + noise) * powf ((grad + noise) / a, beta);
+ fi[i][x + y * level_extent.width] *= value;
+ }
+
+ /* create next level */
+ if (i > 1)
+ {
+ level_extent.width = LEVEL_WIDTH (extent, i - 1);
+ level_extent.height = LEVEL_HEIGHT (extent, i - 1);
+
+ fi[i - 1] = g_new (gfloat,
+ level_extent.width * level_extent.height);
+ }
+ else
+ {
+ fi[0] = FI; /* highest level -> result */
+ }
+
+ if (i > 0)
+ {
+ /* upsample to next level */
+ fattal02_upsample (fi[i ], &LEVEL_EXTENT (extent, i ), fi[i - 1]);
+ fattal02_gaussian_blur (fi[i - 1], &LEVEL_EXTENT (extent, i - 1), fi[i - 1]);
+ }
+ }
+
+ /* Careful not to delete the result memory in fi[0] */
+ for (i = 1; i < levels; ++i)
+ g_free (fi[i]);
+
+ g_free (fi);
+}
+
+/********************************************************************/
+
+
+static int
+fattal02_float_cmp (const void *_a,
+ const void *_b)
+{
+ const gfloat a = *(gfloat *)_a,
+ b = *(gfloat *)_b;
+
+ if (a < b) return -1;
+ if (a > b) return 1;
+ return 0;
+}
+
+
+static void
+fattal02_find_percentiles (const gfloat *array,
+ const guint size,
+ const gfloat min_percent,
+ gfloat *min_value,
+ const gfloat max_percent,
+ gfloat *max_value)
+{
+ guint i;
+ gfloat *sorting;
+
+ g_return_if_fail (min_percent <= max_percent);
+ g_return_if_fail (min_percent >= 0.0f && min_percent <= 1.0f);
+ g_return_if_fail (max_percent >= 0.0f && max_percent <= 1.0f);
+
+ sorting = g_new (gfloat, size);
+ for (i = 0; i < size; ++i)
+ {
+ sorting[i] = array[i];
+ }
+
+ qsort (sorting, size, sizeof (sorting[0]), fattal02_float_cmp);
+
+ *min_value = sorting[(guint)(min_percent * size)];
+ *max_value = sorting[(guint)(max_percent * size)];
+}
+
+/********************************************************************/
+
+static void
+fattal02_tonemap (const gfloat *input, /* Y */
+ const GeglRectangle *extent,
+ gfloat *output, /* L */
+ gfloat alfa,
+ gfloat beta,
+ gfloat noise)
+{
+ gint height = extent->height,
+ width = extent->width,
+ size = height * width;
+ gint x, y, i;
+ gfloat *H, *FI, *Gx, *Gy, *divergence, *U;
+ gint levels;
+ gfloat **pyramid;
+ gfloat **gradient,
+ *averages;
+
+ /* find max & min values, normalize to range 0..100 and take logarithm */
+ {
+ gfloat min_input = G_MAXFLOAT,
+ max_input = G_MINFLOAT;
+
+ for (i = 0; i < size; ++i)
+ {
+ min_input = MIN (min_input, input[i]);
+ max_input = MAX (max_input, input[i]);
+ }
+ g_return_if_fail (min_input <= max_input);
+
+ H = g_new (gfloat, size);
+ for (i = 0; i < size; ++i)
+ {
+ H[i] = log (100.0f * input[i] / max_input + 1e-4f);
+ }
+ }
+
+ GEGL_NOTE (GEGL_DEBUG_PROCESS, "calculating attenuation matrix");
+
+ /* create gaussian pyramids */
+ {
+ gint min_size = MIN (extent->width, extent->height);
+
+ for (levels = 0; min_size / 2 >= MINIMUM_PYRAMID; )
+ {
+ ++levels;
+ min_size /= 2;
+ }
+
+ pyramid = g_new (gfloat*, levels);
+ fattal02_create_gaussian_pyramids (H, extent, pyramid, levels);
+ }
+
+ /* calculate gradients and its average values on pyramid levels */
+ gradient = g_new (gfloat*, levels);
+ averages = g_new (gfloat, levels);
+
+ for (i = 0; i < levels; ++i)
+ {
+ gradient[i] = g_new (gfloat, LEVEL_SIZE (extent, i));
+ averages[i] = fattal02_calculate_gradients (pyramid[i],
+ &LEVEL_EXTENT (extent, i),
+ gradient[i],
+ i);
+ }
+
+ /* calculate fi matrix */
+ FI = g_new (gfloat, size);
+ fattal02_FI_matrix (FI, extent, gradient, averages, levels,
+ alfa, beta, noise);
+
+ /* attenuate gradients */
+ Gx = g_new (gfloat, size);
+ Gy = g_new (gfloat, size);
+
+ for (y = 0; y < extent->height; ++y)
+ {
+ for (x = 0; x < extent->width; ++x)
+ {
+ guint s = (y + 1 == height ? y : y + 1),
+ e = (x + 1 == width ? x : x + 1);
+
+ Gx[x + y * width] = ( H[e + y * width] - H[x + y * width]) *
+ FI[x + y * width];
+ Gy[x + y * width] = ( H[x + s * width] - H[x + y * width]) *
+ FI[x + y * width];
+ }
+ }
+
+ GEGL_NOTE (GEGL_DEBUG_PROCESS, "compressing gradients");
+
+ /* calculate divergence */
+ divergence = g_new (gfloat, size);
+ for (y = 0; y < height; ++y)
+ {
+ for (x = 0; x < width; ++x)
+ {
+ divergence[x + y * width] = Gx[x + y * width] + Gy[x + y * width];
+ if (x > 0) divergence[x + y * width] -= Gx[x - 1 + (y ) * width];
+ if (y > 0) divergence[x + y * width] -= Gy[x + (y - 1) * width];
+ }
+ }
+
+ GEGL_NOTE (GEGL_DEBUG_PROCESS, "recovering image");
+
+ /* solve pde and exponentiate (ie recover compressed image) */
+ U = g_new (gfloat, size);
+ fattal02_solve_pde_multigrid (divergence, extent, U, extent);
+
+ for (i = 0; i < size; ++i)
+ output[i] = expf (U[i]) - 1e-4f;
+
+ {
+ gfloat min, max, range;
+
+ /* remove percentile of min and max values and renormalize */
+ fattal02_find_percentiles (output, size,
+ 0.001f, &min,
+ 0.995f, &max);
+ range = max - min;
+
+ for (i = 0; i < size; ++i)
+ {
+ output[i] = (output[i] - min) / range;
+ if (output[i] <= 0.0f)
+ output[i] = 1e-4f;
+ }
+ }
+
+ /* clean up */
+ g_free (H);
+ for (i = 0; i < levels; ++i)
+ {
+ g_free ( pyramid[i]);
+ g_free (gradient[i]);
+ }
+
+ g_free (pyramid);
+ g_free (gradient);
+ g_free (averages);
+ g_free (FI);
+ g_free (Gx);
+ g_free (Gy);
+ g_free (divergence);
+ g_free (U);
+}
+
+
+static void
+fattal02_prepare (GeglOperation *operation)
+{
+ gegl_operation_set_format (operation, "input", babl_format (OUTPUT_FORMAT));
+ gegl_operation_set_format (operation, "output", babl_format (OUTPUT_FORMAT));
+}
+
+
+static GeglRectangle
+fattal02_get_required_for_output (GeglOperation *operation,
+ const gchar *input_pad,
+ const GeglRectangle *roi)
+{
+ GeglRectangle result = *gegl_operation_source_get_bounding_box (operation,
+ "input");
+ return result;
+}
+
+
+static GeglRectangle
+fattal02_get_cached_region (GeglOperation *operation,
+ const GeglRectangle *roi)
+{
+ return *gegl_operation_source_get_bounding_box (operation, "input");
+}
+
+
+static gboolean
+fattal02_process (GeglOperation *operation,
+ GeglBuffer *input,
+ GeglBuffer *output,
+ const GeglRectangle *result)
+{
+ const GeglChantO *o = GEGL_CHANT_PROPERTIES (operation);
+ gfloat noise;
+
+ const gint pix_stride = 3; /* RGBA */
+ gfloat *lum_in,
+ *lum_out,
+ *pix;
+ gint i;
+
+ g_return_val_if_fail (operation, FALSE);
+ g_return_val_if_fail (input, FALSE);
+ g_return_val_if_fail (output, FALSE);
+ g_return_val_if_fail (result, FALSE);
+
+ g_return_val_if_fail (babl_format_get_n_components (babl_format (OUTPUT_FORMAT)) == pix_stride, FALSE);
+
+ /* Adjust noise floor if not set by the user */
+ if (o->noise == 0.0)
+ {
+ noise = o->alpha * 0.1;
+ }
+ else
+ {
+ noise = o->noise;
+ }
+
+ /* Obtain the pixel data */
+ lum_in = g_new (gfloat, result->width * result->height);
+ lum_out = g_new (gfloat, result->width * result->height);
+
+ gegl_buffer_get (input, 1.0, result, babl_format ("Y float"),
+ lum_in, GEGL_AUTO_ROWSTRIDE);
+
+ pix = g_new (gfloat, result->width * result->height * pix_stride);
+ gegl_buffer_get (input, 1.0, result, babl_format (OUTPUT_FORMAT),
+ pix, GEGL_AUTO_ROWSTRIDE);
+
+ fattal02_tonemap (lum_in, result, lum_out, o->alpha, o->beta, noise);
+
+ for (i = 0; i < result->width * result->height * pix_stride; ++i)
+ {
+ pix[i] = (powf (pix[i] / lum_in[i / pix_stride],
+ o->saturation) *
+ lum_out[i / pix_stride]);
+ }
+
+ gegl_buffer_set (output, result, babl_format (OUTPUT_FORMAT), pix,
+ GEGL_AUTO_ROWSTRIDE);
+ return TRUE;
+}
+
+
+/*
+ */
+static void
+gegl_chant_class_init (GeglChantClass *klass)
+{
+ GeglOperationClass *operation_class;
+ GeglOperationFilterClass *filter_class;
+
+ operation_class = GEGL_OPERATION_CLASS (klass);
+ filter_class = GEGL_OPERATION_FILTER_CLASS (klass);
+
+ filter_class->process = fattal02_process;
+
+ operation_class->prepare = fattal02_prepare;
+ operation_class->get_required_for_output = fattal02_get_required_for_output;
+ operation_class->get_cached_region = fattal02_get_cached_region;
+
+ operation_class->name = "gegl:fattal02";
+ operation_class->categories = "tonemapping";
+ operation_class->description =
+ _("Adapt an image, which may have a high dynamic range, for "
+ "presentation using a low dynamic range. This operator attenuates "
+ "the magnitudes of local image gradients, producing luminance "
+ "within the range 0.0-1.0");
+}
+
+#endif
+
+
diff --git a/tests/compositions/fattal02.xml b/tests/compositions/fattal02.xml
new file mode 100644
index 0000000..de36999
--- /dev/null
+++ b/tests/compositions/fattal02.xml
@@ -0,0 +1,17 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<gegl>
+ <node operation='gegl:fattal02'>
+ <params>
+ <param name='alpha'>0.8</param>
+ <param name='beta'>0.6</param>
+ <param name='saturation'>0.8</param>
+ <param name='noise'>0.0</param>
+ </params>
+ </node>
+ <node operation='gegl:load'>
+ <params>
+ <param name='path'>data/parliament.hdr</param>
+ </params>
+ </node>
+</gegl>
+
diff --git a/tests/compositions/reference/fattal02.png b/tests/compositions/reference/fattal02.png
new file mode 100644
index 0000000..2dec16d
Binary files /dev/null and b/tests/compositions/reference/fattal02.png differ
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
