/* 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
 */

/*

    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 <config.h>
#endif /*HAVE_CONFIG_H*/
#include <vips/intl.h>

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <vips/vips.h>

/*
#define DEBUG
 */

/* Our state.
 */
typedef struct {
	DOUBLEMASK *input;	/* Input mask */
	IMAGE *output;		/* Output lut */
	int lut_size;		/* Number of output elements to generate */

	double **data;		/* Rows of unpacked matrix */
	double *buf;		/* Output buffer */
} State;

/* 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 );
}

/* Free our state.
 */
static void
free_state( State *state )
{
	if( state->data ) {
		int i;

		for( i = 0; i < state->input->ysize; i++ )
			if( state->data[i] ) {
				im_free( state->data[i] );
				state->data[i] = NULL;
			}

		im_free( state->data );
		state->data = NULL;
	}
}

/* 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
invertlut( State *state )
{
	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 );
}

/**
 * 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:
 * 
 *   <tgroup cols='4' align='left' colsep='1' rowsep='1'>
 *     <tbody>
 *       <row>
 *         <entry>4</entry>
 *         <entry>3</entry>
 *       </row>
 *       <row>
 *         <entry>0.1</entry>
 *         <entry>0.2</entry>
 *         <entry>0.3</entry>
 *         <entry>0.1</entry>
 *       </row>
 *       <row>
 *         <entry>0.2</entry>
 *         <entry>0.4</entry>
 *         <entry>0.4</entry>
 *         <entry>0.2</entry>
 *       </row>
 *       <row>
 *         <entry>0.7</entry>
 *         <entry>0.5</entry>
 *         <entry>0.6</entry>
 *         <entry>0.3</entry>
 *       </row>
 *     </tbody>
 *   </tgroup>
 *
 * 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: im_buildlut(), im_invertlut()
 *
 * Returns: 0 on success, -1 on error
 */
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 );
}