From 222a7abdbf9fe067bc38c37b67ab75c7e8da5bb1 Mon Sep 17 00:00:00 2001 From: John Cupitt Date: Wed, 4 Sep 2013 14:58:07 +0100 Subject: [PATCH] start on vips_invertlut() --- libvips/create/invertlut.c | 425 +++++++++++++++++++++++++++++++++++++ 1 file changed, 425 insertions(+) create mode 100644 libvips/create/invertlut.c diff --git a/libvips/create/invertlut.c b/libvips/create/invertlut.c new file mode 100644 index 00000000..1c3586ee --- /dev/null +++ b/libvips/create/invertlut.c @@ -0,0 +1,425 @@ +/* invert a lut + * + * Written on: 5/6/01 + * Modified on : + * + * 7/7/03 JC + * - generate image rather than doublemask (arrg) + * 23/3/10 + * - gtkdoc + * 23/5/13 + * - fix 1 high input matrices + * - fix file output + * 4/9/13 + * - convert to a class + */ + +/* + + This file is part of VIPS. + + VIPS 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 2 of the License, or + (at your option) any later version. + + This program 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 this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA + 02110-1301 USA + + */ + +/* + + These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk + + */ + +#ifdef HAVE_CONFIG_H +#include +#endif /*HAVE_CONFIG_H*/ +#include + +#include +#include +#include + +#include + +#include "pcreate.h" + +/* +#define DEBUG + */ + +/* Our state. + */ +typedef struct _VipsInvertlut { + VipsCreate parent_instance; + + /* Input image. + */ + VipsImage *in; + + /* .. and cast to a matrix. + */ + VipsImage *mat; + + int lut_size; /* Number of output elements to generate */ + + double **data; /* Rows of unpacked matrix */ + double *buf; /* Output buffer */ +} VipsInvertlut; + +typedef VipsCreateClass VipsInvertlutClass; + +G_DEFINE_TYPE( VipsInvertlut, vips_invertlut, VIPS_TYPE_CREATE ); + +static void +vips_invertlut_dispose( GObject *gobject ) +{ + VipsBuildlut *lut = (VipsBuildlut *) gobject; + + if( lut->data ) { + int i; + + for( i = 0; i < lut->input->ysize; i++ ) + VIPS_FREE( lut->data[i] ); + + VIPS_FREE( lut->data ); + } + + VIPS_FREE( lut->buf ); + VIPS_UNREF( lut->mat ); + + G_OBJECT_CLASS( vips_invertlut_parent_class )->dispose( gobject ); +} + +/* Use this to sort our input rows by the first column. + */ +static int +compare( const void *a, const void *b ) +{ + double **r1 = (double **) a; + double **r2 = (double **) b; + + double diff = r1[0][0] - r2[0][0]; + + if( diff > 0 ) + return( 1 ); + else if( diff == 0 ) + return( 0 ); + else + return( -1 ); +} + +static int +vips_invertlut_build_init( VipsInvertlut *lut ) +{ + VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( lut ); + +/* Fill our state. + */ +static int +build_state( State *state, DOUBLEMASK *input, IMAGE *output, int lut_size ) +{ + int x, y, i; + + state->input = input; + state->output = output; + state->lut_size = lut_size; + state->data = NULL; + + if( !(state->buf = im_malloc( NULL, IM_IMAGE_SIZEOF_LINE( output ) )) ) + return( -1 ); + + if( !(state->data = IM_ARRAY( NULL, input->ysize, double * )) ) + return( -1 ); + for( y = 0; y < input->ysize; y++ ) + state->data[y] = NULL; + + for( y = 0; y < input->ysize; y++ ) + if( !(state->data[y] = IM_ARRAY( NULL, input->xsize, double )) ) + return( -1 ); + + for( i = 0, y = 0; y < input->ysize; y++ ) + for( x = 0; x < input->xsize; x++, i++ ) + state->data[y][x] = input->coeff[i]; + + /* Sanity check for data range. + */ + for( y = 0; y < input->ysize; y++ ) + for( x = 0; x < input->xsize; x++ ) + if( state->data[y][x] > 1.0 || + state->data[y][x] < 0.0 ) { + im_error( "im_invertlut", + _( "element (%d, %d) is %g, " + "outside range [0,1]" ), + x, y, state->data[y][x] ); + return( -1 ); + } + + /* Sort by 1st column in input. + */ + qsort( state->data, input->ysize, sizeof( double * ), compare ); + +#ifdef DEBUG + printf( "Input table, sorted by 1st column\n" ); + for( y = 0; y < input->ysize; y++ ) { + printf( "%.4d ", y ); + + for( x = 0; x < input->xsize; x++ ) + printf( "%.9f ", state->data[y][x] ); + + printf( "\n" ); + } +#endif /*DEBUG*/ + + return( 0 ); +} + +static int +vips_invertlut_build_create( VipsBuildlut *lut ) +{ + DOUBLEMASK *input = state->input; + int ysize = input->ysize; + int xsize = input->xsize; + double *buf = state->buf; + int bands = xsize - 1; + double **data = state->data; + int lut_size = state->lut_size; + + int b; + + /* Do each output channel separately. + */ + for( b = 0; b < bands; b++ ) { + /* The first and last lut positions we know real values for. + */ + int first = data[0][b + 1] * (lut_size - 1); + int last = data[ysize - 1][b + 1] * (lut_size - 1); + + int k; + + /* Extrapolate bottom and top segments to (0,0) and (1,1). + */ + for( k = 0; k < first; k++ ) { + /* Have this inside the loop to avoid /0 errors if + * first == 0. + */ + double fac = data[0][0] / first; + + buf[b + k * bands] = k * fac; + } + + for( k = last; k < lut_size; k++ ) { + /* Inside the loop to avoid /0 errors for last == + * (lut_size - 1). + */ + double fac = (1 - data[ysize - 1][0]) / + ((lut_size - 1) - last); + + buf[b + k * bands] = + data[ysize - 1][0] + (k - last) * fac; + } + + /* Interpolate the data sections. + */ + for( k = first; k < last; k++ ) { + /* Where we're at in the [0,1] range. + */ + double ki = (double) k / (lut_size - 1); + + double irange, orange; + int j; + + /* Search for the lowest real value < ki. There may + * not be one: if not, just use 0. Tiny error. + */ + for( j = ysize - 1; j >= 0; j-- ) + if( data[j][b + 1] < ki ) + break; + if( j == -1 ) + j = 0; + + /* Interpolate k as being between row data[j] and row + * data[j + 1]. + */ + irange = data[j + 1][b + 1] - data[j][b + 1]; + orange = data[j + 1][0] - data[j][0]; + + buf[b + k * bands] = data[j][0] + + orange * ((ki - data[j][b + 1]) / irange); + } + } + + return( 0 ); +} + +static int +vips_invertlut_build( VipsObject *object ) +{ + VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( object ); + VipsCreate *create = VIPS_CREATE( object ); + VipsBuildlut *lut = (VipsBuildlut *) object; + + if( VIPS_OBJECT_CLASS( vips_invertlut_parent_class )->build( object ) ) + return( -1 ); + + if( vips_check_matrix( class->nickname, lut->in, &lut->mat ) ) + return( -1 ); + + if( vips_invertlut_build_init( lut ) || + vips_invertlut_build_create( lut ) ) + return( -1 ); + + vips_image_init_fields( create->out, + lut->lut_size, 1, lut->mat->Xsize - 1, + VIPS_FORMAT_DOUBLE, VIPS_CODING_NONE, + VIPS_INTERPRETATION_HISTOGRAM, 1.0, 1.0 ); + if( vips_image_write_line( create->out, 0, (VipsPel *) lut->buf ) ) + return( -1 ); + + return( 0 ); +} + +static void +vips_invertlut_class_init( VipsBuildlutClass *class ) +{ + GObjectClass *gobject_class = G_OBJECT_CLASS( class ); + VipsObjectClass *vobject_class = VIPS_OBJECT_CLASS( class ); + + gobject_class->dispose = vips_invertlut_dispose; + gobject_class->set_property = vips_object_set_property; + gobject_class->get_property = vips_object_get_property; + + vobject_class->nickname = "invertlut"; + vobject_class->description = _( "build an inverted look-up table" ); + vobject_class->build = vips_invertlut_build; + + VIPS_ARG_IMAGE( class, "in", 0, + _( "Input" ), + _( "Matrix of XY coordinates" ), + VIPS_ARGUMENT_REQUIRED_INPUT, + G_STRUCT_OFFSET( VipsBuildlut, in ) ); + +} + +static void +vips_invertlut_init( VipsBuildlut *lut ) +{ +} + +/** + * im_invertlut: + * @input: input mask + * @output: output LUT + * @lut_size: generate this much + * + * Given a mask of target values and real values, generate a LUT which + * will map reals to targets. Handy for linearising images from + * measurements of a colour chart. All values in [0,1]. Piecewise linear + * interpolation, extrapolate head and tail to 0 and 1. + * + * Eg. input like this: + * + * + * + * + * 4 + * 3 + * + * + * 0.1 + * 0.2 + * 0.3 + * 0.1 + * + * + * 0.2 + * 0.4 + * 0.4 + * 0.2 + * + * + * 0.7 + * 0.5 + * 0.6 + * 0.3 + * + * + * + * + * Means a patch with 10% reflectance produces an image with 20% in + * channel 1, 30% in channel 2, and 10% in channel 3, and so on. + * + * Inputs don't need to be sorted (we do that). Generate any precision + * LUT, typically you might ask for 256 elements. + * + * It won't work too well for non-monotonic camera responses + * (we should fix this). Interpolation is simple piecewise linear; we ought to + * do something better really. + * + * See also: vips_buildlut(). + * + * Returns: 0 on success, -1 on error + */ +int +vips_invertlut( VipsImage *in, VipsImage **out, ... ) +{ + va_list ap; + int result; + + va_start( ap, out ); + result = vips_call_split( "invertlut", ap, in, out ); + va_end( ap ); + + return( result ); +} + + + + + + + +int +im_invertlut( DOUBLEMASK *input, IMAGE *output, int lut_size ) +{ + State state; + + if( !input || + input->xsize < 2 || + input->ysize < 1 ) { + im_error( "im_invertlut", "%s", _( "bad input matrix" ) ); + return( -1 ); + } + if( lut_size < 1 || + lut_size > 65536 ) { + im_error( "im_invertlut", "%s", _( "bad lut_size" ) ); + return( -1 ); + } + + im_initdesc( output, + lut_size, 1, input->xsize - 1, + IM_BBITS_DOUBLE, IM_BANDFMT_DOUBLE, + IM_CODING_NONE, IM_TYPE_HISTOGRAM, 1.0, 1.0, 0, 0 ); + if( im_setupout( output ) ) + return( -1 ); + + if( build_state( &state, input, output, lut_size ) || + invertlut( &state ) || + im_writeline( 0, output, (VipsPel *) state.buf ) ) { + free_state( &state ); + return( -1 ); + } + free_state( &state ); + + return( 0 ); +}