[gegl] operations: mantiuk06 tone mapping operator
- From: Martin Nordholts <martinn src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gegl] operations: mantiuk06 tone mapping operator
- Date: Wed, 11 Aug 2010 18:24:12 +0000 (UTC)
commit 62b3e6ac185666b66ec8616447a84ed7d144096c
Author: Danny Robson <danny blubinc net>
Date: Wed Aug 11 20:31:29 2010 +0200
operations: mantiuk06 tone mapping operator
Implements the paper 'A Perceptual Framework for Contrast Processing of
High Dynamic Range Images'. It is a global tone mapping operator, which
constrains the set of contrasts across the full range of spatial
frequencies.
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/mantiuk06.c | 1643 ++++++++++++++++++++++++++++
tests/compositions/mantiuk06.xml | 15 +
tests/compositions/reference/mantiuk06.png | Bin 0 -> 193823 bytes
3 files changed, 1658 insertions(+), 0 deletions(-)
---
diff --git a/operations/common/mantiuk06.c b/operations/common/mantiuk06.c
new file mode 100644
index 0000000..d8d1c7f
--- /dev/null
+++ b/operations/common/mantiuk06.c
@@ -0,0 +1,1643 @@
+/* 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/>.
+ *
+ * Copyright 2010 Danny Robson <danny blubinc net>
+ * (pfstmo) 2007 Grzegorz Krawczyk <krawczyk mpi-sb mpg de>
+ * 2007-2008 Ed Brambley <E J Brambley damtp cam ac uk>
+ * Lebed Dmytry
+ * Radoslaw Mantiuk <radoslaw mantiuk gmail com>
+ * Rafal Mantiuk <mantiuk gmail com>
+ */
+
+#include "config.h"
+#include <glib/gi18n-lib.h>
+
+
+#ifdef GEGL_CHANT_PROPERTIES
+
+gegl_chant_double (contrast, _("Contrast"),
+ 0.0, 1.0, 0.1,
+ _("The amount of contrast compression"))
+gegl_chant_double (saturation, _("Saturation"),
+ 0.0, 2.0, 0.8,
+ _("Global colour saturation factor"))
+gegl_chant_double (detail, _("Detail"),
+ 1.0, 99.0, 1.0,
+ _("Level of emphasis on image gradient details"))
+
+
+#else
+
+#define GEGL_CHANT_TYPE_FILTER
+#define GEGL_CHANT_C_FILE "mantiuk06.c"
+
+#include "gegl-chant.h"
+#include <stdio.h>
+#include <stdlib.h>
+
+#ifdef HAVE_OPENMP
+#define _OMP(x) _Pragma(#x)
+#else
+#define _OMP(x) /* OMP disabled: "#x" */
+#endif
+
+/* Common return codes for operators */
+#define PFSTMO_OK 1 /* Successful */
+#define PFSTMO_ABORTED -1 /* User aborted (from callback) */
+#define PFSTMO_ERROR -2 /* Failed, encountered error */
+
+/* Return codes for the progress_callback */
+#define PFSTMO_CB_CONTINUE 1
+#define PFSTMO_CB_ABORT -1
+
+typedef struct pyramid_s {
+ guint rows;
+ guint cols;
+ gfloat *Gx;
+ gfloat *Gy;
+ struct pyramid_s *next;
+ struct pyramid_s *prev;
+} pyramid_t;
+
+
+#define PYRAMID_MIN_PIXELS 3
+#define LOOKUP_W_TO_R 107
+
+typedef int (*pfstmo_progress_callback)(int progress);
+
+
+static void mantiuk06_contrast_equalization (pyramid_t *pp,
+ const gfloat contrastFactor);
+static void mantiuk06_transform_to_luminance (pyramid_t *pp,
+ gfloat *const x,
+ pfstmo_progress_callback progress,
+ const gboolean bcg,
+ const int itmax,
+ const gfloat tol);
+static gfloat* mantiuk06_matrix_alloc (const guint size);
+static void mantiuk06_matrix_free (gfloat *m);
+static void mantiuk06_matrix_zero (const guint n,
+ gfloat *const m);
+static void mantiuk06_matrix_copy (const guint n,
+ const gfloat *const a,
+ gfloat *const b);
+static void mantiuk06_matrix_subtract (const guint n,
+ const gfloat *const a,
+ gfloat *const b);
+static void mantiuk06_matrix_multiply_const (const guint n,
+ gfloat *const a,
+ const gfloat val);
+static void mantiuk06_matrix_divide (const guint n,
+ const gfloat *const a,
+ gfloat *const b);
+static gfloat mantiuk06_matrix_dot_product (const guint n,
+ const gfloat *const a,
+ const gfloat *const b);
+static void mantiuk06_calculate_and_add_divergence (const int rows,
+ const int cols,
+ const gfloat *const Gx,
+ const gfloat *const Gy,
+ gfloat *const divG);
+static void mantiuk06_pyramid_calculate_divergence_sum (pyramid_t *pyramid,
+ gfloat *divG_sum);
+static void mantiuk06_calculate_scale_factor (const int n,
+ const gfloat *const G,
+ gfloat *const C);
+static void mantiuk06_pyramid_calculate_scale_factor (pyramid_t *pyramid,
+ pyramid_t *pC);
+static void mantiuk06_scale_gradient (const int n,
+ gfloat *const G,
+ const gfloat *const C);
+static void mantiuk06_pyramid_scale_gradient (pyramid_t *pyramid,
+ pyramid_t *pC);
+static void mantiuk06_pyramid_free (pyramid_t *pyramid);
+static pyramid_t* mantiuk06_pyramid_allocate (const int cols,
+ const int rows);
+static void mantiuk06_calculate_gradient (const int cols,
+ const int rows,
+ const gfloat *const lum,
+ gfloat *const Gx,
+ gfloat *const Gy);
+static void mantiuk06_pyramid_calculate_gradient (pyramid_t *pyramid,
+ gfloat *lum);
+static void mantiuk06_solveX (const gint n,
+ const gfloat *const b,
+ gfloat *const x);
+static void mantiuk06_multiplyA (pyramid_t *px,
+ pyramid_t *pyramid,
+ const gfloat *const x,
+ gfloat *const divG_sum);
+static void mantiuk06_linbcg (pyramid_t *pyramid,
+ pyramid_t *pC,
+ const gfloat *const b,
+ gfloat *const x,
+ const int itmax,
+ const gfloat tol,
+ pfstmo_progress_callback progress_cb);
+static void mantiuk06_lincg (pyramid_t *pyramid,
+ pyramid_t *pC,
+ const gfloat *const b,
+ gfloat *const x,
+ const int itmax,
+ const gfloat tol,
+ pfstmo_progress_callback progress_cb);
+static gfloat mantiuk06_lookup_table (const int n,
+ const gfloat *const in_tab,
+ const gfloat *const out_tab,
+ const gfloat val);
+static void mantiuk06_transform_to_R (const int n,
+ gfloat *const G);
+static void mantiuk06_pyramid_transform_to_R (pyramid_t *pyramid);
+static void mantiuk06_transform_to_G (const gint n,
+ gfloat *const R);
+static void mantiuk06_pyramid_transform_to_G (pyramid_t *pyramid);
+static void mantiuk06_pyramid_gradient_multiply (pyramid_t *pyramid,
+ const gfloat val);
+
+
+static const gchar *OUTPUT_FORMAT = "RGBA float";
+
+static gfloat W_table[] =
+{
+ 0.000000, 0.010000, 0.021180, 0.031830, 0.042628,
+ 0.053819, 0.065556, 0.077960, 0.091140, 0.105203,
+ 0.120255, 0.136410, 0.153788, 0.172518, 0.192739,
+ 0.214605, 0.238282, 0.263952, 0.291817, 0.322099,
+ 0.355040, 0.390911, 0.430009, 0.472663, 0.519238,
+ 0.570138, 0.625811, 0.686754, 0.753519, 0.826720,
+ 0.907041, 0.995242, 1.092169, 1.198767, 1.316090,
+ 1.445315, 1.587756, 1.744884, 1.918345, 2.109983,
+ 2.321863, 2.556306, 2.815914, 3.103613, 3.422694,
+ 3.776862, 4.170291, 4.607686, 5.094361, 5.636316,
+ 6.240338, 6.914106, 7.666321, 8.506849, 9.446889,
+ 10.499164, 11.678143, 13.000302, 14.484414, 16.151900,
+ 18.027221, 20.138345, 22.517282, 25.200713, 28.230715,
+ 31.655611, 35.530967, 39.920749, 44.898685, 50.549857,
+ 56.972578, 64.280589, 72.605654, 82.100619, 92.943020,
+ 105.339358, 119.530154, 135.795960, 154.464484, 175.919088,
+ 200.608905, 229.060934, 261.894494, 299.838552, 343.752526,
+ 394.651294, 453.735325, 522.427053, 602.414859, 695.706358,
+ 804.693100, 932.229271, 1081.727632, 1257.276717, 1463.784297,
+ 1707.153398, 1994.498731, 2334.413424, 2737.298517, 3215.770944,
+ 3785.169959, 4464.187290, 5275.653272, 6247.520102, 7414.094945,
+ 8817.590551, 10510.080619
+};
+
+static gfloat R_table[] =
+{
+ 0.000000, 0.009434, 0.018868, 0.028302, 0.037736,
+ 0.047170, 0.056604, 0.066038, 0.075472, 0.084906,
+ 0.094340, 0.103774, 0.113208, 0.122642, 0.132075,
+ 0.141509, 0.150943, 0.160377, 0.169811, 0.179245,
+ 0.188679, 0.198113, 0.207547, 0.216981, 0.226415,
+ 0.235849, 0.245283, 0.254717, 0.264151, 0.273585,
+ 0.283019, 0.292453, 0.301887, 0.311321, 0.320755,
+ 0.330189, 0.339623, 0.349057, 0.358491, 0.367925,
+ 0.377358, 0.386792, 0.396226, 0.405660, 0.415094,
+ 0.424528, 0.433962, 0.443396, 0.452830, 0.462264,
+ 0.471698, 0.481132, 0.490566, 0.500000, 0.509434,
+ 0.518868, 0.528302, 0.537736, 0.547170, 0.556604,
+ 0.566038, 0.575472, 0.584906, 0.594340, 0.603774,
+ 0.613208, 0.622642, 0.632075, 0.641509, 0.650943,
+ 0.660377, 0.669811, 0.679245, 0.688679, 0.698113,
+ 0.707547, 0.716981, 0.726415, 0.735849, 0.745283,
+ 0.754717, 0.764151, 0.773585, 0.783019, 0.792453,
+ 0.801887, 0.811321, 0.820755, 0.830189, 0.839623,
+ 0.849057, 0.858491, 0.867925, 0.877358, 0.886792,
+ 0.896226, 0.905660, 0.915094, 0.924528, 0.933962,
+ 0.943396, 0.952830, 0.962264, 0.971698, 0.981132,
+ 0.990566, 1.000000
+};
+
+
+/* upsample the matrix
+ * upsampled matrix is twice bigger in each direction than data[]
+ * res should be a pointer to allocated memory for bigger matrix
+ * cols and rows are the dimmensions of the output matrix
+ */
+static void
+mantiuk06_matrix_upsample (const gint outCols,
+ const gint outRows,
+ const gfloat *const in,
+ gfloat *const out)
+{
+ const int inRows = outRows/2;
+ const int inCols = outCols/2;
+ gint x, y;
+
+ /* Transpose of experimental downsampling matrix (theoretically the
+ * correct thing to do)
+ */
+
+ const gfloat dx = (gfloat)inCols / ((gfloat)outCols);
+ const gfloat dy = (gfloat)inRows / ((gfloat)outRows);
+ const gfloat factor = 1.0f / (dx*dy); /* This gives a genuine upsampling
+ * matrix, not the transpose of the
+ * downsampling matrix
+ */
+ /* const gfloat factor = 1.0f; */ /* Theoretically, this should be the
+ * best.
+ */
+
+ _OMP (omp parallel for schedule(static))
+ for (y = 0; y < outRows; y++)
+ {
+ const gfloat sy = y * dy;
+ const gint iy1 = ( y * inRows) / outRows;
+ const gint iy2 = MIN (((y+1) * inRows) / outRows, inRows-1);
+
+ for (x = 0; x < outCols; x++)
+ {
+ const gfloat sx = x * dx;
+ const gint ix1 = ( x * inCols) / outCols;
+ const gint ix2 = MIN (((x+1) * inCols) / outCols, inCols-1);
+
+ out[x + y*outCols] = (
+ ((ix1+1) - sx)*((iy1+1 - sy)) * in[ix1 + iy1*inCols] +
+ ((ix1+1) - sx)*(sy+dy - (iy1+1)) * in[ix1 + iy2*inCols] +
+ (sx+dx - (ix1+1))*((iy1+1 - sy)) * in[ix2 + iy1*inCols] +
+ (sx+dx - (ix1+1))*(sy+dx - (iy1+1)) * in[ix2 + iy2*inCols])*factor;
+ }
+ }
+}
+
+
+/* downsample the matrix */
+static void
+mantiuk06_matrix_downsample (const gint inCols,
+ const gint inRows,
+ const gfloat *const data,
+ gfloat *const res)
+{
+ const int outRows = inRows / 2;
+ const int outCols = inCols / 2;
+ gint x, y, i, j;
+
+ const gfloat dx = (gfloat)inCols / ((gfloat)outCols);
+ const gfloat dy = (gfloat)inRows / ((gfloat)outRows);
+
+ /* New downsampling by Ed Brambley:
+ * Experimental downsampling that assumes pixels are square and
+ * integrates over each new pixel to find the average value of the
+ * underlying pixels.
+ *
+ * Consider the original pixels laid out, and the new (larger)
+ * pixels layed out over the top of them. Then the new value for
+ * the larger pixels is just the integral over that pixel of what
+ * shows through; i.e., the values of the pixels underneath
+ * multiplied by how much of that pixel is showing.
+ *
+ * (ix1, iy1) is the coordinate of the top left visible pixel.
+ * (ix2, iy2) is the coordinate of the bottom right visible pixel.
+ * (fx1, fy1) is the fraction of the top left pixel showing.
+ * (fx2, fy2) is the fraction of the bottom right pixel showing.
+ */
+
+ const gfloat normalize = 1.0f/(dx*dy);
+ _OMP (omp parallel for schedule(static))
+ for (y = 0; y < outRows; y++)
+ {
+ const gint iy1 = ( y * inRows) / outRows;
+ const gint iy2 = ((y+1) * inRows) / outRows;
+ const gfloat fy1 = (iy1+1) - y * dy;
+ const gfloat fy2 = (y+1) * dy - iy2;
+
+ for (x = 0; x < outCols; x++)
+ {
+ const gint ix1 = ( x * inCols) / outCols;
+ const gint ix2 = ((x+1) * inCols) / outCols;
+ const gfloat fx1 = (ix1+1) - x * dx;
+ const gfloat fx2 = (x+1) * dx - ix2;
+
+ gfloat pixVal = 0.0f;
+ gfloat factorx, factory;
+
+ for (i = iy1; i <= iy2 && i < inRows; i++)
+ {
+ if (i == iy1)
+ {
+ /* We're just getting the bottom edge of this pixel */
+ factory = fy1;
+ }
+ else if (i == iy2)
+ {
+ /* We're just gettting the top edge of this pixel */
+ factory = fy2;
+ }
+ else
+ {
+ /* We've got the full height of this pixel */
+ factory = 1.0f;
+ }
+
+ for (j = ix1; j <= ix2 && j < inCols; j++)
+ {
+ if (j == ix1)
+ {
+ /* We've just got the right edge of this pixel */
+ factorx = fx1;
+ }
+ else if (j == ix2)
+ {
+ /* We've just got the left edge of this pixel */
+ factorx = fx2;
+ }
+ else
+ {
+ /* We've got the full width of this pixel */
+ factorx = 1.0f;
+ }
+
+ pixVal += data[j + i*inCols] * factorx * factory;
+ }
+ }
+
+ /* Normalize by the area of the new pixel */
+ res[x + y * outCols] = pixVal * normalize;
+ }
+ }
+}
+
+
+/* return = a - b */
+static inline void
+mantiuk06_matrix_subtract (const guint n,
+ const gfloat *const a,
+ gfloat *const b)
+{
+ guint i;
+
+ _OMP (omp parallel for schedule(static))
+ for (i = 0; i < n; i++)
+ b[i] = a[i] - b[i];
+}
+
+/* copy matix a to b, return = a */
+static inline void
+mantiuk06_matrix_copy (const guint n,
+ const gfloat *const a,
+ gfloat *const b)
+{
+ memcpy (b, a, sizeof (gfloat) * n);
+}
+
+/* multiply matrix a by scalar val */
+static inline void
+mantiuk06_matrix_multiply_const (const guint n,
+ gfloat *const a,
+ const gfloat val)
+{
+ guint i;
+ _OMP (omp parallel for schedule(static))
+ for (i = 0; i < n; i++)
+ a[i] *= val;
+}
+
+/* b = a[i] / b[i] */
+static inline void
+mantiuk06_matrix_divide (const guint n,
+ const gfloat *const a,
+ gfloat *const b)
+{
+ guint i;
+ _OMP (omp parallel for schedule(static))
+ for (i = 0; i < n; i++)
+ b[i] = a[i] / b[i];
+}
+
+
+static inline gfloat *
+mantiuk06_matrix_alloc (guint size)
+{
+ return g_new (gfloat, size);
+}
+
+
+static inline void
+mantiuk06_matrix_free (gfloat *m)
+{
+ g_return_if_fail (m);
+ g_free (m);
+}
+
+/* multiply vector by vector (each vector should have one dimension equal to 1) */
+static inline gfloat
+mantiuk06_matrix_dot_product (const guint n,
+ const gfloat *const a,
+ const gfloat *const b)
+{
+ gfloat val = 0;
+ guint j;
+
+ _OMP (omp parallel for reduction(+:val) schedule(static))
+ for (j = 0; j < n; j++)
+ val += a[j] * b[j];
+
+ return val;
+}
+
+/* set zeros for matrix elements */
+static inline void
+mantiuk06_matrix_zero (guint n,
+ gfloat *m)
+{
+ bzero (m, n * sizeof (gfloat));
+}
+
+/* calculate divergence of two gradient maps (Gx and Gy)
+ * divG(x,y) = Gx(x,y) - Gx(x-1,y) + Gy(x,y) - Gy(x,y-1)
+ */
+static inline void
+mantiuk06_calculate_and_add_divergence (const gint cols,
+ const gint rows,
+ const gfloat *const Gx,
+ const gfloat *const Gy,
+ gfloat *const divG)
+{
+ gint ky, kx;
+
+ _OMP (omp parallel for schedule(static))
+ for (ky = 0; ky < rows; ky++)
+ {
+ for (kx = 0; kx<cols; kx++)
+ {
+ gfloat divGx, divGy;
+ const int idx = kx + ky*cols;
+
+ if (kx == 0)
+ divGx = Gx[idx];
+ else
+ divGx = Gx[idx] - Gx[idx-1];
+
+ if (ky == 0)
+ divGy = Gy[idx];
+ else
+ divGy = Gy[idx] - Gy[idx - cols];
+
+ divG[idx] += divGx + divGy;
+ }
+ }
+}
+
+/* calculate the sum of divergences for the all pyramid level. the smaller
+ * divergence map is upsamled and added to the divergence map for the higher
+ * level of pyramid.
+ * temp is a temporary matrix of size (cols, rows), assumed to already be
+ * allocated
+ */
+static void
+mantiuk06_pyramid_calculate_divergence_sum (pyramid_t *pyramid,
+ gfloat *divG_sum)
+{
+ gfloat *temp = mantiuk06_matrix_alloc (pyramid->rows*pyramid->cols);
+
+ /* Find the coarsest pyramid, and the number of pyramid levels */
+ int levels = 1;
+ while (pyramid->next != NULL)
+ {
+ levels++;
+ pyramid = pyramid->next;
+ }
+
+ /* For every level, we swap temp and divG_sum. So, if there are an odd
+ * number of levels...
+ */
+ if (levels % 2)
+ {
+ gfloat *const dummy = divG_sum;
+ divG_sum = temp;
+ temp = dummy;
+ }
+
+ /* Add them all together */
+ while (pyramid != NULL)
+ {
+ gfloat *dummy;
+
+ /* Upsample or zero as needed */
+ if (pyramid->next != NULL)
+ {
+ mantiuk06_matrix_upsample (pyramid->cols,
+ pyramid->rows,
+ divG_sum,
+ temp);
+ }
+ else
+ {
+ mantiuk06_matrix_zero (pyramid->rows * pyramid->cols, temp);
+ }
+
+ /* Add in the (freshly calculated) divergences */
+ mantiuk06_calculate_and_add_divergence (pyramid->cols,
+ pyramid->rows,
+ pyramid->Gx,
+ pyramid->Gy,
+ temp);
+
+ /* Rather than copying, just switch round the pointers: we know we get
+ * them the right way round at the end.
+ */
+ dummy = divG_sum;
+ divG_sum = temp;
+ temp = dummy;
+
+ pyramid = pyramid->prev;
+ }
+
+ mantiuk06_matrix_free (temp);
+}
+
+/* calculate scale factors (Cx,Cy) for gradients (Gx,Gy)
+ * C is equal to EDGE_WEIGHT for gradients smaller than GFIXATE or
+ * 1.0 otherwise
+ */
+static inline void
+mantiuk06_calculate_scale_factor (const gint n,
+ const gfloat *const G,
+ gfloat *const C)
+{
+ const gfloat detectT = 0.001f;
+ const gfloat a = 0.038737;
+ const gfloat b = 0.537756;
+
+ gint i;
+
+ _OMP (omp parallel for schedule(static))
+ for (i = 0; i < n; i++)
+ {
+#if 1
+ const gfloat g = MAX (detectT, fabsf (G[i]));
+ C[i] = 1.0f / (a * powf (g, b));
+#else
+ const gfloat GFIXATE = 0.10f,
+ EDGE_WEIGHT = 0.01f;
+
+ if (fabsf (G[i]) < GFIXATE)
+ C[i] = 1.0f / EDGE_WEIGHT;
+ else
+ C[i] = 1.0f;
+#endif
+ }
+}
+
+/* calculate scale factor for the whole pyramid */
+static void
+mantiuk06_pyramid_calculate_scale_factor (pyramid_t *pyramid,
+ pyramid_t *pC)
+{
+ while (pyramid != NULL)
+ {
+ const int size = pyramid->rows * pyramid->cols;
+ mantiuk06_calculate_scale_factor (size, pyramid->Gx, pC->Gx);
+ mantiuk06_calculate_scale_factor (size, pyramid->Gy, pC->Gy);
+ pyramid = pyramid->next;
+ pC = pC->next;
+ }
+}
+
+/* Scale gradient (Gx and Gy) by C (Cx and Cy)
+ * G = G / C
+ */
+static inline void
+mantiuk06_scale_gradient (const gint n,
+ gfloat *const G,
+ const gfloat *const C)
+{
+ gint i;
+ _OMP (omp parallel for schedule(static))
+ for (i = 0; i < n; i++)
+ G[i] *= C[i];
+}
+
+/* scale gradients for the whole one pyramid with the use of (Cx,Cy) from the
+ * other pyramid
+ */
+static void
+mantiuk06_pyramid_scale_gradient (pyramid_t *pyramid,
+ pyramid_t *pC)
+{
+ while (pyramid != NULL)
+ {
+ const int size = pyramid->rows * pyramid->cols;
+ mantiuk06_scale_gradient (size, pyramid->Gx, pC->Gx);
+ mantiuk06_scale_gradient (size, pyramid->Gy, pC->Gy);
+ pyramid = pyramid->next;
+ pC = pC->next;
+ }
+}
+
+
+/* free memory allocated for the pyramid */
+static void
+mantiuk06_pyramid_free (pyramid_t *pyramid)
+{
+ while (pyramid)
+ {
+ pyramid_t *const next = pyramid->next;
+
+ if (pyramid->Gx != NULL)
+ {
+ g_free (pyramid->Gx);
+ pyramid->Gx = NULL;
+ }
+
+ if (pyramid->Gy != NULL)
+ {
+ g_free (pyramid->Gy);
+ pyramid->Gy = NULL;
+ }
+
+ pyramid->prev = NULL;
+ pyramid->next = NULL;
+ g_free (pyramid);
+ pyramid = next;
+ }
+}
+
+
+/* allocate memory for the pyramid */
+static pyramid_t*
+mantiuk06_pyramid_allocate (int cols,
+ int rows)
+{
+ pyramid_t *level = NULL;
+ pyramid_t *pyramid = NULL;
+ pyramid_t *prev = NULL;
+
+ while (rows >= PYRAMID_MIN_PIXELS && cols >= PYRAMID_MIN_PIXELS)
+ {
+ guint size;
+ level = g_new (pyramid_t, 1);
+ memset (level, 0, sizeof (pyramid_t));
+
+ level->rows = rows;
+ level->cols = cols;
+ size = level->rows * level->cols;
+ level->Gx = mantiuk06_matrix_alloc (size);
+ level->Gy = mantiuk06_matrix_alloc (size);
+
+ level->prev = prev;
+ if (prev != NULL)
+ prev->next = level;
+ prev = level;
+
+ if (pyramid == NULL)
+ pyramid = level;
+
+ rows /= 2;
+ cols /= 2;
+ }
+
+ return pyramid;
+}
+
+
+/* calculate gradients */
+static inline void
+mantiuk06_calculate_gradient (const gint cols,
+ const gint rows,
+ const gfloat *const lum,
+ gfloat *const Gx,
+ gfloat *const Gy)
+{
+ gint ky, kx;
+
+ _OMP (omp parallel for schedule(static))
+ for (ky = 0; ky < rows; ky++)
+ {
+ for (kx = 0; kx < cols; kx++)
+ {
+ const gint idx = kx + ky*cols;
+
+ if (kx == (cols - 1))
+ Gx[idx] = 0;
+ else
+ Gx[idx] = lum[idx+1] - lum[idx];
+
+ if (ky == (rows - 1))
+ Gy[idx] = 0;
+ else
+ Gy[idx] = lum[idx + cols] - lum[idx];
+ }
+ }
+}
+
+
+/* calculate gradients for the pyramid
+ * lum_temp gets overwritten!
+ */
+static void
+mantiuk06_pyramid_calculate_gradient (pyramid_t *pyramid,
+ gfloat *lum_temp)
+{
+ gfloat *temp = mantiuk06_matrix_alloc ((pyramid->rows / 2) *
+ (pyramid->cols / 2));
+ gfloat *const temp_saved = temp;
+
+ mantiuk06_calculate_gradient (pyramid->cols,
+ pyramid->rows,
+ lum_temp,
+ pyramid->Gx,
+ pyramid->Gy);
+
+ pyramid = pyramid->next;
+
+ while (pyramid)
+ {
+ gfloat *dummy;
+ mantiuk06_matrix_downsample (pyramid->prev->cols,
+ pyramid->prev->rows,
+ lum_temp,
+ temp);
+ mantiuk06_calculate_gradient (pyramid->cols,
+ pyramid->rows,
+ temp,
+ pyramid->Gx,
+ pyramid->Gy);
+
+ dummy = lum_temp;
+ lum_temp = temp;
+ temp = dummy;
+
+ pyramid = pyramid->next;
+ }
+
+ mantiuk06_matrix_free (temp_saved);
+}
+
+
+/* x = -0.25 * b */
+static inline void
+mantiuk06_solveX (const gint n,
+ const gfloat *const b,
+ gfloat *const x)
+{
+ gint i;
+
+ _OMP (omp parallel for schedule(static))
+ for (i = 0; i < n; i++)
+ x[i] = -0.25f * b[i];
+}
+
+/* divG_sum = A * x = sum (divG (x))
+ * memory for the temporary pyramid px should be allocated
+ */
+static inline void
+mantiuk06_multiplyA (pyramid_t *px,
+ pyramid_t *pC,
+ const gfloat *const x,
+ gfloat *const divG_sum)
+{
+ /* use divG_sum as a temp variable */
+ mantiuk06_matrix_copy (pC->rows*pC->cols, x, divG_sum);
+ mantiuk06_pyramid_calculate_gradient (px, divG_sum);
+ /* scale gradients by Cx,Cy from main pyramid */
+ mantiuk06_pyramid_scale_gradient (px, pC);
+ /* calculate the sum of divergences */
+ mantiuk06_pyramid_calculate_divergence_sum (px, divG_sum);
+}
+
+
+/* bi-conjugate linear equation solver
+ * overwrites pyramid!
+ */
+static void
+mantiuk06_linbcg (pyramid_t *pyramid,
+ pyramid_t *pC,
+ const gfloat *const b,
+ gfloat *const x,
+ const gint itmax,
+ const gfloat tol,
+ pfstmo_progress_callback progress_cb)
+{
+ const guint rows = pyramid->rows,
+ cols = pyramid->cols,
+ n = rows*cols;
+ const gfloat tol2 = tol*tol;
+
+ gfloat *const z = mantiuk06_matrix_alloc (n),
+ *const zz = mantiuk06_matrix_alloc (n),
+ *const p = mantiuk06_matrix_alloc (n),
+ *const pp = mantiuk06_matrix_alloc (n),
+ *const r = mantiuk06_matrix_alloc (n),
+ *const rr = mantiuk06_matrix_alloc (n),
+ *const x_save = mantiuk06_matrix_alloc (n);
+
+ const gfloat bnrm2 = mantiuk06_matrix_dot_product (n, b, b);
+
+ gfloat err2, bkden, saved_err2, ierr2, percent_sf;
+ gint iter = 0, num_backwards = 0, num_backwards_ceiling = 3;
+ gboolean reset = TRUE;
+
+ mantiuk06_multiplyA (pyramid, pC, x, r); /* r = A*x = divergence (x) */
+ mantiuk06_matrix_subtract (n, b, r); /* r = b - r */
+ err2 = mantiuk06_matrix_dot_product (n, r, r); /* err2 = r.r */
+
+ mantiuk06_multiplyA (pyramid, pC, r, rr); /* rr = A*r */
+
+ bkden = 0;
+ saved_err2 = err2;
+ mantiuk06_matrix_copy (n, x, x_save);
+
+ ierr2 = err2;
+ percent_sf = 100.0f /logf (tol2 * bnrm2 / ierr2);
+
+ for (; iter < itmax; iter++)
+ {
+ guint i;
+ gfloat bknum, ak, old_err2;
+
+ if (progress_cb != NULL)
+ progress_cb ((int) (logf (err2 / ierr2) * percent_sf));
+
+ mantiuk06_solveX (n, r, z); /* z = ~A (-1) * r = -0.25 * r */
+ mantiuk06_solveX (n, rr, zz); /* zz = ~A (-1) * rr = -0.25 * rr */
+
+ bknum = mantiuk06_matrix_dot_product (n, z, rr);
+
+ if (reset)
+ {
+ reset = FALSE;
+ mantiuk06_matrix_copy (n, z, p);
+ mantiuk06_matrix_copy (n, zz, pp);
+ }
+ else
+ {
+ const gfloat bk = bknum / bkden; /* beta = ... */
+
+ _OMP (omp parallel for schedule(static))
+ for (i = 0; i < n; i++)
+ {
+ p[i] = z[i] + bk * p[i];
+ pp[i] = zz[i] + bk * pp[i];
+ }
+ }
+
+ bkden = bknum; /* numerator becomes the dominator for the next iteration */
+
+ mantiuk06_multiplyA (pyramid, pC, p, z); /* z = A* p = divergence (p) */
+ mantiuk06_multiplyA (pyramid, pC, pp, zz); /* zz = A*pp = divergence (pp) */
+
+ ak = bknum / mantiuk06_matrix_dot_product (n, z, pp); /* alfa = ... */
+
+ _OMP (omp parallel for schedule(static))
+ for (i = 0 ; i < n ; i++ )
+ {
+ r[i] -= ak * z[i]; /* r = r - alfa * z */
+ rr[i] -= ak * zz[i]; /* rr = rr - alfa * zz */
+ }
+
+ old_err2 = err2;
+ err2 = mantiuk06_matrix_dot_product (n, r, r);
+
+ /* Have we gone unstable? */
+ if (err2 > old_err2)
+ {
+ /* Save where we've got to if it's the best yet */
+ if (num_backwards == 0 && old_err2 < saved_err2)
+ {
+ saved_err2 = old_err2;
+ mantiuk06_matrix_copy (n, x, x_save);
+ }
+
+ num_backwards++;
+ }
+ else
+ {
+ num_backwards = 0;
+ }
+
+ _OMP (omp parallel for schedule(static))
+ for (i = 0 ; i < n ; i++ )
+ x[i] += ak * p[i]; /* x = x + alfa * p */
+
+ if (num_backwards > num_backwards_ceiling)
+ {
+ /* Reset */
+ reset = TRUE;
+ num_backwards = 0;
+
+ /* Recover saved value */
+ mantiuk06_matrix_copy (n, x_save, x);
+
+ /* r = Ax */
+ mantiuk06_multiplyA (pyramid, pC, x, r);
+
+ /* r = b - r */
+ mantiuk06_matrix_subtract (n, b, r);
+
+ /* err2 = r.r */
+ err2 = mantiuk06_matrix_dot_product (n, r, r);
+ saved_err2 = err2;
+
+ /* rr = A*r */
+ mantiuk06_multiplyA (pyramid, pC, r, rr);
+ }
+
+ if (err2 / bnrm2 < tol2)
+ break;
+ }
+
+ /* Use the best version we found */
+ if (err2 > saved_err2)
+ {
+ err2 = saved_err2;
+ mantiuk06_matrix_copy (n, x_save, x);
+ }
+
+ if (err2/bnrm2 > tol2)
+ {
+ /* Not converged */
+ if (progress_cb != NULL)
+ progress_cb ((int) (logf (err2 / ierr2) * percent_sf));
+ if (iter == itmax)
+ g_warning ("mantiuk06: Warning: "
+ "Not converged (hit maximum iterations), "
+ "error = %g (should be below %g).",
+ sqrtf (err2 / bnrm2), tol);
+ else
+ g_warning ("mantiuk06: Warning: "
+ "Not converged (going unstable), "
+ "error = %g (should be below %g).",
+ sqrtf (err2 / bnrm2), tol);
+ }
+ else if (progress_cb != NULL)
+ progress_cb (100);
+
+ mantiuk06_matrix_free (x_save);
+ mantiuk06_matrix_free (p);
+ mantiuk06_matrix_free (pp);
+ mantiuk06_matrix_free (z);
+ mantiuk06_matrix_free (zz);
+ mantiuk06_matrix_free (r);
+ mantiuk06_matrix_free (rr);
+}
+
+
+/* conjugate linear equation solver
+ * overwrites pyramid!
+ */
+static void
+mantiuk06_lincg (pyramid_t *pyramid,
+ pyramid_t *pC,
+ const gfloat *const b,
+ gfloat *const x,
+ const int itmax,
+ const gfloat tol,
+ pfstmo_progress_callback progress_cb)
+{
+ const int rows = pyramid->rows,
+ cols = pyramid->cols,
+ n = rows*cols;
+ int iter = 0, num_backwards = 0, num_backwards_ceiling = 3;
+ const gfloat tol2 = tol*tol;
+
+ gfloat *const x_save = mantiuk06_matrix_alloc (n),
+ *const r = mantiuk06_matrix_alloc (n),
+ *const p = mantiuk06_matrix_alloc (n),
+ *const Ap = mantiuk06_matrix_alloc (n);
+
+ /* bnrm2 = ||b|| */
+ const gfloat bnrm2 = mantiuk06_matrix_dot_product (n, b, b);
+ gfloat rdotr, irdotr, saved_rdotr, percent_sf;
+
+ /* r = b - Ax */
+ mantiuk06_multiplyA (pyramid, pC, x, r);
+ mantiuk06_matrix_subtract (n, b, r);
+ rdotr = mantiuk06_matrix_dot_product (n, r, r); /* rdotr = r.r */
+
+ /* p = r */
+ mantiuk06_matrix_copy (n, r, p);
+
+ /* Setup initial vector */
+ saved_rdotr = rdotr;
+ mantiuk06_matrix_copy (n, x, x_save);
+
+ irdotr = rdotr;
+ percent_sf = 100.0f / logf (tol2 * bnrm2 / irdotr);
+ for (; iter < itmax; iter++)
+ {
+ gint i;
+ gfloat alpha, old_rdotr;
+
+ if (progress_cb != NULL) {
+ gint ret = progress_cb ((gint) (logf (rdotr / irdotr) * percent_sf));
+ if (ret == PFSTMO_CB_ABORT && iter > 0 ) /* User requested abort */
+ break;
+ }
+
+ /* Ap = A p */
+ mantiuk06_multiplyA (pyramid, pC, p, Ap);
+
+ /* alpha = r.r / (p . Ap) */
+ alpha = rdotr / mantiuk06_matrix_dot_product (n, p, Ap);
+
+ /* r = r - alpha Ap */
+ _OMP (omp parallel for schedule(static))
+ for (i = 0; i < n; i++)
+ r[i] -= alpha * Ap[i];
+
+ /* rdotr = r.r */
+ old_rdotr = rdotr;
+ rdotr = mantiuk06_matrix_dot_product (n, r, r);
+
+ /* Have we gone unstable? */
+ if (rdotr > old_rdotr)
+ {
+ /* Save where we've got to */
+ if (num_backwards == 0 && old_rdotr < saved_rdotr)
+ {
+ saved_rdotr = old_rdotr;
+ mantiuk06_matrix_copy (n, x, x_save);
+ }
+
+ num_backwards++;
+ }
+ else
+ {
+ num_backwards = 0;
+ }
+
+ /* x = x + alpha p */
+ _OMP (omp parallel for schedule(static))
+ for (i = 0; i < n; i++)
+ x[i] += alpha * p[i];
+
+
+ /* Exit if we're done */
+ if (rdotr/bnrm2 < tol2)
+ break;
+
+ if (num_backwards > num_backwards_ceiling)
+ {
+ /* Reset */
+ num_backwards = 0;
+ mantiuk06_matrix_copy (n, x_save, x);
+
+ /* r = Ax */
+ mantiuk06_multiplyA (pyramid, pC, x, r);
+
+ /* r = b - r */
+ mantiuk06_matrix_subtract (n, b, r);
+
+ /* rdotr = r.r */
+ rdotr = mantiuk06_matrix_dot_product (n, r, r);
+ saved_rdotr = rdotr;
+
+ /* p = r */
+ mantiuk06_matrix_copy (n, r, p);
+ }
+ else
+ {
+ /* p = r + beta p */
+ const gfloat beta = rdotr/old_rdotr;
+
+ _OMP (omp parallel for schedule(static))
+ for (i = 0; i < n; i++)
+ p[i] = r[i] + beta*p[i];
+ }
+ }
+
+ /* Use the best version we found */
+ if (rdotr > saved_rdotr)
+ {
+ rdotr = saved_rdotr;
+ mantiuk06_matrix_copy (n, x_save, x);
+ }
+
+ if (rdotr/bnrm2 > tol2)
+ {
+ /* Not converged */
+ if (progress_cb != NULL)
+ progress_cb ((int) (logf (rdotr / irdotr) * percent_sf));
+ if (iter == itmax)
+ g_warning ("mantiuk06: Warning: "
+ "Not converged (hit maximum iterations), "
+ "error = %g (should be below %g).",
+ sqrtf (rdotr/bnrm2), tol);
+ else
+ g_warning ("mantiuk06: Warning: "
+ "Not converged (going unstable), "
+ "error = %g (should be below %g).",
+ sqrtf (rdotr/bnrm2), tol);
+ }
+ else if (progress_cb != NULL)
+ progress_cb (100);
+
+ mantiuk06_matrix_free (x_save);
+ mantiuk06_matrix_free (p);
+ mantiuk06_matrix_free (Ap);
+ mantiuk06_matrix_free (r);
+}
+
+
+/* in_tab and out_tab should contain inccreasing float values */
+static inline gfloat
+mantiuk06_lookup_table (const gint n,
+ const gfloat *const in_tab,
+ const gfloat *const out_tab,
+ const gfloat val)
+{
+ gint j;
+
+ if (G_UNLIKELY (val < in_tab[0]))
+ return out_tab[0];
+
+ for (j = 1; j < n; j++)
+ if (val < in_tab[j])
+ {
+ const gfloat dd = (val - in_tab[j-1]) / (in_tab[j] - in_tab[j-1]);
+ return out_tab[j-1] + (out_tab[j] - out_tab[j-1]) * dd;
+ }
+
+ return out_tab[n-1];
+}
+
+
+/* transform gradient (Gx,Gy) to R */
+static inline void
+mantiuk06_transform_to_R (const gint n,
+ gfloat *const G)
+{
+ gint j;
+
+ _OMP (omp parallel for schedule(static))
+ for (j = 0; j < n; j++)
+ {
+ /* G to W */
+ const gfloat absG = fabsf (G[j]);
+ int sign;
+ if (G[j] < 0)
+ sign = -1;
+ else
+ sign = 1;
+ G[j] = (powf (10, absG) - 1.0f) * sign;
+
+ /* W to RESP */
+ if (G[j] < 0)
+ sign = -1;
+ else
+ sign = 1;
+
+ G[j] = sign * mantiuk06_lookup_table (LOOKUP_W_TO_R,
+ W_table,
+ R_table,
+ fabsf (G[j]));
+ }
+}
+
+/* transform gradient (Gx,Gy) to R for the whole pyramid */
+static inline void
+mantiuk06_pyramid_transform_to_R (pyramid_t *pyramid)
+{
+ while (pyramid != NULL)
+ {
+ const int size = pyramid->rows * pyramid->cols;
+ mantiuk06_transform_to_R (size, pyramid->Gx);
+ mantiuk06_transform_to_R (size, pyramid->Gy);
+ pyramid = pyramid->next;
+ }
+}
+
+/* transform from R to G */
+static inline void
+mantiuk06_transform_to_G (const gint n,
+ gfloat *const R)
+{
+ gint j;
+
+ _OMP (omp parallel for schedule(static))
+ for (j = 0; j < n; j++){
+ /* RESP to W */
+ gint sign;
+ if (R[j] < 0)
+ sign = -1;
+ else
+ sign = 1;
+ R[j] = sign * mantiuk06_lookup_table (LOOKUP_W_TO_R,
+ R_table,
+ W_table,
+ fabsf (R[j]));
+
+ /* W to G */
+ if (R[j] < 0)
+ sign = -1;
+ else
+ sign = 1;
+ R[j] = log10f (fabsf (R[j]) + 1.0f) * sign;
+
+ }
+}
+
+/* transform from R to G for the pyramid */
+static inline void
+mantiuk06_pyramid_transform_to_G (pyramid_t *pyramid)
+{
+ while (pyramid != NULL)
+ {
+ mantiuk06_transform_to_G (pyramid->rows*pyramid->cols, pyramid->Gx);
+ mantiuk06_transform_to_G (pyramid->rows*pyramid->cols, pyramid->Gy);
+ pyramid = pyramid->next;
+ }
+}
+
+/* multiply gradient (Gx,Gy) values by float scalar value for the
+ * whole pyramid
+ */
+static inline void
+mantiuk06_pyramid_gradient_multiply (pyramid_t *pyramid,
+ const gfloat val)
+{
+ while (pyramid != NULL)
+ {
+ mantiuk06_matrix_multiply_const (pyramid->rows * pyramid->cols,
+ pyramid->Gx, val);
+ mantiuk06_matrix_multiply_const (pyramid->rows * pyramid->cols,
+ pyramid->Gy, val);
+ pyramid = pyramid->next;
+ }
+}
+
+
+static int
+mantiuk06_sort_float (const void *const v1,
+ const void *const v2)
+{
+ if (*((gfloat*)v1) < *((gfloat*)v2))
+ return -1;
+
+ if (G_LIKELY (*((gfloat*)v1) > *((gfloat*)v2)))
+ return 1;
+
+ return 0;
+}
+
+
+/* transform gradients to luminance */
+static void
+mantiuk06_transform_to_luminance (pyramid_t *pp,
+ gfloat *const x,
+ pfstmo_progress_callback progress,
+ const gboolean bcg,
+ const gint itmax,
+ const gfloat tol)
+{
+ gfloat *b;
+ pyramid_t *pC = mantiuk06_pyramid_allocate (pp->cols, pp->rows);
+ /* calculate (Cx,Cy) */
+ mantiuk06_pyramid_calculate_scale_factor (pp, pC);
+ /* scale small gradients by (Cx,Cy); */
+ mantiuk06_pyramid_scale_gradient (pp, pC);
+
+ b = mantiuk06_matrix_alloc (pp->cols * pp->rows);
+ /* calculate the sum of divergences (equal to b) */
+ mantiuk06_pyramid_calculate_divergence_sum (pp, b);
+
+ /* calculate luminances from gradients */
+ if (bcg)
+ mantiuk06_linbcg (pp, pC, b, x, itmax, tol, progress);
+ else
+ mantiuk06_lincg (pp, pC, b, x, itmax, tol, progress);
+
+ mantiuk06_matrix_free (b);
+ mantiuk06_pyramid_free (pC);
+}
+
+
+struct hist_data
+{
+ gfloat size;
+ gfloat cdf;
+ gint index;
+};
+
+static int
+mantiuk06_hist_data_order (const void *const v1,
+ const void *const v2)
+{
+ if (((struct hist_data*) v1)->size < ((struct hist_data*) v2)->size)
+ return -1;
+
+ if (((struct hist_data*) v1)->size > ((struct hist_data*) v2)->size)
+ return 1;
+
+ return 0;
+}
+
+
+static int
+mantiuk06_hist_data_index (const void *const v1,
+ const void *const v2)
+{
+ return ((struct hist_data*) v1)->index - ((struct hist_data*) v2)->index;
+}
+
+
+static void
+mantiuk06_contrast_equalization (pyramid_t *pp,
+ const gfloat contrastFactor )
+{
+ gint i, idx;
+ struct hist_data *hist;
+ gint total_pixels = 0;
+
+ /* Count sizes */
+ pyramid_t *l = pp;
+ while (l != NULL)
+ {
+ total_pixels += l->rows * l->cols;
+ l = l->next;
+ }
+
+ /* Allocate memory */
+ hist = g_new (struct hist_data, total_pixels);
+
+ /* Build histogram info */
+ l = pp;
+ idx = 0;
+ while (l != NULL)
+ {
+ const int pixels = l->rows*l->cols;
+ const int offset = idx;
+ gint c;
+
+ _OMP (omp parallel for schedule(static))
+ for (c = 0; c < pixels; c++)
+ {
+ hist[c+offset].size = sqrtf (l->Gx[c] * l->Gx[c] +
+ l->Gy[c] * l->Gy[c]);
+ hist[c+offset].index = c + offset;
+ }
+ idx += pixels;
+ l = l->next;
+ }
+
+ /* Generate histogram */
+ qsort (hist, total_pixels, sizeof (struct hist_data),
+ mantiuk06_hist_data_order);
+
+ /* Calculate cdf */
+ {
+ const gfloat norm = 1.0f / (gfloat) total_pixels;
+ _OMP (omp parallel for schedule(static))
+ for (i = 0; i < total_pixels; i++)
+ hist[i].cdf = ((gfloat) i) * norm;
+ }
+
+ /* Recalculate in terms of indexes */
+ qsort (hist, total_pixels,
+ sizeof (struct hist_data),
+ mantiuk06_hist_data_index);
+
+ /*Remap gradient magnitudes */
+ l = pp;
+ idx = 0;
+ while (l != NULL )
+ {
+ gint c;
+ const int pixels = l->rows*l->cols;
+ const int offset = idx;
+
+ _OMP (omp parallel for schedule(static))
+ for (c = 0; c < pixels; c++)
+ {
+ const gfloat scale = contrastFactor *
+ hist[c+offset].cdf /
+ hist[c+offset].size;
+ l->Gx[c] *= scale;
+ l->Gy[c] *= scale;
+ }
+ idx += pixels;
+ l = l->next;
+ }
+
+ g_free (hist);
+}
+
+
+/* tone mapping */
+static int
+mantiuk06_contmap (const int c,
+ const int r,
+ gfloat *const rgb,
+ gfloat *const Y,
+ const gfloat contrastFactor,
+ const gfloat saturationFactor,
+ const gboolean bcg,
+ const int itmax,
+ const gfloat tol,
+ pfstmo_progress_callback progress)
+{
+ const guint n = c*r;
+ guint j;
+
+ /* Normalize */
+ gfloat Ymax = Y[0],
+ clip_min;
+
+ for (j = 1; j < n; j++)
+ Ymax = MAX (Y[j], Ymax);
+
+ clip_min = 1e-7f * Ymax;
+ _OMP (omp parallel for schedule(static))
+ for (j = 0; j < n * 4; j++)
+ if (G_UNLIKELY (rgb[j] < clip_min)) rgb[j] = clip_min;
+
+ _OMP (omp parallel for schedule(static))
+ for (j = 0; j < n; j++)
+ if (G_UNLIKELY ( Y[j] < clip_min)) Y[j] = clip_min;
+
+ _OMP (omp parallel for schedule(static))
+ for (j = 0; j < n; j++)
+ {
+ rgb[j * 4 + 0] /= Y[j];
+ rgb[j * 4 + 1] /= Y[j];
+ rgb[j * 4 + 2] /= Y[j];
+ Y[j] = log10f (Y[j]);
+ }
+
+ {
+ /* create pyramid */
+ pyramid_t *pp = mantiuk06_pyramid_allocate (c,r);
+ gfloat *tY = mantiuk06_matrix_alloc (n);
+
+ /* copy Y to tY */
+ mantiuk06_matrix_copy (n, Y, tY);
+ /* calculate gradients for pyramid, destroys tY */
+ mantiuk06_pyramid_calculate_gradient (pp,tY);
+ mantiuk06_matrix_free (tY);
+ /* transform gradients to R */
+ mantiuk06_pyramid_transform_to_R (pp);
+
+ /* Contrast map */
+ if (contrastFactor > 0.0f)
+ /* Contrast mapping */
+ mantiuk06_pyramid_gradient_multiply (pp, contrastFactor);
+ else
+ /* Contrast equalization */
+ mantiuk06_contrast_equalization (pp, -contrastFactor);
+
+ /* transform R to gradients */
+ mantiuk06_pyramid_transform_to_G (pp);
+ /* transform gradients to luminance Y */
+ mantiuk06_transform_to_luminance (pp, Y, progress, bcg, itmax, tol);
+ mantiuk06_pyramid_free (pp);
+ }
+
+ /* Renormalize luminance */
+ {
+ const gfloat CUT_MARGIN = 0.1f;
+ gfloat *temp = mantiuk06_matrix_alloc (n),
+ trim, delta, l_min, l_max;
+
+ /* copy Y to temp */
+ mantiuk06_matrix_copy (n, Y, temp);
+ /* sort temp in ascending order */
+ qsort (temp, n, sizeof (gfloat), mantiuk06_sort_float);
+
+ /* const float median = (temp[(int)((n-1)/2)] + temp[(int)((n-1)/2+1)]) * 0.5f; */
+ /* calculate median */
+ trim = (n - 1) * CUT_MARGIN * 0.01f;
+ delta = trim - floorf (trim);
+ l_min = temp[(int)floorf (trim)] * delta +
+ temp[(int) ceilf (trim)] * (1.0f - delta);
+
+ trim = (n - 1) * (100.0f - CUT_MARGIN) * 0.01f;
+ delta = trim - floorf (trim);
+ l_max = temp[(int)floorf (trim)] * delta +
+ temp[(int) ceilf (trim)] * (1.0f - delta);
+
+ mantiuk06_matrix_free (temp);
+ {
+ const gfloat disp_dyn_range = 2.3f;
+ _OMP (omp parallel for schedule(static))
+ for (j = 0; j < n; j++)
+ /* x scaled */
+ Y[j] = ( Y[j] - l_min) /
+ (l_max - l_min) *
+ disp_dyn_range - disp_dyn_range;
+ }
+
+ /* Transform to linear scale RGB */
+ _OMP (omp parallel for schedule(static))
+ for (j = 0; j < n; j++)
+ {
+ Y[j] = powf (10,Y[j]);
+
+ rgb[j * 4 + 0] = powf (rgb[j * 4 + 0], saturationFactor) * Y[j];
+ rgb[j * 4 + 1] = powf (rgb[j * 4 + 1], saturationFactor) * Y[j];
+ rgb[j * 4 + 2] = powf (rgb[j * 4 + 2], saturationFactor) * Y[j];
+ }
+ }
+
+ return PFSTMO_OK;
+}
+
+
+static void
+mantiuk06_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
+mantiuk06_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
+mantiuk06_get_cached_region (GeglOperation *operation,
+ const GeglRectangle *roi)
+{
+ return *gegl_operation_source_get_bounding_box (operation, "input");
+}
+
+
+static gboolean
+mantiuk06_process (GeglOperation *operation,
+ GeglBuffer *input,
+ GeglBuffer *output,
+ const GeglRectangle *result)
+{
+ const GeglChantO *o = GEGL_CHANT_PROPERTIES (operation);
+ const gint pix_stride = 4; /* RGBA */
+ gfloat *lum,
+ *pix;
+
+ 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);
+
+ /* Obtain the pixel data */
+ lum = g_new (gfloat, result->width * result->height),
+ gegl_buffer_get (input, 1.0, result, babl_format ("Y float"),
+ lum, 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);
+
+ mantiuk06_contmap (result->width, result->height, pix, lum,
+ o->contrast, o->saturation, FALSE, 200, 1e-3, NULL);
+
+ /* Cleanup and set the output */
+ gegl_buffer_set (output, result, babl_format (OUTPUT_FORMAT), pix,
+ GEGL_AUTO_ROWSTRIDE);
+ g_free (pix);
+ g_free (lum);
+
+ 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 = mantiuk06_process;
+
+ operation_class->prepare = mantiuk06_prepare;
+ operation_class->get_required_for_output = mantiuk06_get_required_for_output;
+ operation_class->get_cached_region = mantiuk06_get_cached_region;
+
+ operation_class->name = "gegl:mantiuk06";
+ 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 constrains "
+ "contrasts across multiple spatial frequencies, producing "
+ "luminance within the range 0.0-1.0");
+}
+
+#endif
+
diff --git a/tests/compositions/mantiuk06.xml b/tests/compositions/mantiuk06.xml
new file mode 100644
index 0000000..114dcb3
--- /dev/null
+++ b/tests/compositions/mantiuk06.xml
@@ -0,0 +1,15 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<gegl>
+ <node operation='gegl:mantiuk06'>
+ <params>
+ <param name='contrast'>0.3</param>
+ <param name='saturation'>0.9</param>
+ <param name='detail'>2.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/mantiuk06.png b/tests/compositions/reference/mantiuk06.png
new file mode 100644
index 0000000..8f4061c
Binary files /dev/null and b/tests/compositions/reference/mantiuk06.png differ
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
