Commit de6649b1 authored by Edouard Gomez's avatar Edouard Gomez

interpolation: rework code structure and API

Since multiple interpolators framework was introduced, inlining
was not possible anymore. So no reason to keep everything in a
header.

- Keep only types (struct and enum), plus two functions as public
  API in header
- Move the code to a proper file
- Modified "public" API a bit preferring passing interpolation entry
  directly instead of its id.
parent a193c2ca
......@@ -29,6 +29,7 @@ FILE(GLOB SOURCE_FILES
"common/imageio_rgbe.c"
"common/imageio_tiff.c"
"common/imageio_rawspeed.cc"
"common/interpolation.c"
"common/metadata.c"
"common/mipmap_cache.c"
"common/styles.c"
......
/* --------------------------------------------------------------------------
This file is part of darktable,
copyright (c) 2012 Edouard Gomez <ed.gomez@free.fr>
darktable 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 3 of the License, or
(at your option) any later version.
darktable 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 darktable. If not, see <http://www.gnu.org/licenses/>.
* ------------------------------------------------------------------------*/
#include "common/interpolation.h"
#include "control/conf.h"
#include <math.h>
#include <stddef.h>
#include <stdint.h>
#include <glib.h>
#include <assert.h>
/* --------------------------------------------------------------------------
* Functions
* ------------------------------------------------------------------------*/
static inline float
_dt_interpolation_func_bilinear(float width, float t)
{
float r;
t = fabsf(t);
if (t>1.f) {
r = 0.f;
} else {
r = 1.f - t;
}
return r;
}
static inline float
_dt_interpolation_func_bicubic(float width, float t)
{
float r;
t = fabsf(t);
if (t>=2.f) {
r = 0.f;
} else if (t>1.f && t<2.f) {
float t2 = t*t;
r = 0.5f*(t*(-t2 + 5.f*t - 8.f) + 4.f);
} else {
float t2 = t*t;
r = 0.5f*(t*(3.f*t2 - 5.f*t) + 2.f);
}
return r;
}
#define DT_LANCZOS_EPSILON (1e-9f)
#if 0
// Canonic version
static inline float
_dt_interpolation_func_lanczos(float width, float t)
{
float r;
if (t<-width || t>width) {
r = 0.f;
} else if (t>-DT_LANCZOS_EPSILON && t<DT_LANCZOS_EPSILON) {
r = 1.f;
} else {
r = width*sinf(M_PI*t)*sinf(M_PI*t/width)/(M_PI*M_PI*t*t);
}
return r;
}
#else
/* Fast lanczos version, no calls to math.h functions, too accurate, too slow
*
* Based on a forum entry at
* http://devmaster.net/forums/topic/4648-fast-and-accurate-sinecosine/
*
* Apart the fast sine function approximation, the only trick is to compute:
* sin(pi.t) = sin(a.pi + r.pi) where t = a + r = trunc(t) + r
* = sin(a.pi).cos(r.pi) + sin(r.pi).cos(a.pi)
* = 0*cos(r.pi) + sin(r.pi).cos(a.pi)
* = sign.sin(r.pi) where sign = 1 if the a is even
* = -1 if the a is odd
*
* Of course we know that lanczos func will only be called for
* the range -width < t < width so we can additionally avoid the
* range check. */
// Valid for [-pi pi] only
static inline float
_dt_sinf_fast(float t)
{
static const float a = 4/(M_PI*M_PI);
static const float p = 0.225f;
t = a*t*(M_PI - fabsf(t));
return t*(p*(fabsf(t) - 1) + 1);
}
static inline float
_dt_interpolation_func_lanczos(float width, float t)
{
/* Compute a value for sinf(pi.t) in [-pi pi] for which the value will be
* correct */
int a = (int)t;
float r = t - (float)a;
// Compute the correct sign for sinf(pi.r)
union { float f; uint32_t i; } sign;
sign.i = ((a&1)<<31) | 0x3f800000;
return (DT_LANCZOS_EPSILON + width*sign.f*_dt_sinf_fast(M_PI*r)*_dt_sinf_fast(M_PI*t/width))/(DT_LANCZOS_EPSILON + M_PI*M_PI*t*t);
}
#endif
#undef DT_LANCZOS_EPSILON
/* --------------------------------------------------------------------------
* Interpolators
* ------------------------------------------------------------------------*/
static const struct dt_interpolation dt_interpolator[] =
{
{DT_INTERPOLATION_BILINEAR, "bilinear", 1, &_dt_interpolation_func_bilinear},
{DT_INTERPOLATION_BICUBIC, "bicubic", 2, &_dt_interpolation_func_bicubic},
{DT_INTERPOLATION_LANCZOS2, "lanczos2", 2, &_dt_interpolation_func_lanczos},
{DT_INTERPOLATION_LANCZOS3, "lanczos3", 3, &_dt_interpolation_func_lanczos},
};
/* --------------------------------------------------------------------------
* Kernel utility method
* ------------------------------------------------------------------------*/
static inline float
_dt_interpolation_compute_kernel(
const struct dt_interpolation* itor,
float* kernel,
float t)
{
/* Find closest integer position and then offset that to match first
* filtered sample position */
t = t - (float)((int)t) + (float)itor->width - 1.f;
// Will hold kernel norm
float norm = 0.f;
// Compute the raw kernel
for (int i=0; i<2*itor->width; i++) {
float tap = itor->func((float)itor->width, t);
norm += tap;
kernel[i] = tap;
t -= 1.f;
}
return norm;
}
/* --------------------------------------------------------------------------
* Interpolation function (see usage in iop/lens.c and iop/clipping.c)
* ------------------------------------------------------------------------*/
float
dt_interpolation_compute_sample(
const struct dt_interpolation* itor,
const float* in,
const float x, const float y,
const int samplestride, const int linestride)
{
assert(itor->width < 4);
float kernelh[6];
float kernelv[6];
// Compute both horizontal and vertical kernels
float normh = _dt_interpolation_compute_kernel(itor, kernelh, x);
float normv = _dt_interpolation_compute_kernel(itor, kernelv, y);
// Go to top left pixel
in = in - (itor->width-1)*(samplestride + linestride);
// Apply the kernel
float s = 0.f;
for (int i=0; i<2*itor->width; i++) {
float h = 0.0f;
for (int j=0; j<2*itor->width; j++) {
h += kernelh[j]*in[j*samplestride];
}
s += kernelv[i]*h;
in += linestride;
}
return s/(normh*normv);
}
const struct dt_interpolation*
dt_interpolation_new(
enum dt_interpolation_type type)
{
const struct dt_interpolation* itor = NULL;
if (type == DT_INTERPOLATION_USERPREF) {
// Find user preferred interpolation method
gchar* uipref = dt_conf_get_string("plugins/lighttable/export/pixel_interpolator");
for (int i=DT_INTERPOLATION_FIRST; uipref && i<DT_INTERPOLATION_LAST; i++) {
if (!strcmp(uipref, dt_interpolator[i].name)) {
// Found the one
itor = &dt_interpolator[i];
break;
}
}
g_free(uipref);
/* In the case the search failed (!uipref or name not found),
* prepare later search pass with default fallback */
type = DT_INTERPOLATION_DEFAULT;
}
if (!itor) {
// Did not find the userpref one or we've been asked for a specific one
for (int i=DT_INTERPOLATION_FIRST; i<DT_INTERPOLATION_LAST; i++) {
if (dt_interpolator[i].id == type) {
itor = &dt_interpolator[i];
break;
}
if (dt_interpolator[i].id == DT_INTERPOLATION_DEFAULT) {
itor = &dt_interpolator[i];
}
}
}
return itor;
}
......@@ -19,222 +19,65 @@
#ifndef INTERPOLATION_H
#define INTERPOLATION_H
#include <assert.h>
/* --------------------------------------------------------------------------
* Types
* ------------------------------------------------------------------------*/
/** Interpolation function */
typedef float (*dt_interpolation_func)(float width, float t);
/** available interpolations */
enum dt_interpolation
/** Available interpolations */
enum dt_interpolation_type
{
DT_INTERPOLATION_FIRST=0,
DT_INTERPOLATION_BILINEAR=DT_INTERPOLATION_FIRST,
DT_INTERPOLATION_BICUBIC,
DT_INTERPOLATION_LANCZOS2,
DT_INTERPOLATION_LANCZOS3,
DT_INTERPOLATION_LAST,
DT_INTERPOLATION_FIRST=0, /**< Helper for easy iteration on interpolators */
DT_INTERPOLATION_BILINEAR=DT_INTERPOLATION_FIRST, /**< Bilinear interpolation (aka tent filter) */
DT_INTERPOLATION_BICUBIC, /**< Bicubic interpolation (with -0.5 parameter) */
DT_INTERPOLATION_LANCZOS2, /**< Lanczos interpolation (with 2 lobes) */
DT_INTERPOLATION_LANCZOS3, /**< Lanczos interpolation (with 3 lobes) */
DT_INTERPOLATION_LAST, /**< Helper for easy iteration on interpolators */
DT_INTERPOLATION_DEFAULT=DT_INTERPOLATION_BILINEAR,
DT_INTERPOLATION_USERPREF /**< can be specified so that user setting is chosen */
};
/** Interpolation description */
struct dt_interpolation_desc
{
enum dt_interpolation id;
const char* name;
int width;
dt_interpolation_func func;
};
/* --------------------------------------------------------------------------
* Functions
* ------------------------------------------------------------------------*/
static inline float
_dt_interpolation_func_bilinear(float width, float t)
{
float r;
t = fabsf(t);
if (t>1.f) {
r = 0.f;
} else {
r = 1.f - t;
}
return r;
}
static inline float
_dt_interpolation_func_bicubic(float width, float t)
{
float r;
t = fabsf(t);
if (t>=2.f) {
r = 0.f;
} else if (t>1.f && t<2.f) {
float t2 = t*t;
r = 0.5f*(t*(-t2 + 5.f*t - 8.f) + 4.f);
} else {
float t2 = t*t;
r = 0.5f*(t*(3.f*t2 - 5.f*t) + 2.f);
}
return r;
}
#define DT_LANCZOS_EPSILON (1e-9f)
/** Interpolation function */
typedef float (*dt_interpolation_func)(float width, float t);
#if 0
// Canonic version
static inline float
_dt_interpolation_func_lanczos(float width, float t)
/** Interpolation structure */
struct dt_interpolation
{
float r;
enum dt_interpolation_type id; /**< Id such as defined by the dt_interpolation_type */
const char* name; /**< internal name */
int width; /**< Half width of its kernel support */
dt_interpolation_func func; /**< Kernel function */
};
if (t<-width || t>width) {
r = 0.f;
} else if (t>-DT_LANCZOS_EPSILON && t<DT_LANCZOS_EPSILON) {
r = 1.f;
} else {
r = width*sinf(M_PI*t)*sinf(M_PI*t/width)/(M_PI*M_PI*t*t);
}
return r;
}
#else
/* Fast lanczos version, no calls to math.h functions, too accurate, too slow
/** Compute a single interpolated sample.
*
* Based on a forum entry at
* http://devmaster.net/forums/topic/4648-fast-and-accurate-sinecosine/
* This function computes a single interpolated sample. Implied costs are:
* <ul>
* <li>Horizontal filtering kernel computation</li>
* <li>Vertical filtering kernel computation</li>
* <li>Sample computation</li>
* </ul>
*
* Apart the fast sine function approximation, the only trick is to compute:
* sin(pi.t) = sin(a.pi + r.pi) where t = a + r = trunc(t) + r
* = sin(a.pi).cos(r.pi) + sin(r.pi).cos(a.pi)
* = 0*cos(r.pi) + sin(r.pi).cos(a.pi)
* = sign.sin(r.pi) where sign = 1 if the a is even
* = -1 if the a is odd
* @param in Input samples
* @param itor interpolator to be used
* @param x X-Coordinate of the requested sample
* @param y Y-Coordinate of the requested sample
* @param samplestride Stride in bytes for a sample
* @param linestride Stride in bytes for complete line
*
* Of course we know that lanczos func will only be called for
* the range -width < t < width so we can additionally avoid the
* range check. */
// Valid for [-pi pi] only
static inline float
_dt_sinf_fast(float t)
{
static const float a = 4/(M_PI*M_PI);
static const float p = 0.225f;
t = a*t*(M_PI - fabsf(t));
return t*(p*(fabsf(t) - 1) + 1);
}
static inline float
_dt_interpolation_func_lanczos(float width, float t)
{
/* Compute a value for sinf(pi.t) in [-pi pi] for which the value will be
* correct */
int a = (int)t;
float r = t - (float)a;
// Compute the correct sign for sinf(pi.r)
union { float f; uint32_t i; } sign;
sign.i = ((a&1)<<31) | 0x3f800000;
return (DT_LANCZOS_EPSILON + width*sign.f*_dt_sinf_fast(M_PI*r)*_dt_sinf_fast(M_PI*t/width))/(DT_LANCZOS_EPSILON + M_PI*M_PI*t*t);
}
#endif
#undef DT_LANCZOS_EPSILON
/* --------------------------------------------------------------------------
* Interpolators
* ------------------------------------------------------------------------*/
static const struct dt_interpolation_desc dt_interpolator[] =
{
{DT_INTERPOLATION_BILINEAR, "bilinear", 1, &_dt_interpolation_func_bilinear},
{DT_INTERPOLATION_BICUBIC, "bicubic", 2, &_dt_interpolation_func_bicubic},
{DT_INTERPOLATION_LANCZOS2, "lanczos2", 2, &_dt_interpolation_func_lanczos},
{DT_INTERPOLATION_LANCZOS3, "lanczos3", 3, &_dt_interpolation_func_lanczos},
};
/* --------------------------------------------------------------------------
* Kernel utility method
* ------------------------------------------------------------------------*/
static inline float
_dt_interpolation_compute_kernel(enum dt_interpolation itype, float* kernel, float t)
{
const struct dt_interpolation_desc* interpolator = &dt_interpolator[itype];
/* Find closest integer position and then offset that to match first
* filtered sample position */
t = t - (float)((int)t) + (float)interpolator->width - 1.f;
// Will hold kernel norm
float norm = 0.f;
// Compute the raw kernel
for (int i=0; i<2*interpolator->width; i++) {
float tap = interpolator->func((float)interpolator->width, t);
norm += tap;
kernel[i] = tap;
t -= 1.f;
}
return norm;
}
/* --------------------------------------------------------------------------
* Interpolation function (see usage in iop/lens.c and iop/clipping.c)
* ------------------------------------------------------------------------*/
static inline float
dt_interpolation_compute(const float* in, const float x, const float y, enum dt_interpolation itype, const int samplestride, const int linestride)
{
const struct dt_interpolation_desc* interpolator = &dt_interpolator[itype];
assert(interpolator->width < 4);
float kernelh[6];
float kernelv[6];
// Compute both horizontal and vertical kernels
float normh = _dt_interpolation_compute_kernel(itype, kernelh, x);
float normv = _dt_interpolation_compute_kernel(itype, kernelv, y);
// Go to top left pixel
in = in - (interpolator->width-1)*(samplestride + linestride);
// Apply the kernel
float s = 0.f;
for (int i=0; i<2*interpolator->width; i++) {
float h = 0.0f;
for (int j=0; j<2*interpolator->width; j++) {
h += kernelh[j]*in[j*samplestride];
}
s += kernelv[i]*h;
in += linestride;
}
return s/(normh*normv);
}
static inline enum dt_interpolation
dt_interpolation_get_type()
{
enum dt_interpolation itype = DT_INTERPOLATION_DEFAULT;
gchar* uipref = dt_conf_get_string("plugins/lighttable/export/pixel_interpolator");
for (int i=DT_INTERPOLATION_FIRST; uipref && i<DT_INTERPOLATION_LAST; i++) {
if (!strcmp(uipref, dt_interpolator[i].name)) {
itype = dt_interpolator[i].id;
break;
}
}
g_free(uipref);
return itype;
}
* @return computed sample
*/
float
dt_interpolation_compute_sample(
const struct dt_interpolation* itor,
const float* in,
const float x,
const float y,
const int samplestride,
const int linestride);
/** Get an interpolator from type
* @param type Interpolator to search for
* @return requested interpolator or default if not found (this function can't fail)
*/
const struct dt_interpolation*
dt_interpolation_new(
enum dt_interpolation_type type);
#endif /* INTERPOLATION_H */
......@@ -350,8 +350,7 @@ void modify_roi_in(struct dt_iop_module_t *self, struct dt_dev_pixelpipe_iop_t *
// adjust roi_in to minimally needed region
// +/-1 stands for imprecision in transform
enum dt_interpolation itype = dt_interpolation_get_type();
const struct dt_interpolation_desc* interpolation = &dt_interpolator[itype];
const struct dt_interpolation* interpolation = dt_interpolation_new(DT_INTERPOLATION_USERPREF);
roi_in->x = aabb_in[0] - interpolation->width - 1;
roi_in->y = aabb_in[1] - interpolation->width - 1;
roi_in->width = aabb_in[2]-aabb_in[0]+2*(interpolation->width+1);
......@@ -402,11 +401,10 @@ void process (struct dt_iop_module_t *self, dt_dev_pixelpipe_iop_t *piece, void
}
else
{
enum dt_interpolation itype = dt_interpolation_get_type();
const struct dt_interpolation_desc* interpolation = &dt_interpolator[itype];
const struct dt_interpolation* interpolation = dt_interpolation_new(DT_INTERPOLATION_USERPREF);
#ifdef _OPENMP
#pragma omp parallel for schedule(static) default(none) shared(d,ivoid,ovoid,roi_in,roi_out,itype,interpolation)
#pragma omp parallel for schedule(static) default(none) shared(d,ivoid,ovoid,roi_in,roi_out,interpolation)
#endif
// (slow) point-by-point transformation.
// TODO: optimize with scanlines and linear steps between?
......@@ -446,7 +444,7 @@ void process (struct dt_iop_module_t *self, dt_dev_pixelpipe_iop_t *piece, void
{
const float *in = ((float *)ivoid) + ch*(roi_in->width*jj+ii);
for(int c=0; c<3; c++,in++)
out[c] = dt_interpolation_compute(in, po[0], po[1], itype, ch, ch_width);
out[c] = dt_interpolation_compute_sample(interpolation, in, po[0], po[1], ch, ch_width);
}
else for(int c=0; c<3; c++) out[c] = 0.0f;
}
......@@ -479,10 +477,10 @@ process_cl (struct dt_iop_module_t *self, dt_dev_pixelpipe_iop_t *piece, cl_mem
}
else
{
enum dt_interpolation itype = dt_interpolation_get_type();
const struct dt_interpolation* interpolation = dt_interpolation_new(DT_INTERPOLATION_USERPREF);
// no opencl support for higher level interpolation yet
if(itype != DT_INTERPOLATION_BILINEAR) return FALSE;
if(interpolation->id != DT_INTERPOLATION_BILINEAR) return FALSE;
int roi[2] = { roi_in->x, roi_in->y };
int roo[2] = { roi_out->x, roi_out->y };
......
......@@ -145,11 +145,10 @@ process (struct dt_iop_module_t *self, dt_dev_pixelpipe_iop_t *piece, void *ivoi
d->tmpbuf2 = (float *)dt_alloc_align(16, d->tmpbuf2_len);
}
enum dt_interpolation itype = dt_interpolation_get_type();
const struct dt_interpolation_desc* interpolation = &dt_interpolator[itype];
const struct dt_interpolation* interpolation = dt_interpolation_new(DT_INTERPOLATION_USERPREF);
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(roi_out, roi_in, in, d, ovoid, modifier, itype, interpolation) schedule(static)
#pragma omp parallel for default(none) shared(roi_out, roi_in, in, d, ovoid, modifier, interpolation) schedule(static)
#endif
for (int y = 0; y < roi_out->height; y++)
{
......@@ -167,7 +166,7 @@ process (struct dt_iop_module_t *self, dt_dev_pixelpipe_iop_t *piece, void *ivoi
if(ii >= (interpolation->width-1) && jj >= (interpolation->width-1) && ii < roi_in->width-interpolation->width && jj < roi_in->height-interpolation->width)
{
const float *inp = in + ch*(roi_in->width*jj+ii) + c;
buf[c] = dt_interpolation_compute(inp, pi0, pi1, itype, ch, ch_width);
buf[c] = dt_interpolation_compute_sample(interpolation, inp, pi0, pi1, ch, ch_width);
}
else buf[c] = 0.0f;
}
......@@ -238,11 +237,10 @@ process (struct dt_iop_module_t *self, dt_dev_pixelpipe_iop_t *piece, void *ivoi
d->tmpbuf2 = (float *)dt_alloc_align(16, d->tmpbuf2_len);
}
enum dt_interpolation itype = dt_interpolation_get_type();
const struct dt_interpolation_desc* interpolation = &dt_interpolator[itype];
const struct dt_interpolation* interpolation = dt_interpolation_new(DT_INTERPOLATION_USERPREF);
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(roi_in, roi_out, d, ovoid, modifier, itype, interpolation) schedule(static)
#pragma omp parallel for default(none) shared(roi_in, roi_out, d, ovoid, modifier, interpolation) schedule(static)
#endif
for (int y = 0; y < roi_out->height; y++)
{
......@@ -260,7 +258,7 @@ process (struct dt_iop_module_t *self, dt_dev_pixelpipe_iop_t *piece, void *ivoi
if(ii >= (interpolation->width-1) && jj >= (interpolation->width-1) && ii < roi_in->width-interpolation->width && jj < roi_in->height-interpolation->width)
{
const float *inp = d->tmpbuf + ch*(roi_in->width*jj+ii)+c;
out[c] = dt_interpolation_compute(inp, pi0, pi1, itype, ch, ch_width);
out[c] = dt_interpolation_compute_sample(interpolation, inp, pi0, pi1, ch, ch_width);
}
else out[c] = 0.0f;
}
......@@ -326,10 +324,10 @@ process_cl (struct dt_iop_module_t *self, dt_dev_pixelpipe_iop_t *piece, cl_mem
return TRUE;
}
enum dt_interpolation itype = dt_interpolation_get_type();
const struct dt_interpolation* interpolation = dt_interpolation_new(DT_INTERPOLATION_USERPREF);
// no opencl support for higher level interpolation yet
if(itype != DT_INTERPOLATION_BILINEAR) return FALSE;
if(interpolation->id != DT_INTERPOLATION_BILINEAR) return FALSE;
tmpbuf = (float *)dt_alloc_align(16, tmpbuflen);
if(tmpbuf == NULL) goto error;
......@@ -583,9 +581,7 @@ void modify_roi_in(struct dt_iop_module_t *self, struct dt_dev_pixelpipe_iop_t *
}
}
enum dt_interpolation itype = dt_interpolation_get_type();
const struct dt_interpolation_desc* interpolation = &dt_interpolator[itype];
const struct dt_interpolation* interpolation = dt_interpolation_new(DT_INTERPOLATION_USERPREF);
roi_in->x = fmaxf(0.0f, xm-interpolation->width);
roi_in->y = fmaxf(0.0f, ym-interpolation->width);
roi_in->width = fminf(orig_w-roi_in->x, xM - roi_in->x + interpolation->width);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment