From 5b4c56c53f90f6754cbd5644fbc67b0b607a6056 Mon Sep 17 00:00:00 2001 From: John Cupitt Date: Sat, 19 Sep 2009 10:36:34 +0000 Subject: [PATCH] stuff --- libvips/arithmetic/maths.c | 371 +++++++++++++++++++++++++++++++++++++ 1 file changed, 371 insertions(+) create mode 100644 libvips/arithmetic/maths.c diff --git a/libvips/arithmetic/maths.c b/libvips/arithmetic/maths.c new file mode 100644 index 00000000..8e202ef1 --- /dev/null +++ b/libvips/arithmetic/maths.c @@ -0,0 +1,371 @@ +/* maths.c - +, -, *, / of two images + * + * Copyright: 1990, N. Dessipris. + * + * Author: Nicos Dessipris + * Written on: 02/05/1990 + * Modified on: + * 29/4/93 J.Cupitt + * - now works for partial images + * 1/7/93 JC + * - adapted for partial v2 + * 9/5/95 JC + * - simplified: now just handles 10 cases (instead of 50), using + * im_clip2*() to help + * - now uses im_wrapmany() rather than im_generate() + * 31/5/96 JC + * - SWAP() removed, *p++ removed + * 27/9/04 + * - im__cast_and_call() now matches bands as well + * - ... so 1 band + 4 band image -> 4 band image + * 8/12/06 + * - add liboil support + * 18/8/08 + * - revise upcasting system + * - im__cast_and_call() no longer sets bbits for you + * - add gtkdoc comments + * - remove separate complex case, just double size + * 11/9/09 + * - im__cast_and_call() becomes im__arith_binary() + * - more of operation scaffold moved inside + * 18/9/09 + * - redne as math.c, now generates -, *, / as well + */ + +/* + + 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 +#endif /*HAVE_CONFIG_H*/ +#include + +#include +#include +#include +#include + +#include +#include + +#ifdef HAVE_LIBOIL +#include +#endif /*HAVE_LIBOIL*/ + +#ifdef WITH_DMALLOC +#include +#endif /*WITH_DMALLOC*/ + +#define LOOP( IN, OUT ) { \ + IN *p1 = (IN *) in[0]; \ + IN *p2 = (IN *) in[1]; \ + OUT *q = (OUT *) out; \ + \ + for( x = 0; x < sz; x++ ) \ + q[x] = p1[x] + p2[x]; \ +} + +static void +add_buffer( PEL **in, PEL *out, int width, IMAGE *im ) +{ + /* Complex just doubles the size. + */ + const int sz = width * im->Bands * (im_iscomplex( im ) ? 2 : 1); + + int x; + + /* Add all input types. Keep types here in sync with bandfmt_add[] + * below. + */ + switch( im->BandFmt ) { + case IM_BANDFMT_UCHAR: LOOP( unsigned char, unsigned short ); break; + case IM_BANDFMT_CHAR: LOOP( signed char, signed short ); break; + case IM_BANDFMT_USHORT: LOOP( unsigned short, unsigned int ); break; + case IM_BANDFMT_SHORT: LOOP( signed short, signed int ); break; + case IM_BANDFMT_UINT: LOOP( unsigned int, unsigned int ); break; + case IM_BANDFMT_INT: LOOP( signed int, signed int ); break; + + case IM_BANDFMT_FLOAT: + case IM_BANDFMT_COMPLEX: +#ifdef HAVE_LIBOIL + oil_add_f32( (float *) out, + (float *) in[0], (float *) in[1], sz ); +#else /*!HAVE_LIBOIL*/ + LOOP( float, float ); +#endif /*HAVE_LIBOIL*/ + break; + + case IM_BANDFMT_DOUBLE: + case IM_BANDFMT_DPCOMPLEX: + LOOP( double, double ); + break; + + default: + assert( 0 ); + } +} + +/* Save a bit of typing. + */ +#define UC IM_BANDFMT_UCHAR +#define C IM_BANDFMT_CHAR +#define US IM_BANDFMT_USHORT +#define S IM_BANDFMT_SHORT +#define UI IM_BANDFMT_UINT +#define I IM_BANDFMT_INT +#define F IM_BANDFMT_FLOAT +#define X IM_BANDFMT_COMPLEX +#define D IM_BANDFMT_DOUBLE +#define DX IM_BANDFMT_DPCOMPLEX + +/* For two integer types, the "largest", ie. one which can represent the + * full range of both. + */ +static int bandfmt_largest[6][6] = { + /* UC C US S UI I */ +/* UC */ { UC, S, US, S, UI, I }, +/* C */ { S, C, I, S, I, I }, +/* US */ { US, I, US, I, UI, I }, +/* S */ { S, S, I, S, I, I }, +/* UI */ { UI, I, UI, I, UI, I }, +/* I */ { I, I, I, I, I, I } +}; + +/* For two formats, find one which can represent the full range of both. + */ +static VipsBandFmt +im__format_common( IMAGE *in1, IMAGE *in2 ) +{ + if( im_iscomplex( in1 ) || im_iscomplex( in2 ) ) { + /* What kind of complex? + */ + if( in1->BandFmt == IM_BANDFMT_DPCOMPLEX || + in2->BandFmt == IM_BANDFMT_DPCOMPLEX ) + /* Output will be DPCOMPLEX. + */ + return( IM_BANDFMT_DPCOMPLEX ); + else + return( IM_BANDFMT_COMPLEX ); + + } + else if( im_isfloat( in1 ) || im_isfloat( in2 ) ) { + /* What kind of float? + */ + if( in1->BandFmt == IM_BANDFMT_DOUBLE || + in2->BandFmt == IM_BANDFMT_DOUBLE ) + return( IM_BANDFMT_DOUBLE ); + else + return( IM_BANDFMT_FLOAT ); + } + else + /* Must be int+int -> int. + */ + return( bandfmt_largest[in1->BandFmt][in2->BandFmt] ); +} + +/* Make an n-band image. Input 1 or n bands. + */ +int +im__bandup( IMAGE *in, IMAGE *out, int n ) +{ + IMAGE *bands[256]; + int i; + + if( in->Bands == n ) + return( im_copy( in, out ) ); + if( in->Bands != 1 ) { + im_error( "im__bandup", _( "not one band or %d bands" ), n ); + return( -1 ); + } + if( n > 256 || n < 1 ) { + im_error( "im__bandup", "%s", _( "bad bands" ) ); + return( -1 ); + } + + for( i = 0; i < n; i++ ) + bands[i] = in; + + return( im_gbandjoin( bands, out, n ) ); +} + +/* The common part of most binary arithmetic, relational and boolean + * operators. We: + * + * - check in and out + * - cast in1 and in2 up to a common format + * - cast the common format to the output format with the supplied table + * - equalise bands + * - run the supplied buffer operation + */ +int +im__arith_binary( const char *name, + IMAGE *in1, IMAGE *in2, IMAGE *out, + int format_table[10], + im_wrapmany_fn fn, void *a ) +{ + VipsBandFmt fmt; + IMAGE *t[5]; + + if( im_piocheck( in1, out ) || + im_pincheck( in2 ) || + im_check_bands_1orn( name, in1, in2 ) || + im_check_uncoded( name, in1 ) || + im_check_uncoded( name, in2 ) ) + return( -1 ); + + if( im_cp_descv( out, in1, in2, NULL ) ) + return( -1 ); + + /* What number of bands will we write? + */ + out->Bands = IM_MAX( in1->Bands, in2->Bands ); + + /* What output type will we write? int, float or complex. + */ + out->BandFmt = format_table[im__format_common( in1, in2 )]; + out->Bbits = im_bits_of_fmt( out->BandFmt ); + + if( im_open_local_array( out, t, 4, "type cast:1", "p" ) ) + return( -1 ); + + /* Cast our input images up to a common type. + */ + fmt = im__format_common( in1, in2 ); + if( im_clip2fmt( in1, t[0], fmt ) || + im_clip2fmt( in2, t[1], fmt ) ) + return( -1 ); + + /* Force bands up to the same as out. + */ + if( im__bandup( t[0], t[2], out->Bands ) || + im__bandup( t[1], t[3], out->Bands ) ) + return( -1 ); + + /* And process! + */ + t[4] = NULL; + if( im_wrapmany( t + 2, out, fn, out, a ) ) + return( -1 ); + + return( 0 ); +} + +/* Type promotion for addition. Sign and value preserving. Make sure these + * match the case statement in add_buffer() above. + */ +static int bandfmt_add[10] = { +/* UC C US S UI I F X D DX */ + US, S, UI, I, UI, I, F, X, D, DX +}; + +/** + * im_add: + * @in1: input #IMAGE 1 + * @in2: input #IMAGE 2 + * @out: output #IMAGE + * + * This operation calculates @in1 + @in2 and writes the result to @out. + * The images must be the same size. They may have any format. + * + * If the number of bands differs, one of the images + * must have one band. In this case, an n-band image is formed from the + * one-band image by joining n copies of the one-band image together, and then + * the two n-band images are operated upon. + * + * The two input images are cast up to the smallest common type (see table + * Smallest common format in + * arithmetic), then the + * following table is used to determine the output type: + * + * + * im_add() type promotion + * + * + * + * input type + * output type + * + * + * + * + * uchar + * ushort + * + * + * char + * short + * + * + * ushort + * uint + * + * + * short + * int + * + * + * uint + * uint + * + * + * int + * int + * + * + * float + * float + * + * + * double + * double + * + * + * complex + * complex + * + * + * double complex + * double complex + * + * + * + *
+ * + * In other words, the output type is just large enough to hold the whole + * range of possible values. + * + * See also: im_subtract(), im_lintra(). + * + * Returns: 0 on success, -1 on error + */ +int +im_add( IMAGE *in1, IMAGE *in2, IMAGE *out ) +{ + return( im__arith_binary( "im_add", + in1, in2, out, + bandfmt_add, + (im_wrapmany_fn) add_buffer, NULL ) ); +}