/* @(#) Map an image through another image, acting as a LUT (Look Up Table). 
 * @(#) The lut may have any type, and the output image will be that type.
 * @(#) 
 * @(#) The input image must be an unsigned integer types, that is, it must
 * @(#) be one of IM_BANDFMT_UCHAR, IM_BANDFMT_USHORT or IM_BANDFMT_UINT.
 * @(#) 
 * @(#) If the input is IM_BANDFMT_UCHAR, then the LUT must have 256 elements, 
 * @(#) in other words, lut->Xsize * lut->Ysize == 256.
 * @(#)  
 * @(#) If the input is IM_BANDFMT_USHORT or IM_BANDFMT_UINT, then the lut 
 * @(#) may have any number of elements, and input pels whose value is 
 * @(#) greater than lut->Xsize * lut->Ysize are mapped with the last LUT 
 * @(#) element.
 * @(#) 
 * @(#) If LUT has one band, then all bands of input pass through it. If LUT
 * @(#) has same number of bands as input, then each band is LUTed
 * @(#) separately. If input has one band, then LUT may have many bands and
 * @(#) the output will have the same number of bands as the LUT.
 * @(#)
 * @(#) int 
 * @(#) im_maplut( in, out, lut )
 * @(#) IMAGE *in, *out, *lut;
 * @(#)
 * @(#) Returns 0 on success and -1 on error
 * @(#)
 *
 * Modified:
 * 18/6/93 JC
 *	- oops! im_incheck() added for LUT image
 * 	- some ANSIfication
 * 15/7/93 JC
 *	- adapted for partial v2
 *	- ANSIfied
 *	- now does complex LUTs too
 * 10/3/94 JC
 *	- more helpful error messages, slight reformatting
 * 24/8/94 JC
 *	- now allows non-uchar image input
 * 7/10/94 JC
 *	- uses im_malloc(), IM_NEW() etc.
 * 13/3/95 JC
 *	- now takes a private copy of LUT, so user can im_close() LUT image
 *	  after im_maplut() without fear of coredumps
 * 23/6/95 JC
 *	- lut may now have many bands if image has just one band
 * 3/3/01 JC
 *	- small speed ups
 * 30/6/04
 *	- heh, 1 band image + 3 band lut + >8bit output has been broken for 9
 *	years :-)
 */

/*

    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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  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 <string.h>

#include <vips/vips.h>

#ifdef WITH_DMALLOC
#include <dmalloc.h>
#endif /*WITH_DMALLOC*/

/* Struct we carry for LUT operations.
 */
typedef struct {
	int fmt;		/* LUT image BandFmt */
	int nb;			/* Number of bands in lut */
	int es;			/* IM_IMAGE_SIZEOF_ELEMENT() for lut image */
	int sz;			/* Number of elements in minor dimension */
	int clp;		/* Value we clip against */
	PEL **table;		/* Lut converted to 2d array */
	int overflow;		/* Number of overflows for non-uchar lut */
} LutInfo;

/* Print overflows, if any.
 */
static int
end_lut( LutInfo *st )
{
	if( st->overflow ) {
		im_warn( "im_maplut", _( "%d overflows detected" ), 
			st->overflow );
		st->overflow = 0;
	}

	return( 0 );
}

/* Build a lut table.
 */
static LutInfo *
build_luts( IMAGE *out, IMAGE *lut )
{
	LutInfo *st = IM_NEW( out, LutInfo );
	int i, x;
	PEL *q;

	if( !st )
                return( NULL );

	/* Make luts. We unpack the LUT image into a C 2D array to speed
	 * processing.
	 */
	st->fmt = lut->BandFmt;
	st->es = IM_IMAGE_SIZEOF_ELEMENT( lut );
	st->nb = lut->Bands;
	st->sz = lut->Xsize * lut->Ysize;
	st->clp = st->sz - 1;
	st->overflow = 0;
	st->table = NULL;
	if( im_add_evalend_callback( out, 
		(im_callback_fn) end_lut, st, NULL ) ) 
		return( NULL );

	/* Attach tables.
	 */
	if( !(st->table = IM_ARRAY( out, lut->Bands, PEL * )) ) 
                return( NULL );
	for( i = 0; i < lut->Bands; i++ )
		if( !(st->table[i] = IM_ARRAY( out, st->sz * st->es, PEL )) )
			return( NULL );

	/* Scan LUT and fill table.
	 */
	q = (PEL *) lut->data;
	for( x = 0; x < st->sz; x++ )
		for( i = 0; i < st->nb; i++ ) {
			memcpy( st->table[i] + x * st->es, q, st->es );
			q += st->es;
		}
	
	return( st );
}

/* Our sequence value: the region this sequence is using, and local stats.
 */
typedef struct {
	REGION *ir;		/* Input region */
	int overflow;		/* Number of overflows */
} Seq;

/* Destroy a sequence value.
 */
static int
stop_maplut( Seq *seq, IMAGE *in, LutInfo *st )
{
	/* Add to global stats.
	 */
	st->overflow += seq->overflow;
	
	if( seq->ir ) {
		im_region_free( seq->ir );
		seq->ir = NULL;
	}

	return( 0 );
}

/* Our start function.
 */
static void *
start_maplut( IMAGE *out, IMAGE *in, LutInfo *st )
{
	Seq *seq = IM_NEW( out, Seq );
	 
	if( !seq )
		 return( NULL );

	/* Init!
	 */
	seq->ir = NULL;
	seq->overflow = 0;

	if( !(seq->ir = im_region_create( in )) ) 
		return( NULL );

	return( (void *) seq );
}

/* Map through n non-complex luts.
 */
#define loop(OUT) { \
	int b = st->nb; \
	\
	for( y = to; y < bo; y++ ) { \
		for( z = 0; z < b; z++ ) { \
			PEL *p = (PEL *) IM_REGION_ADDR( ir, le, y ); \
			OUT *q = (OUT *) IM_REGION_ADDR( or, le, y ); \
			OUT *tlut = (OUT *) st->table[z]; \
			\
			for( x = z; x < ne; x += b ) \
				q[x] = tlut[p[x]]; \
		} \
	} \
}

/* Map through n complex luts.
 */
#define loopc(OUT) { \
	int b = in->Bands; \
	\
	for( y = to; y < bo; y++ ) { \
		for( z = 0; z < b; z++ ) { \
			PEL *p = (PEL *) IM_REGION_ADDR( ir, le, y ) + z; \
			OUT *q = (OUT *) IM_REGION_ADDR( or, le, y ) + z*2; \
			OUT *tlut = (OUT *) st->table[z]; \
			\
			for( x = 0; x < ne; x += b ) { \
				int n = p[x]*2; \
				\
				q[0] = tlut[n]; \
				q[1] = tlut[n + 1]; \
				q += b*2; \
			} \
		} \
	} \
}

#define loopg(IN,OUT) { \
	int b = st->nb; \
	\
	for( y = to; y < bo; y++ ) { \
		for( z = 0; z < b; z++ ) { \
			IN *p = (IN *) IM_REGION_ADDR( ir, le, y ); \
			OUT *q = (OUT *) IM_REGION_ADDR( or, le, y ); \
			OUT *tlut = (OUT *) st->table[z]; \
			\
			for( x = z; x < ne; x += b ) { \
				int index = p[x]; \
				\
				if( index > st->clp ) { \
					index = st->clp; \
					seq->overflow++; \
				} \
				\
				q[x] = tlut[index]; \
			} \
		} \
	} \
}

#define loopcg(IN,OUT) { \
	int b = in->Bands; \
	\
	for( y = to; y < bo; y++ ) { \
		for( z = 0; z < b; z++ ) { \
			IN *p = (IN *) IM_REGION_ADDR( ir, le, y ) + z; \
			OUT *q = (OUT *) IM_REGION_ADDR( or, le, y ) + z*2; \
			OUT *tlut = (OUT *) st->table[z]; \
			\
			for( x = 0; x < ne; x += b ) { \
				int index = p[x]; \
				\
				if( index > st->clp ) { \
					index = st->clp; \
					seq->overflow++; \
				} \
				\
				q[0] = tlut[ index*2 ]; \
				q[1] = tlut[ index*2+1 ]; \
				\
				q += b*2; \
			} \
		} \
	} \
}

/* Map image through one non-complex lut.
 */
#define loop1(OUT) { \
	OUT *tlut = (OUT *) st->table[0]; \
	\
	for( y = to; y < bo; y++ ) { \
		OUT *q = (OUT *) IM_REGION_ADDR( or, le, y ); \
		PEL *p = (PEL *) IM_REGION_ADDR( ir, le, y ); \
		\
		for( x = 0; x < ne; x++ ) \
			q[x] = tlut[p[x]]; \
	} \
}

/* Map image through one complex lut.
 */
#define loop1c(OUT) { \
	OUT *tlut = (OUT *) st->table[0]; \
	\
	for( y = to; y < bo; y++ ) { \
		OUT *q = (OUT *) IM_REGION_ADDR( or, le, y ); \
		PEL *p = (PEL *) IM_REGION_ADDR( ir, le, y ); \
		\
		for( x = 0; x < ne; x++ ) { \
			int n = p[x] * 2; \
			\
			q[0] = tlut[n]; \
			q[1] = tlut[n + 1]; \
			q += 2; \
		} \
	} \
}

/* As above, but the input image may be any unsigned integer type. We have to
 * index the lut carefully, and record the number of overflows we detect.
 */
#define loop1g(IN,OUT) { \
	OUT *tlut = (OUT *) st->table[0]; \
	\
	for( y = to; y < bo; y++ ) { \
		OUT *q = (OUT *) IM_REGION_ADDR( or, le, y ); \
		IN *p = (IN *) IM_REGION_ADDR( ir, le, y ); \
		\
		for( x = 0; x < ne; x++ ) { \
			int index = p[x]; \
			\
			if( index > st->clp ) { \
				index = st->clp; \
				seq->overflow++; \
			} \
			\
			q[x] = tlut[index]; \
		} \
	} \
}

#define loop1cg(IN,OUT) { \
	OUT *tlut = (OUT *) st->table[0]; \
	\
	for( y = to; y < bo; y++ ) { \
		OUT *q = (OUT *) IM_REGION_ADDR( or, le, y ); \
		IN *p = (IN *) IM_REGION_ADDR( ir, le, y ); \
		\
		for( x = 0; x < ne; x++ ) { \
			int index = p[x]; \
			\
			if( index > st->clp ) { \
				index = st->clp; \
				seq->overflow++; \
			} \
			\
			q[0] = tlut[index*2]; \
			q[1] = tlut[index*2 + 1]; \
			q += 2; \
		} \
	} \
}

/* Map 1-band image through a many-band non-complex lut.
 */
#define loop1m(OUT) { \
	OUT **tlut = (OUT **) st->table; \
	\
	for( y = to; y < bo; y++ ) { \
		OUT *q = (OUT *) IM_REGION_ADDR( or, le, y ); \
		PEL *p = (PEL *) IM_REGION_ADDR( ir, le, y ); \
		\
		for( i = 0, x = 0; x < np; x++ ) { \
			int n = p[x]; \
			\
			for( z = 0; z < st->nb; z++, i++ ) \
				q[i] = tlut[z][n]; \
		} \
	} \
}

/* Map 1-band image through many-band complex lut.
 */
#define loop1cm(OUT) { \
	OUT **tlut = (OUT **) st->table; \
	\
	for( y = to; y < bo; y++ ) { \
		OUT *q = (OUT *) IM_REGION_ADDR( or, le, y ); \
		PEL *p = (PEL *) IM_REGION_ADDR( ir, le, y ); \
		\
		for( x = 0; x < np; x++ ) { \
			int n = p[x] * 2; \
			\
			for( z = 0; z < st->nb; z++ ) { \
				q[0] = tlut[z][n]; \
				q[1] = tlut[z][n+1]; \
				q += 2; \
			} \
		} \
	} \
}

/* Map 1-band uint or ushort image through a many-band non-complex LUT.
 */
#define loop1gm(IN,OUT) { \
	OUT **tlut = (OUT **) st->table; \
	\
	for( y = to; y < bo; y++ ) { \
		IN *p = (IN *) IM_REGION_ADDR( ir, le, y ); \
		OUT *q = (OUT *) IM_REGION_ADDR( or, le, y ); \
		\
		for( i = 0, x = 0; x < np; x++ ) { \
			int n = p[x]; \
			\
			if( n > st->clp ) { \
				n = st->clp; \
				seq->overflow++; \
			} \
			\
			for( z = 0; z < st->nb; z++, i++ ) \
				q[i] = tlut[z][n]; \
		} \
	} \
}

/* Map 1-band uint or ushort image through a many-band complex LUT.
 */
#define loop1cgm(IN,OUT) { \
	OUT **tlut = (OUT **) st->table; \
	\
	for( y = to; y < bo; y++ ) { \
		IN *p = (IN *) IM_REGION_ADDR( ir, le, y ); \
		OUT *q = (OUT *) IM_REGION_ADDR( or, le, y ); \
		\
		for( x = 0; x < np; x++ ) { \
			int n = p[x]; \
			\
			if( n > st->clp ) { \
				n = st->clp; \
				seq->overflow++; \
			} \
			\
			for( z = 0; z < st->nb; z++ ) { \
				q[0] = tlut[z][n*2]; \
				q[1] = tlut[z][n*2 + 1]; \
				q += 2; \
			} \
		} \
	} \
}

/* Switch for input types. Has to be uint type!
 */
#define inner_switch( UCHAR, GEN, OUT ) \
	switch( ir->im->BandFmt ) { \
	case IM_BANDFMT_UCHAR:		UCHAR( OUT ); break; \
	case IM_BANDFMT_USHORT:		GEN( unsigned short, OUT ); break; \
	case IM_BANDFMT_UINT:		GEN( unsigned int, OUT ); break; \
	default: \
		im_error( "im_maplut", _( "bad input file" ) ); \
		return( -1 ); \
	}

/* Switch for LUT types. One function for non-complex images, a
 * variant for complex ones. Another pair as well, in case the input is not
 * uchar.
 */
#define outer_switch( UCHAR_F, UCHAR_FC, GEN_F, GEN_FC ) \
	switch( st->fmt ) { \
	case IM_BANDFMT_UCHAR:		inner_switch( UCHAR_F, GEN_F, \
					unsigned char ); break; \
	case IM_BANDFMT_CHAR:		inner_switch( UCHAR_F, GEN_F, \
					char ); break; \
	case IM_BANDFMT_USHORT:		inner_switch( UCHAR_F, GEN_F, \
					unsigned short ); break; \
	case IM_BANDFMT_SHORT:		inner_switch( UCHAR_F, GEN_F, \
					short ); break; \
	case IM_BANDFMT_UINT:		inner_switch( UCHAR_F, GEN_F, \
					unsigned int ); break; \
	case IM_BANDFMT_INT:		inner_switch( UCHAR_F, GEN_F, \
					int ); break; \
	case IM_BANDFMT_FLOAT:		inner_switch( UCHAR_F, GEN_F, \
					float ); break; \
	case IM_BANDFMT_DOUBLE:		inner_switch( UCHAR_F, GEN_F, \
					double ); break; \
	case IM_BANDFMT_COMPLEX:	inner_switch( UCHAR_FC, GEN_FC, \
					float ); break; \
	case IM_BANDFMT_DPCOMPLEX:	inner_switch( UCHAR_FC, GEN_FC, \
					double ); break; \
	default: \
		im_error( "im_maplut", _( "bad lut file" ) ); \
		return( -1 ); \
	}

/* Do a map.
 */
static int 
gen_map( REGION *or, Seq *seq, IMAGE *in, LutInfo *st )
{
	REGION *ir = seq->ir;
	Rect *r = &or->valid;
	int le = r->left;
	int to = r->top;
	int bo = IM_RECT_BOTTOM(r);
	int np = r->width;		/* Pels across region */
	int ne = IM_REGION_N_ELEMENTS( or );		/* Number of elements */
	int x, y, z, i;

	/* Get input ready.
	 */
	if( im_prepare( ir, r ) )
		return( -1 );

	/* Process!
	 */
	if( st->nb == 1 )
		/* One band lut.
		 */
		outer_switch( loop1, loop1c, loop1g, loop1cg ) 
	else 
		/* Many band lut.
		 */
		if( ir->im->Bands == 1 )
			/* ... but 1 band input.
			 */
			outer_switch( loop1m, loop1cm, loop1gm, loop1cgm ) 
		else
			outer_switch( loop, loopc, loopg, loopcg )

	return( 0 );
}

int 
im_maplut( IMAGE *in, IMAGE *out, IMAGE *lut )
{
	LutInfo *st;

	/* Check lut.
	 */
	if( lut->Coding != IM_CODING_NONE ) {
                im_error( "im_maplut", _( "lut is not uncoded" ) );
                return( -1 );
	}
	if( lut->Xsize * lut->Ysize > 100000 ) {
                im_error( "im_maplut", _( "lut seems very large!" ) );
                return( -1 );
	}

	/* Check input output. Old-style IO from lut, for simplicity.
	 */
	if( im_piocheck( in, out ) || im_incheck( lut ) )
		return( -1 );

	/* Check args.
	 */
        if( in->Coding != IM_CODING_NONE ) {
                im_error( "im_maplut", _( "input is not uncoded" ) );
                return( -1 );
	}
        if( !im_isuint( in ) ) {
                im_error( "im_maplut", _( "input is not some unsigned "
			"integer type" ) );
                return( -1 );
	}
	if( in->Bands != 1 && lut->Bands != 1 && lut->Bands != in->Bands ) {
                im_error( "im_maplut", _( "lut should have 1 band, or same "
			"number of bands as input, or any number of bands "
			"if input has 1 band" ) );
                return( -1 );
	}
	if( in->BandFmt == IM_BANDFMT_UCHAR && 
		lut->Xsize * lut->Ysize != 256 ) {
		im_error( "im_maplut", _( "input is uchar and lut does not "
			"have 256 elements" ) );
		return( -1 );
	}

	/* Prepare the output header.
	 */
        if( im_cp_descv( out, in, lut, NULL ) )
                return( -1 );

	/* Force output to be the same type as lut.
	 */
	out->Bbits = lut->Bbits;
	out->BandFmt = lut->BandFmt;

	/* Output has same number of bands as LUT, unless LUT has 1 band, in
	 * which case output has same number of bands as input.
	 */
	if( lut->Bands != 1 )
		out->Bands = lut->Bands;

	/* Make tables.
	 */
	if( !(st = build_luts( out, lut )) )
		return( -1 );

	/* Set demand hints.
	 */
	if( im_demand_hint( out, IM_THINSTRIP, in, NULL ) )
		return( -1 );

	/* Process!
	 */
        if( im_generate( out, start_maplut, gen_map, stop_maplut, in, st ) )
                return( -1 );

        return( 0 );
}