here you go. i ran it against their matlab implementation using octave, the output seemed to match. let me know how you go with it or if you need any help deciphering it!
you may need to clean up the makefile, i was also compiling another test case in this directory. -jo On Thu, Aug 11, 2016 at 4:26 PM, Maximilian Trescher <maximil...@trescherpost.de> wrote: > Hi Johannes, > > Well, I talked to Wolfgang and just started reading this paper: > http://www.di.ens.fr/~aubry/llf.html > to see how they implement the local laplacian. > If you already have code, I would be interested of course and I would try to > incorporate it in the new monochrome module. > > Cheers > Maximilian > > > Am 11.08.2016 um 15:53 schrieb johannes hanika: >> >> heya, >> >> if you have any show-off result images, i'd be interested.. >> >> also, if you need laplacian pyramid code, i think i have a working >> implementation of a local laplacian now (it's monochrome though). >> >> cheers, >> jo >> >> On Wed, Aug 10, 2016 at 10:50 PM, Wolfgang Mader >> <wolfgang_ma...@brain-frog.de> wrote: >>> >>> I finally have something to show. >>> >>> You can grab the code under >>> https://github.com/wmader/darktable >>> the branch in called color2gray. >>> >>> The iop is designed to work with multiple operators. Right now, there are >>> two >>> implemented. The first "lightness" simple keeps L and sets a = b = 0. >>> >>> The second one is called "apparent grayscale" and is more interesting. >>> This >>> method uses a model for the Helmholtz-Kohlrausch (HK) Effect. This effect >>> describes the difference in the perceived brightness of an stimulus >>> depending >>> on its color. The model is described in [1,2]. Basically, the >>> chromaticity of >>> a color is measured as an angle with respect to a gray-point. This angle >>> is >>> the input for the HK-Effect model. There is one free parameter, the >>> ambient >>> luminance, which the paper says should be 20. As there is no sound >>> reasoning >>> for this value, I implemented a slider for it. The effect, however, is >>> minor. >>> >>> I used a test image from [2] to test my implementation against the result >>> of >>> their published gimp-plugin. As they use d65 and we use d50, I patch >>> relevant >>> functions it dt to use d65 for testing. Then, I get matching results. >>> >>> While I implemented the HK model, [2] takes this one step further by >>> enhancing >>> the contrast for areas in the grayscale image where the respective >>> contrast in >>> the color image is larger. I had not yet time to implement this, and for >>> family reasons, I will not have the time in the near future. Fortunately, >>> Maximilian, also on this list, approached me, as he is also interested in >>> the >>> topic. He will take up the task of implementing the chromatic contrast, >>> which >>> is based on Laplacian pyramids. >>> >>> So, check out the code, give it a try. If you have any questions, please >>> let >>> me know, I try to answer in a timely manner. And as I already mentioned, >>> this >>> is my first contribution. Therefore, if the code does not live up to your >>> expectations, let me know, and I try to fix it. >>> >>> Best, >>> Wolfgang >>> >>> P.s.: If anyone is interested in the test images, I am happy to send. I >>> assume, however, that this list does not accept attachments. >>> >>> [1] Simple Estimation Methods for the Helmholtz–Kohlrausch Effect, Y. >>> Nayatani, CIE Journal, 5:2, 1986 >>> >>> [2] Apparent Greyscale: A Simple and Fast Conversion to Perceptually >>> Accurate >>> Images and Video, K. Smith, P. Landes, J. Thollot, K. Myszkowski, >>> EUROGRAPHICS, 27:2, 2008 >>> >>> >>> >>> On Friday, July 1, 2016 8:18:19 PM CEST Wolfgang Mader wrote: >>>> >>>> Hi, >>>> >>>> thanks for asking. I defensed my Ph.D. theses last week and >>>> overestimated >>>> the free time I would have while preparing the defense. Therefore, I >>>> will >>>> only find time for darktable starting with the coming week. However, I >>>> come >>>> back to the list once there are results such that I do not have make >>>> promises I can but live up to. >>>> >>>> Sorry, >>>> Wolfgang >>>> >>>> On June 30, 2016 3:49:19 PM GMT+02:00, Moritz Moeller >>> >>> <virtualr...@gmail.com> wrote: >>>>> >>>>> Hey Wolfang, >>>>> >>>>> On 09/05/16 11:46, Wolfgang Mader wrote: >>>>>> >>>>>> Thanks for the literature. I will have a look at is at soon as >>>>> >>>>> >>>>> possible, and >>>>> >>>>>> will come back to you, once there is something do discuss. >>>>> >>>>> >>>>> I'll have some B&W images to develop soon. Anything I could test >>>>> already? :) >>>>> >>>>> .mm >>>>> >>>>> >>>>> ___________________________________________________________________________ >>>>> darktable developer mailing list >>>>> to unsubscribe send a mail to >>>>> darktable-dev+unsubscr...@lists.darktable.org >>> >>> >> >> ___________________________________________________________________________ >> darktable developer mailing list >> to unsubscribe send a mail to >> darktable-dev+unsubscr...@lists.darktable.org >> > ___________________________________________________________________________ darktable developer mailing list to unsubscribe send a mail to darktable-dev+unsubscr...@lists.darktable.org
#pragma once
#include <string.h>
#include <stdint.h>
#include <stdlib.h>
#include <assert.h>
#include <stdio.h>
#include <math.h>
// downsample width/height to given level
static inline int dl(int size, const int level)
{
for(int l=0;l<level;l++)
size = (size-1)/2+1;
return size;
}
// needs a boundary of 1px around i,j or else it will crash.
static inline float ll_expand_gaussian(
const float *const coarse,
const int i,
const int j,
const int wd,
const int ht)
{
assert(i > 0);
assert(i < wd-1);
assert(j > 0);
assert(j < ht-1);
float c = 0.0f;
const int cw = (wd-1)/2+1;
const float a = 0.4f;
const float w[5] = {1./4.-a/2., 1./4., a, 1./4., 1./4.-a/2.};
const int ind = (j/2)*cw+i/2;
switch((i&1) + 2*(j&1))
{
case 0: // both are even, 3x3 stencil
for(int ii=-1;ii<=1;ii++)
for(int jj=-1;jj<=1;jj++)
c += coarse[ind+cw*jj+ii]*w[2*jj+2]*w[2*ii+2];
break;
case 1: // i is odd, 2x3 stencil
for(int ii=0;ii<=1;ii++)
for(int jj=-1;jj<=1;jj++)
c += coarse[ind+cw*jj+ii]*w[2*jj+2]*w[2*ii+1];
break;
case 2: // j is odd, 3x2 stencil
for(int ii=-1;ii<=1;ii++)
for(int jj=0;jj<=1;jj++)
c += coarse[ind+cw*jj+ii]*w[2*jj+1]*w[2*ii+2];
break;
default: // case 3: // both are odd, 2x2 stencil
for(int ii=0;ii<=1;ii++)
for(int jj=0;jj<=1;jj++)
c += coarse[ind+cw*jj+ii]*w[2*jj+1]*w[2*ii+1];
break;
}
assert(c==c);
return 4.0f * c;
}
// helper to fill in one pixel boundary by copying it
static inline void ll_fill_boundary(
float *const input,
const int wd,
const int ht)
{
for(int j=1;j<ht-1;j++) input[j*wd] = input[j*wd+1];
for(int j=1;j<ht-1;j++) input[j*wd+wd-1] = input[j*wd+wd-2];
memcpy(input, input+wd, sizeof(float)*wd);
memcpy(input+wd*(ht-1), input+wd*(ht-2), sizeof(float)*wd);
}
// static inline float ll_get_laplacian(
static inline void gauss_expand(
const float *const input, // coarse input
float *const fine, // upsampled, blurry output
const int wd, // fine res
const int ht)
{
#ifdef _OPENMP
#pragma omp parallel for default(none) schedule(static) collapse(2)
#endif
for(int j=1;j<ht-1;j++)
for(int i=1;i<wd-1;i++)
fine[j*wd+i] = ll_expand_gaussian(input, i, j, wd, ht);
ll_fill_boundary(fine, wd, ht);
}
static inline void gauss_reduce(
const float *const input, // fine input buffer
float *const coarse, // coarse scale, blurred input buf
const int wd, // fine res
const int ht)
{
// blur, store only coarse res
const int cw = (wd-1)/2+1, ch = (ht-1)/2+1;
memset(coarse, 0, sizeof(float)*cw*ch);
const float a = 0.4f;
const float w[5] = {1./4.-a/2., 1./4., a, 1./4., 1./4.-a/2.};
// direct 5x5 stencil only on required pixels:
for(int j=1;j<ch-1;j++) for(int i=1;i<cw-1;i++)
for(int jj=-2;jj<=2;jj++) for(int ii=-2;ii<=2;ii++)
{
coarse[j*cw+i] += input[(2*j+jj)*wd+2*i+ii] * w[ii+2] * w[jj+2];
assert(
coarse[j*cw+i] ==
coarse[j*cw+i]);
}
ll_fill_boundary(coarse, cw, ch);
}
static inline void rgb_to_yuv(
const float *const rgb,
float *const yuv)
{
const float M[] = {
0.299f, 0.587f, 0.114f,
-0.14713f, -0.28886f, 0.436f,
0.615f, -0.51499f, -0.10001f};
yuv[0] = yuv[1] = yuv[2] = 0.0f;
for(int i=0;i<3;i++)
for(int j=0;j<3;j++)
yuv[i] += M[3*i+j]*rgb[j];
}
static inline float rgb_to_y(
const float *const rgb)
{
const float M[] = {
0.299f, 0.587f, 0.114f,
-0.14713f, -0.28886f, 0.436f,
0.615f, -0.51499f, -0.10001f};
float y = 0.0f;
for(int j=0;j<3;j++)
y += M[j]*rgb[j];
// XXX try shadow lifting: if(y < .5f) return .5f * powf(2.0f*y, 1./1.8);
return y;
}
static inline void yuv_to_rgb(
const float *const yuv,
float *const rgb)
{
const float Mi[] = {
1.0f, 0.0f, 1.13983f,
1.0f, -0.39465f, -0.58060f,
1.0f, 2.03211f, 0.0f};
rgb[0] = rgb[1] = rgb[2] = 0.0f;
for(int i=0;i<3;i++)
for(int j=0;j<3;j++)
rgb[i] += Mi[3*i+j]*yuv[j];
}
// allocate output buffer with monochrome brightness channel from input, padded
// up by max_supp on all four sides, dimensions written to wd2 ht2
static inline float *ll_pad_input(
const float *const input,
const int wd,
const int ht,
const int max_supp,
int *wd2,
int *ht2)
{
const int stride = 3;
*wd2 = 2*max_supp + wd;
*ht2 = 2*max_supp + ht;
// float *out = (float *)dt_alloc_align(64, *wd2**ht2*sizeof(float));
float *out = (float *)malloc(*wd2**ht2*sizeof(float));
for(int j=0;j<ht;j++)
{
for(int i=0;i<max_supp;i++)
// out[(j+max_supp)**wd2+i] = powf(MAX(0.0f, input[stride*wd*j+0]), 1./2.2);
out[(j+max_supp)**wd2+i] = rgb_to_y(input+stride*wd*j+0);
for(int i=0;i<wd;i++)
// out[(j+max_supp)**wd2+i+max_supp] = powf(MAX(0.0f, input[stride*(wd*j+i)+0]), 1./2.2);
out[(j+max_supp)**wd2+i+max_supp] = rgb_to_y(input+stride*(wd*j+i)+0);
for(int i=wd+max_supp;i<*wd2;i++)
// out[(j+max_supp)**wd2+i] = powf(MAX(0.0f, input[stride*(j*wd+wd-1)+0]), 1./2.2);
out[(j+max_supp)**wd2+i] = rgb_to_y(input+stride*(j*wd+wd-1)+0);
}
for(int j=0;j<max_supp;j++)
memcpy(out + *wd2*j, out+max_supp**wd2, sizeof(float)**wd2);
for(int j=max_supp+ht;j<*ht2;j++)
memcpy(out + *wd2*j, out + *wd2*(max_supp+ht-1), sizeof(float)**wd2);
return out;
}
static inline float ll_laplacian(
const float *const coarse, // coarse res gaussian
const float *const fine, // fine res gaussian
const int i, // fine index
const int j,
const int wd, // fine width
const int ht) // fine height
{
const float c = ll_expand_gaussian(coarse, i, j, wd, ht);
return fine[j*wd+i] - c;
}
static inline float ll_curve(
const float x,
const float g,
const float sigma,
const float shadows,
const float highlights,
const float clarity)
{
// this is in the original matlab code instead:
// I_remap=fact*(I-ref).*exp(-(I-ref).*(I-ref)./(2*sigma*sigma));
// also they add up the laplacian to the ones of the base image
// return g + 0.7 * (x-g) + clarity * (x - g) * expf(-(x-g)*(x-g)/(2.0*sigma*sigma));
// XXX highlights does weird things, causing halos in dark regions
// XXX shadows seems to compress highlight contrast (???)
// XXX need to adjust this to bias shad/hi depending on g!
// const float compress = .2;// + powf(g, .4);
// if(x > g+sigma)
// // return g + highlights * (x-g) + clarity * (x - g) * expf(-(x-g)*(x-g)/(2.0*sigma*sigma));
// return g+sigma + highlights * (x-g-sigma) + clarity * (x - g) * expf(-(x-g)*(x-g)/(2.0*sigma*sigma));
// else if(x < g-sigma)
// // return g + shadows * (x-g) + clarity * (x - g) * expf(-(x-g)*(x-g)/(2.0*sigma*sigma));
// return g-sigma + shadows * (x-g+sigma) + clarity * (x - g) * expf(-(x-g)*(x-g)/(2.0*sigma*sigma));
// else
return g + (x-g) + clarity * (x - g) * expf(-(x-g)*(x-g)/(2.0*sigma*sigma));
#if 0
// XXX DEBUG: shad/hi curve needs to be something smart along these lines:
// if(x < .5f)
//{
// const float g2 = powf(g, 1./shadows);
// return powf(x, g);
// return x + g * shadows;//- g + g2;
//}
// centered value
const float c = x-g;
if(c > sigma) return g + sigma + highlights * (c-sigma);
if(c < -sigma) return g - sigma + shadows * (c+sigma);
// const float beta = 0.1;//(1.0-x)*(1.0-x);
// if(c > sigma) return g + sigma + beta * (c-sigma);
// if(c < -sigma) return g - sigma + beta * (c+sigma);
// TODO: paper says to blend this in to avoid boosting noise:
// t = smoothstep(x, 1%, 2%)
const float t = MIN(MAX(0.0f, (x-0.01)/(0.02-0.01)), 1.0f);
assert(t == t);
const float delta = fabsf(c)/sigma;
assert(delta == delta);
const float f_d = t * powf(delta, 1.-clarity) + (1.0f-t)*delta;
assert(f_d == f_d);
assert(g + copysignf(sigma * f_d, c) == g + copysignf(sigma * f_d, c));
return g + copysignf(sigma * f_d, c);
#endif
}
static inline void local_laplacian(
const float *const input, // input buffer in some Labx or yuvx format
float *const out, // output buffer with colour
const int wd, // width and
const int ht, // height of the input buffer
const float sigma, // user param: separate shadows/midtones/highlights
const float shadows, // user param: lift shadows
const float highlights, // user param: compress highlights
const float clarity) // user param: increase clarity/local contrast
{
// XXX TODO: the paper says level 5 is good enough, too?
#define num_levels 8
#define num_gamma 5
const int max_supp = 1<<(num_levels-1);
int w, h;
float *padded[num_levels] = {0};
padded[0] = ll_pad_input(input, wd, ht, max_supp, &w, &h);
// allocate pyramid pointers for padded input
for(int l=1;l<num_levels;l++)
padded[l] = (float *)malloc(sizeof(float)*dl(w,l)*dl(h,l));
// allocate pyramid pointers for output
float *output[num_levels] = {0};
for(int l=0;l<num_levels;l++)
output[l] = (float *)malloc(sizeof(float)*dl(w,l)*dl(h,l));
// create gauss pyramid of padded input, write coarse directly to output
for(int l=1;l<num_levels-1;l++)
gauss_reduce(padded[l-1], padded[l], dl(w,l-1), dl(h,l-1));
gauss_reduce(padded[num_levels-2], output[num_levels-1], dl(w,num_levels-2), dl(h,num_levels-2));
// XXX DEBUG show only detail coeffs
// memset(output[num_levels-1], 0, sizeof(float)*dl(w,num_levels-1)*dl(h,num_levels-1));
// evenly sample brightness [0,1]:
float gamma[num_gamma] = {0.0f};
// for(int k=0;k<num_gamma;k++) gamma[k] = powf(k/(num_gamma-1.0f), 2.0); // XXX DEBUG
// for(int k=0;k<num_gamma;k++) gamma[k] = (k+.5f)/num_gamma;
for(int k=0;k<num_gamma;k++) gamma[k] = k/(num_gamma-1.0f);
// XXX FIXME: don't need to alloc all the memory at once!
// XXX FIXME: accumulate into output pyramid one by one?
// allocate memory for intermediate laplacian pyramids
float *buf[num_gamma][num_levels] = {0};
for(int k=0;k<num_gamma;k++) for(int l=0;l<num_levels;l++)
buf[k][l] = (float *)malloc(sizeof(float)*dl(w,l)*dl(h,l));
// XXX TODO: the paper says remapping only level 3 not 0 does the trick, too:
for(int k=0;k<num_gamma;k++)
{ // process images
#pragma omp parallel for default(none) schedule(dynamic) collapse(2) shared(w,h,k,buf,gamma,padded)
for(int j=0;j<h;j++) for(int i=0;i<w;i++)
buf[k][0][w*j+i] = ll_curve(padded[0][w*j+i], gamma[k], sigma, shadows, highlights, clarity);
// create gaussian pyramids
for(int l=1;l<num_levels;l++)
gauss_reduce(buf[k][l-1], buf[k][l], dl(w,l-1), dl(h,l-1));
}
// assemple output pyramid coarse to fine
for(int l=num_levels-2;l >= 0; l--)
{
const int pw = dl(w,l), ph = dl(h,l);
// XXX i don't think we should manipulate gaussians at all (=>halos) :(
#if 0 // TODO: at a certain level, apply shadow/highlight compression!
if(l==3)
{
const int cw = dl(w,l+1), ch = dl(h,l+1);
#pragma omp parallel for default(none) schedule(static) collapse(2) shared(w,h,output,l,gamma,padded)
for(int j=0;j<ch;j++) for(int i=0;i<cw;i++)
{
// FIXME: this remapping is rubbish
float x = output[l+1][j*cw+i];
// x = 0.5f + 0.2 * (x-0.5f);
// if(x < .5f) x = .5f * powf(MAX(2.0*x,0.0), 1./shadows);
// else x = .5f + .5f * powf(MAX(2.0*(x-.5f),0.0), highlights);
// if(x < .5f) x = .5f * (2.0f*x - powf(1.0-2.0f*x, 2*shadows+1));
if(x< .5f) x = (x - powf(1.-2.*x, 2*(int)shadows+1))/3.+1./3.;
// else x = .5f * (2.0f*x - powf(1.0-2.0f*x, 2*highlights+1));
output[l+1][j*cw+i] = x;
}
}
#endif
gauss_expand(output[l+1], output[l], pw, ph);
// go through all coefficients in the upsampled gauss buffer:
#pragma omp parallel for default(none) schedule(static) collapse(2) shared(w,h,buf,output,l,gamma,padded)
for(int j=1;j<ph-1;j++) for(int i=1;i<pw-1;i++)
{
const float v = padded[l][j*pw+i];
// const float v = output[l][j*pw+i];
#if 1
int hi = 1;
for(;hi<num_gamma-1 && gamma[hi] <= v;hi++);
int lo = hi-1;
const float a = MIN(MAX((v - gamma[lo])/(gamma[hi]-gamma[lo]), 0.0f), 1.0f);
// const float a = (v - gamma[lo])/(gamma[hi]-gamma[lo]);
// printf("val %g %g %g a=%g\n", gamma[lo], v, gamma[hi], a);
#else
int hi = 25, lo = 25;
float a = 0.0f;
#endif
// printf("interpolating %d %d %g (%g)\n", lo, hi, a, v);
// ??? matlab has this: const float li = ll_laplacian(padded[l+1], padded[l], i, j, pw, ph);
const float l0 = ll_laplacian(buf[lo][l+1], buf[lo][l], i, j, pw, ph);
const float l1 = ll_laplacian(buf[hi][l+1], buf[hi][l], i, j, pw, ph);
assert(a == a);
assert(l0 == l0);
assert(l1 == l1);
// XXX is this actually trying to increase contrast of image vs. quantized gammas?
// float laplace = l1 * (1.0f-a) + l0 * a;
float laplace = l0 * (1.0f-a) + l1 * a;
// float laplace = l1;// * (1.0f-a) + l1 * a;
// looks weird: if(laplace < 0.0f) laplace *= 0.2f;
// output[l][j*pw+i] += l0 * (1.0f-a) + l1 * a; // ???
output[l][j*pw+i] += laplace;
}
ll_fill_boundary(output[l], pw, ph);
}
#if 0 // shad/hi by subtracting ll result
#pragma omp parallel for default(none) schedule(static) collapse(2) shared(w,h,output)
for(int j=0;j<ht;j++) for(int i=0;i<wd;i++)
{
const int shad = 1;
float x = output[0][(j+max_supp)*w+max_supp+i];
// else x = .5f * (2.0f*x - powf(1.0-2.0f*x, 2*highlights+1));
// output[0][j*wd+i] = x;
float yuv[3];
rgb_to_yuv(input+3*(j*wd+i), yuv);
const float ydetail = yuv[0] - x;
if(x < .5f) x = (x - powf(1.-2.*x, 2*(int)shad+1))/3.+1./3.;
yuv[0] = x + ydetail;//output[0][(j+max_supp)*w+max_supp+i];
yuv_to_rgb(yuv, out+3*(j*wd+i));
}
#else // output result of ll
#pragma omp parallel for default(none) schedule(dynamic) collapse(2) shared(w,output,buf)
for(int j=0;j<ht;j++) for(int i=0;i<wd;i++)
{
float yuv[3];
rgb_to_yuv(input+3*(j*wd+i), yuv);
yuv[0] = output[0][(j+max_supp)*w+max_supp+i];
yuv_to_rgb(yuv, out+3*(j*wd+i));
// out[3*(j*wd+i)+0] = output[0][(j+max_supp)*w+max_supp+i];
// out[3*(j*wd+i)+1] = buf[num_gamma/2][0][(j+max_supp)*w+max_supp+i]; // XXX DEBUG copy processed channel w/o pyramid
// out[3*(j*wd+i)+2] = input[3*(j*wd+i)+0]; // XXX DEBUG copy original channel
// out[3*(j*wd+i)+2] = input[3*(j*wd+i)+2];
}
#endif
// free all buffers!
for(int l=0;l<num_levels;l++)
{
free(padded[l]);
free(output[l]);
for(int k=0;k<num_gamma;k++)
free(buf[k][l]);
}
#undef num_levels
#undef num_gamma
}
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))
#include "locallaplacian.h"
#include "pfm.h"
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <stdint.h>
#include <string.h>
#include <float.h>
#include <math.h>
int main(int argc, char *argv[])
{
pfm_t in = {0};
pfm_read(&in, argv[1]);
float m = FLT_MAX, M = -FLT_MAX;
for(int k=0;k<in.wd*in.ht;k++) for(int c=0;c<3;c++)
{
in.px[3*k+c] = MIN(1.0, MAX(0.0, in.px[3*k+c]));
m = fminf(m, in.px[3*k+c]);
M = fmaxf(M, in.px[3*k+c]);
}
fprintf(stderr, "range: %g %g\n", m, M);
pfm_t out = {0};
out.wd = in.wd;
out.ht = in.ht;
out.px = (float *)malloc(sizeof(float)*3*out.wd*out.ht);
local_laplacian(
in.px, out.px, in.wd, in.ht,
// sigma, shadows, highlights, clarity
// 0.1, 1.0, 1.0, 5.0); // <= increase clarity
0.2, 1.0, 1.0, -1.0); // <= soften
// 0.10, 1.0, 1.0, .2); // <= lift shadows?
// 0.01, 1.0, 0.0, 0.0); // <= reconstruct h
pfm_write(&out, "locallaplacian.pfm");
#if 0
int wd2, ht2;
float *padded = ll_pad_input(
in.px, in.wd, in.ht, 128, &wd2, &ht2);
out.wd = wd2;
out.ht = ht2;
out.px = (float *)malloc(sizeof(float)*3*out.wd*out.ht);
for(int k=0;k<out.wd*out.ht;k++)
for(int c=0;c<3;c++)
out.px[3*k+c] = padded[k];
pfm_write(&out, "debug.pfm");
float *coarse = (float *)malloc(sizeof(float)*wd2*ht2);
gauss_reduce(padded, coarse, wd2, ht2);
out.wd = wd2/2;
out.ht = ht2/2;
for(int k=0;k<out.wd*out.ht;k++)
for(int c=0;c<3;c++)
out.px[3*k+c] = coarse[k];
pfm_write(&out, "coarse.pfm");
gauss_expand(coarse, padded, wd2, ht2);
out.wd = wd2;
out.ht = ht2;
for(int k=0;k<out.wd*out.ht;k++)
for(int c=0;c<3;c++)
out.px[3*k+c] = padded[k];
pfm_write(&out, "expanded.pfm");
#endif
pfm_clear(&in);
pfm_clear(&out);
exit(0);
}
Makefile
Description: Binary data
#pragma once
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <sys/mman.h>
#include <unistd.h>
typedef struct pfm_t
{
int fd; // file descriptor
void *data; // mapped data
size_t data_size; // size of file
char filename[1024]; // input filename
int wd, ht; // width and height
float *px; // start of pixel data
float exposure; // scale value
float clip; // white point
}
pfm_t;
static inline void pfm_init(
pfm_t *p)
{
memset(p, 0, sizeof(*p));
}
// read pfm file
static inline int pfm_read(
pfm_t *p,
const char *filename)
{
p->fd = -1;
p->data = 0;
if(filename != p->filename)
(void)strncpy(p->filename, filename, sizeof(p->filename));
FILE *f = fopen(p->filename, "rb");
if(!f) return 1;
int w, h;
int r = fscanf(f, "PF\n%d %d\n%*[^\n]", &w, &h);
if(r != 2) { fclose(f); return 1; }
fgetc(f);
if(p->px)
{
if(p->wd != w || p->ht != h) { fclose(f); return 1; }
}
else
{
p->wd = w;
p->ht = h;
p->px = (float *)malloc(sizeof(float)*3*w*h);
}
if(!p->px) { fclose(f); return 1; }
r = fread(p->px, sizeof(float)*3, w*h, f);
if(r != w*h) { fclose(f); return 1; }
fclose(f);
p->clip = 1.0f;
return 0;
}
static inline void pfm_clear(
pfm_t *pfm)
{
if(pfm->data) munmap(pfm->data, pfm->data_size);
if(pfm->fd > 2)
{
close(pfm->fd);
pfm->px = 0;
}
pfm->fd = -1;
free(pfm->px);
pfm->px = 0;
}
// we're likely going out of core when loading all raws at
// the same time:
static inline int pfm_map(
pfm_t *pfm,
const char *filename)
{
pfm->data = 0;
pfm->fd = open(filename, O_RDONLY);
(void)strncpy(pfm->filename, filename, sizeof(pfm->filename));
if(pfm->fd == -1) return 1;
pfm->data_size = lseek(pfm->fd, 0, SEEK_END);
if(pfm->data_size < 100)
{
pfm_clear(pfm);
return 2;
}
lseek(pfm->fd, 0, SEEK_SET);
// this will cause segfaults in case anybody else is writing it while we have it mapped
pfm->data = mmap(0, pfm->data_size, PROT_READ, MAP_SHARED, pfm->fd, 0);
// get pfm header for faster grabbing later on.
// no error handling is done, no comments supported.
char *endptr;
pfm->wd = strtol(pfm->data+3, &endptr, 10);
pfm->ht = strtol(endptr, &endptr, 10);
endptr++; // remove newline
while(endptr < (char *)pfm->data + 100 && *endptr != '\n') endptr++;
pfm->px = (float*)(++endptr); // remove second newline
// while writing make sure the pixel data is 16-byte aligned for sse.
// achieve this by padding up the idiotic scale factor line in the header with additional 0s
if((size_t)pfm->px & 0xf)
fprintf(stderr, "[pfm_map] `%s' pixel buffer not SSE aligned!\n", filename);
if((size_t)pfm->px & 0x3)
fprintf(stderr, "[pfm_map] `%s' pixel buffer not float aligned!\n", filename);
// sanity check:
if(pfm->wd*pfm->ht*3*sizeof(float) > pfm->data_size - (endptr - (char *)pfm->data))
{
pfm_clear(pfm);
return 3;
}
pfm->clip = 1.0f;
return 0;
}
static inline int pfm_write(
pfm_t *p,
const char *filename)
{
FILE *f = fopen(filename, "wb");
if(!f) return 1;
char header[1024];
snprintf(header, 1024, "PF\n%d %d\n-1.0", p->wd, p->ht);
size_t len = strlen(header);
fprintf(f, "PF\n%d %d\n-1.0", p->wd, p->ht);
ssize_t off = 0;
while((len + 1 + off) & 0xf) off++;
while(off-- > 0) fprintf(f, "0");
fprintf(f, "\n");
fwrite(p->px, sizeof(float)*3, p->wd*p->ht, f);
fclose(f);
return 0;
}
