Menu
▾
▴
[E] /hgrepo/e/en/enblend/enblend: Restore out-of-RGB-cube pixels in ...
|
From: Enblend <enb...@li...> - 2012-01-02 09:27:36
|
branch: details: http://enblend.hg.sourceforge.net/hgweb/enblend/enblend/hgrepo/e/en/enblend/enblend/rev/cf8082d24ef7 changeset: 764:cf8082d24ef7 user: Chris <cs...@us...> date: Mon Jan 02 09:55:21 2012 +0100 description: Restore out-of-RGB-cube pixels in CIECAM02 mode. Implement a "high fidelity" highlight and shadow recovery in CIECAM02 mode. Only deep shadows with a lightness (`J') less than zero are left alone. diffstat: NEWS | 6 + VERSION | 2 +- configure.in | 22 + doc/enblend.info | 0 doc/enfuse.info | 0 doc/versenblend.texi | 4 +- doc/versenfuse.texi | 4 +- src/Makefile.am | 16 +- src/enblend.cc | 2 + src/enfuse.cc | 2 + src/fixmath.h | 206 +++++++++++++++- src/minimizer.h | 638 +++++++++++++++++++++++++++++++++++++++++++++++++++ 12 files changed, 884 insertions(+), 18 deletions(-) diffs (truncated from 1077 to 500 lines): diff -r 04a83babcfc5 -r cf8082d24ef7 NEWS --- a/NEWS Mon Jan 02 09:44:34 2012 +0100 +++ b/NEWS Mon Jan 02 09:55:21 2012 +0100 @@ -25,6 +25,12 @@ blending and the output image is assigned the input images' color profile. +- Enblend and Enfuse exploit a new feature in LittleCMS Version 2.x + called "Unbounded CMM". Thereby, the hue and saturation of extreme + shadows and highlights can be preserved much longer before pure + black or white are reached. See: + http://www.littlecms.com/CIC18_UnboundedCMM.pdf + - Assign different profiles to profile-free input images with option "--fallback-profile" instead of being tied to sRGB. diff -r 04a83babcfc5 -r cf8082d24ef7 VERSION --- a/VERSION Mon Jan 02 09:44:34 2012 +0100 +++ b/VERSION Mon Jan 02 09:55:21 2012 +0100 @@ -1,1 +1,1 @@ -4.1-1864fd62ad72 +4.1-04a83babcfc5 diff -r 04a83babcfc5 -r cf8082d24ef7 configure.in --- a/configure.in Mon Jan 02 09:44:34 2012 +0100 +++ b/configure.in Mon Jan 02 09:55:21 2012 +0100 @@ -20,6 +20,10 @@ AC_C_BIGENDIAN # Checks for libraries. +AC_CHECK_LIB([m], [sqrt]) +AC_CHECK_LIB([gslcblas], [cblas_dgemm]) +AC_CHECK_LIB([gsl], [gsl_blas_dgemm]) + AC_CHECK_LIB(z, gzopen, [], AC_MSG_NOTICE([Compiling without libz.]), []) @@ -200,16 +204,30 @@ AC_MSG_ERROR([Boost "tribool" header file is required to compile Enblend.])) AC_CHECK_HEADER(boost/math/special_functions.hpp, [], AC_MSG_ERROR([Boost "special_functions" header file is required to compile Enblend.])) +AC_CHECK_HEADER(boost/optional.hpp, [], + AC_MSG_ERROR([Boost "optional" header file is required to compile Enblend.])) AC_CHECK_HEADER(boost/pool/pool.hpp, [], AC_MSG_ERROR([Boost "pool" header file is required to compile Enblend.])) AC_CHECK_HEADER(boost/random/mersenne_twister.hpp, [], AC_MSG_ERROR([Boost "mersenne_twister" header file is required to compile Enblend.])) +AC_CHECK_HEADER(boost/random/uniform_real.hpp, [], + AC_MSG_ERROR([Boost "uniform_real" header file is required to compile Enblend.])) +AC_CHECK_HEADER(boost/random/variate_generator.hpp, [], + AC_MSG_ERROR([Boost "variate_generator" header file is required to compile Enblend.])) AC_CHECK_HEADER(boost/scoped_ptr.hpp, [], AC_MSG_ERROR([Boost "scoped_ptr" header file is required to compile Enblend.])) AC_CHECK_HEADER(boost/static_assert.hpp, [], AC_MSG_ERROR([Boost "static_assert" header file is required to compile Enblend.])) AC_CHECK_HEADER(boost/unordered_set.hpp, [], AC_MSG_ERROR([Boost "unordered_set" header file is required to compile Enblend.])) +AC_CHECK_HEADER(gsl/gsl_errno.h, [], + AC_MSG_ERROR([GNU Scientific Library (GSL) header file "gsl_errno" is required to compile Enblend.])) +AC_CHECK_HEADER(gsl/gsl_min.h, [], + AC_MSG_ERROR([GNU Scientific Library (GSL) header file "gsl_min" is required to compile Enblend.])) +AC_CHECK_HEADER(gsl/gsl_multimin.h, [], + AC_MSG_ERROR([GNU Scientific Library (GSL) header file "gsl_multimin" is required to compile Enblend.])) +AC_CHECK_HEADER(gsl/gsl_vector.h, [], + AC_MSG_ERROR([GNU Scientific Library (GSL) header file "gsl_vector" is required to compile Enblend.])) AC_ARG_WITH([boost-filesystem], [AS_HELP_STRING([--with-boost-filesystem], @@ -280,6 +298,10 @@ AC_HEADER_STDBOOL AC_C_CONST AC_C_INLINE +if test "$ac_cv_c_inline" != no; then + AC_DEFINE(HAVE_INLINE, 1, [Define if functions can be declared `inline'.]) + AC_SUBST(HAVE_INLINE) +fi AC_TYPE_OFF_T AC_TYPE_SIGNAL diff -r 04a83babcfc5 -r cf8082d24ef7 doc/enblend.info Binary file doc/enblend.info has changed diff -r 04a83babcfc5 -r cf8082d24ef7 doc/enfuse.info Binary file doc/enfuse.info has changed diff -r 04a83babcfc5 -r cf8082d24ef7 doc/versenblend.texi --- a/doc/versenblend.texi Mon Jan 02 09:44:34 2012 +0100 +++ b/doc/versenblend.texi Mon Jan 02 09:55:21 2012 +0100 @@ -1,4 +1,4 @@ @set UPDATED 19 December 2011 @set UPDATED-MONTH December 2011 -@set EDITION 4.1-1864fd62ad72 -@set VERSION 4.1-1864fd62ad72 +@set EDITION 4.1-04a83babcfc5 +@set VERSION 4.1-04a83babcfc5 diff -r 04a83babcfc5 -r cf8082d24ef7 doc/versenfuse.texi --- a/doc/versenfuse.texi Mon Jan 02 09:44:34 2012 +0100 +++ b/doc/versenfuse.texi Mon Jan 02 09:55:21 2012 +0100 @@ -1,4 +1,4 @@ @set UPDATED 19 December 2011 @set UPDATED-MONTH December 2011 -@set EDITION 4.1-1864fd62ad72 -@set VERSION 4.1-1864fd62ad72 +@set EDITION 4.1-04a83babcfc5 +@set VERSION 4.1-04a83babcfc5 diff -r 04a83babcfc5 -r cf8082d24ef7 src/Makefile.am --- a/src/Makefile.am Mon Jan 02 09:44:34 2012 +0100 +++ b/src/Makefile.am Mon Jan 02 09:55:21 2012 +0100 @@ -16,10 +16,12 @@ filenameparse.h filenameparse.cc \ filespec.h filespec.cc \ self_test.h self_test.cc \ - tiff_message.h tiff_message.cc + tiff_message.h tiff_message.cc \ + minimizer.h enblend_LDFLAGS = $(AM_LDFLAGS) $(OPENGL_CFLAGS) -enblend_LDADD = layer_selection/liblayersel.a vigra_impex/libvigra_impex.a $(OPENGL_LIBS) @EXTRA_LIBS@ -enblend_CXXFLAGS = $(AM_CXXFLAGS) $(OPENGL_CFLAGS) \ +enblend_LDADD = layer_selection/liblayersel.a vigra_impex/libvigra_impex.a \ + $(GSL_LIBS) $(OPENGL_LIBS) @EXTRA_LIBS@ +enblend_CXXFLAGS = $(AM_CXXFLAGS) $(GSL_CFLAGS) $(OPENGL_CFLAGS) \ -D_GNU_SOURCE -D_FILE_OFFSET_BITS=64 \ -I${top_srcdir}/include -I${top_srcdir}/src/layer_selection @@ -30,10 +32,12 @@ filenameparse.h filenameparse.cc \ filespec.h filespec.cc \ self_test.h self_test.cc \ - tiff_message.h tiff_message.cc + tiff_message.h tiff_message.cc \ + minimizer.h enfuse_LDFLAGS = $(AM_LDFLAGS) -enfuse_LDADD = layer_selection/liblayersel.a vigra_impex/libvigra_impex.a @EXTRA_LIBS@ -enfuse_CXXFLAGS = $(AM_CXXFLAGS) \ +enfuse_LDADD = layer_selection/liblayersel.a vigra_impex/libvigra_impex.a \ + $(GSL_LIBS) @EXTRA_LIBS@ +enfuse_CXXFLAGS = $(AM_CXXFLAGS) $(GSL_CFLAGS) \ -D_GNU_SOURCE -D_FILE_OFFSET_BITS=64 \ -I${top_srcdir}/include -I${top_srcdir}/src/layer_selection diff -r 04a83babcfc5 -r cf8082d24ef7 src/enblend.cc --- a/src/enblend.cc Mon Jan 02 09:44:34 2012 +0100 +++ b/src/enblend.cc Mon Jan 02 09:55:21 2012 +0100 @@ -1620,6 +1620,8 @@ sig.initialize(); + gsl_set_error_handler_off(); + TIFFSetWarningHandler(tiff_warning); TIFFSetErrorHandler(tiff_error); diff -r 04a83babcfc5 -r cf8082d24ef7 src/enfuse.cc --- a/src/enfuse.cc Mon Jan 02 09:44:34 2012 +0100 +++ b/src/enfuse.cc Mon Jan 02 09:55:21 2012 +0100 @@ -1597,6 +1597,8 @@ sig.initialize(); + gsl_set_error_handler_off(); + TIFFSetWarningHandler(tiff_warning); TIFFSetErrorHandler(tiff_error); diff -r 04a83babcfc5 -r cf8082d24ef7 src/fixmath.h --- a/src/fixmath.h Mon Jan 02 09:44:34 2012 +0100 +++ b/src/fixmath.h Mon Jan 02 09:55:21 2012 +0100 @@ -28,6 +28,9 @@ #include <cmath> #include <boost/assign/list_of.hpp> +#include <boost/random/mersenne_twister.hpp> +#include <boost/random/uniform_real.hpp> +#include <boost/random/variate_generator.hpp> #ifdef _WIN32 #include <boost/math/special_functions.hpp> @@ -40,6 +43,8 @@ #include "vigra/numerictraits.hxx" #include "vigra/utilities.hxx" +#include "minimizer.h" + #define MAXIMUM_LIGHTNESS 100.0 // J #define MAXIMUM_CHROMA 120.0 // C @@ -96,11 +101,10 @@ static inline void -rgb_to_jch(cmsHTRANSFORM rgb2xyz_transform, cmsHANDLE ciecam_transform, - const double* rgb, cmsJCh* jch) +rgb_to_jch(const double* rgb, cmsJCh* jch) { std::vector<double> xyz(3U); - cmsDoTransform(rgb2xyz_transform, rgb, &xyz[0], 1U); + cmsDoTransform(InputToXYZTransform, rgb, &xyz[0], 1U); const cmsCIEXYZ scaled_xyz = {XYZ_SCALE * xyz[0], XYZ_SCALE * xyz[1], XYZ_SCALE * xyz[2]}; cmsCIECAM02Forward(CIECAMTransform, &scaled_xyz, jch); @@ -109,8 +113,7 @@ static inline void -jch_to_rgb(cmsHANDLE ciecam_transform, cmsHTRANSFORM xyz2rgb_transform, - const cmsJCh* jch, double* rgb) +jch_to_rgb(const cmsJCh* jch, double* rgb) { cmsCIEXYZ scaled_xyz; cmsCIECAM02Reverse(CIECAMTransform, jch, &scaled_xyz); @@ -127,6 +130,78 @@ } +struct extra_minimizer_parameter +{ + extra_minimizer_parameter(const cmsJCh& out_of_box_jch) : jch(out_of_box_jch) + {jch_to_rgb(&jch, bad_rgb);} + + cmsJCh jch; + double bad_rgb[3]; +}; + + +inline double +delta_e_of_rgb(const double* rgb_1, const double* rgb_2) +{ + cmsCIELab lab_1; + cmsCIELab lab_2; + + cmsDoTransform(InputToLabTransform, rgb_1, &lab_1, 1); + cmsDoTransform(InputToLabTransform, rgb_2, &lab_2, 1); + + return cmsCMCdeltaE(&lab_1, &lab_2, 2.0, 1.0); +} + + +inline double +out_of_box_penalty(std::vector<double>::const_iterator first, + std::vector<double>::const_iterator last) +{ + const double infinite_badness = 100.0; + double result = 0.0; + + for (std::vector<double>::const_iterator x = first; x != last; ++x) { + if (*x > 1.0) { + result += *x * infinite_badness; + } else if (*x < 0.0) { + result += (1.0 - *x) * infinite_badness; + } + } + + return result; +} + + +inline double +delta_e_cost(const cmsJCh* jch, const extra_minimizer_parameter* parameter) +{ + std::vector<double> rgb(3U); + jch_to_rgb(jch, &rgb[0]); + + return delta_e_of_rgb(parameter->bad_rgb, &rgb[0]) + out_of_box_penalty(rgb.begin(), rgb.end()); +} + + +double +delta_e_min_cost(double luminance, void* data) +{ + const extra_minimizer_parameter* parameter = static_cast<const extra_minimizer_parameter*>(data); + const cmsJCh jch = {luminance, parameter->jch.C, parameter->jch.h}; + + return delta_e_cost(&jch, parameter); +} + + +double +delta_e_multimin_cost(const gsl_vector* x, void* data) +{ + const extra_minimizer_parameter* parameter = static_cast<const extra_minimizer_parameter*>(data); + const cmsJCh jch = {gsl_vector_get(x, 0), gsl_vector_get(x, 1), parameter->jch.h}; + + return delta_e_cost(&jch, parameter); +} + + /** A functor for converting scalar pixel values to the number representation used * for pyramids. These are either fixed-point integers or floating-point numbers. */ @@ -350,7 +425,7 @@ (scale * NumericTraits<SrcComponentType>::toRealPromote(v.blue())); cmsJCh jch; - rgb_to_jch(InputToXYZTransform, CIECAMTransform, &rgb[0], &jch); + rgb_to_jch(&rgb[0], &jch); // convert cylindrical 'JCh' to cartesian, but reuse (yikes!) the cylindrical structure const double theta = radian_of_degree(jch.h); @@ -383,6 +458,57 @@ } +static inline void +bracket_minimum(const gsl_function& cost, double& x_initial, double x_lower, double x_upper) +{ + const double y_minimum_bound = + std::min(cost.function(x_lower, cost.params), cost.function(x_upper, cost.params)); + double y_initial = cost.function(x_initial, cost.params); + + if (y_initial < y_minimum_bound) { + return; + } + + const unsigned maximum_tries = 100U; + unsigned i = 0U; + + boost::uniform_real<> distribution(std::max(0.001, 1.001 * x_lower), 0.999 * x_upper); + boost::mt19937 seed(1000003); // fixed seed for reproducibility + boost::variate_generator<boost::mt19937&, boost::uniform_real<> > random(seed, distribution); + + while (y_initial >= y_minimum_bound && i < maximum_tries) { + x_initial = random(); + y_initial = cost.function(x_initial, cost.params); + ++i; +#ifdef OPENMP +#pragma omp critical +#endif + std::cout << + "+ highlight recovery -- bracket minimum: x = " << x_initial << ", y = " << y_initial << + std::endl; + } +} + + +static inline double +lightness_guess(const cmsJCh& jch) +{ + return std::min(0.984 * jch.J, // heuristic function with fitted parameter + 0.995 * MAXIMUM_LIGHTNESS); // backstop such that our guess is less than the maximum +} + + +#define LOFI_ERROR (0.5 / 256.0) +#define HIFI_ERROR (0.5 / 65536.0) +#define SUPER_HIFI_ERROR (0.5 / 16777216.0) +#define CIECAM_OPTIMIZER_ERROR HIFI_ERROR + +#define LOFI_GOAL 1.0 +#define HIFI_GOAL 0.5 +#define SUPER_HIFI_GOAL 0.0 +#define CIECAM_OPTIMIZER_GOAL HIFI_GOAL + + /** Fixed point converter that uses ICC profile transformation */ template <typename DestVectorType, typename PyramidVectorType, int PyramidIntegerBits, int PyramidFractionBits> class ConvertJCHPyramidToVectorFunctor { @@ -401,6 +527,10 @@ inline DestVectorType operator()(const PyramidVectorType& v) const { cmsJCh jch = {cf(v.red()), cf(v.green()), cf(v.blue())}; if (jch.J <= 0.0) { +#ifdef OPENMP +#pragma omp critical +#endif + std::cout << "+ unrecoverable dark shadow: J = " << jch.J << "\n" << std::endl; // Lasciate ogne speranza, voi ch'intrate. return DestVectorType(0, 0, 0); } @@ -416,7 +546,69 @@ jch.C = chroma; std::vector<double> rgb(3U); - jch_to_rgb(CIECAMTransform, XYZToInputTransform, &jch, &rgb[0]); + jch_to_rgb(&jch, &rgb[0]); + + if (rgb[0] < 0.0 || rgb[1] < 0.0 || rgb[2] < 0.0) { + extra_minimizer_parameter extra(jch); + gsl_multimin_function cost = {delta_e_multimin_cost, 2U, &extra}; + const MinimizerMultiDimensionSimplex::array_type initial = boost::assign::list_of(jch.J)(jch.C); + const MinimizerMultiDimensionSimplex::array_type step = boost::assign::list_of(1.0)(1.2); + MinimizerMultiDimensionSimplex2Randomized optimizer(cost, initial, step); + + std::vector<double> initial_rgb(3U); // for debug print only + std::copy(rgb.begin(), rgb.end(), initial_rgb.begin()); + + optimizer.set_absolute_error(CIECAM_OPTIMIZER_ERROR)->set_goal(CIECAM_OPTIMIZER_GOAL); + optimizer.run(); + + MinimizerMultiDimensionSimplex::array_type minimum_parameter(2U); + optimizer.x_minimum(minimum_parameter.begin()); + + jch.J = minimum_parameter[0]; + jch.C = minimum_parameter[1]; + jch_to_rgb(&jch, &rgb[0]); +#ifdef OPENMP +#pragma omp critical +#endif + std::cout << + "+ shadow recovery: ini J = " << initial[0] << ", C = " << initial[1] << ", ini RGB = (" << + initial_rgb[0] << ", " << initial_rgb[1] << ", " << initial_rgb[2] << ")\n" << + "+ shadow recovery: opt J = " << minimum_parameter[0] << ", C = " << minimum_parameter[1] << + " after " << optimizer.number_of_iterations() << " iterations\n" << + "+ shadow recovery: opt delta-E = " << optimizer.f_minimum() << + ", opt RGB = (" << rgb[0] << ", " << rgb[1] << ", " << rgb[2] << ")\n" << + std::endl; + } else if (rgb[0] > 1.0 || rgb[1] > 1.0 || rgb[2] > 1.0) { + extra_minimizer_parameter extra(jch); + gsl_function cost = {delta_e_min_cost, &extra}; + double j_initial = lightness_guess(jch); + + bracket_minimum(cost, j_initial, 0.0, MAXIMUM_LIGHTNESS); + GoldenSectionMinimizer1D optimizer(cost, j_initial, 0.0, MAXIMUM_LIGHTNESS); + + const double initial_j = jch.J; // for debug print only + std::vector<double> initial_rgb(3U); // for debug print only + std::copy(rgb.begin(), rgb.end(), initial_rgb.begin()); + + optimizer.set_absolute_error(CIECAM_OPTIMIZER_ERROR)-> + set_goal(CIECAM_OPTIMIZER_GOAL)-> + set_maximum_number_of_iterations(50U); + optimizer.run(); + + jch.J = optimizer.x_minimum(); + jch_to_rgb(&jch, &rgb[0]); +#ifdef OPENMP +#pragma omp critical +#endif + std::cout << + "+ highlight recovery: ini J = " << initial_j << ", {C = " << jch.C << "}, ini RGB = (" << + initial_rgb[0] << ", " << initial_rgb[1] << ", " << initial_rgb[2] << ")\n" << + "+ highlight recovery: opt J = " << optimizer.x_minimum() << + " after " << optimizer.number_of_iterations() << " iterations\n" << + "+ highlight recovery: opt delta-E = " << optimizer.f_minimum() << + ", opt RGB = (" << rgb[0] << ", " << rgb[1] << ", " << rgb[2] << ")\n" << + std::endl; + } limit_sequence(rgb.begin(), rgb.end(), 0.0, 1.0); diff -r 04a83babcfc5 -r cf8082d24ef7 src/minimizer.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/minimizer.h Mon Jan 02 09:55:21 2012 +0100 @@ -0,0 +1,638 @@ +/* + * Copyright (C) 2011 Dr. Christoph L. Spiel + * + * This file is part of Enblend. + * + * Enblend 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. + * + * Enblend 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 Enblend; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + +#ifndef MINIMIZER_H_INCLUDED +#define MINIMIZER_H_INCLUDED + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <vector> + +#include <boost/optional.hpp> + +#include <gsl/gsl_errno.h> +#include <gsl/gsl_min.h> +#include <gsl/gsl_multimin.h> +#include <gsl/gsl_vector.h> + + +class Minimizer +{ + enum {ITERATIONS_PER_DIMENSION = 100U}; + +public: + Minimizer(size_t dimension) : + dimension_(dimension), + maximum_iteration_(ITERATIONS_PER_DIMENSION * dimension), iteration_(0U), + f_goal_(boost::none), absolute_error_(sqrt(std::numeric_limits<double>::epsilon())) + {} + + Minimizer(const Minimizer& minimizer) : + dimension_(minimizer.dimension_), + maximum_iteration_(minimizer.maximum_iteration_), iteration_(minimizer.iteration_), + f_goal_(minimizer.f_goal_), absolute_error_(minimizer.absolute_error_) + {} + + Minimizer& operator=(const Minimizer& minimizer) + { + if (this != &minimizer) + { + dimension_ = minimizer.dimension_; + maximum_iteration_ = minimizer.maximum_iteration_; |