From 641c7fa430e2413d7843d37a79c96873f2b740ab Mon Sep 17 00:00:00 2001 From: John Cupitt Date: Thu, 20 Aug 2009 08:06:25 +0000 Subject: [PATCH] stuff --- ChangeLog | 3 +- TODO | 5 +- libvips/arithmetic/im_divide.c | 248 +++++++++++++++++++----------- libvips/arithmetic/im_multiply.c | 67 ++++---- libvips/arithmetic/im_remainder.c | 66 ++++---- 5 files changed, 236 insertions(+), 153 deletions(-) diff --git a/ChangeLog b/ChangeLog index cb1d52d5..2104e2a5 100644 --- a/ChangeLog +++ b/ChangeLog @@ -21,7 +21,8 @@ - add and use im_check_uncoded() and friends - matlab load handles column-major and plane-separated images (thanks Mikhail) - JPEG save allows "none" for profile, meaning don't attach a profile -- saner, simpler, faster typecasting for im_add(), im_subtract() +- saner, simpler, faster typecasting for im_add(), im_subtract(), + im_multiply(), im_divide() - im_measure() allows sel == NULL, meaning all patches 25/3/09 started 7.18.0 diff --git a/TODO b/TODO index 3ec0ac03..d84ea669 100644 --- a/TODO +++ b/TODO @@ -1,4 +1,4 @@ - +- revising im_remainder() - 1-bit PNG read is broken? @@ -11,9 +11,6 @@ -- im_multiply() should do proper complex multiplication, since im_divide does - - - im__cast_and_call() no longer does out->Bbits = im_bits_of_fmt( out->BandFmt ); diff --git a/libvips/arithmetic/im_divide.c b/libvips/arithmetic/im_divide.c index be9fcc23..bdd9e189 100644 --- a/libvips/arithmetic/im_divide.c +++ b/libvips/arithmetic/im_divide.c @@ -1,13 +1,4 @@ -/* @(#) Divide two images - * @(#) Images must have the same no of bands and can be of any type - * @(#) No check for overflow is carried out. - * @(#) - * @(#) int - * @(#) im_divide(in1, in2, out) - * @(#) IMAGE *in1, *in2, *out; - * @(#) - * @(#) Returns 0 on success and -1 on error - * @(#) +/* im_divide.c * * Copyright: 1990, N. Dessipris. * @@ -29,6 +20,9 @@ * - updated for 1 band $op n band image -> n band image case * 8/12/06 * - add liboil support + * 18/8/08 + * - revise upcasting system + * - add gtkdoc comments */ /* @@ -83,57 +77,66 @@ #ifdef USE_MODARG_DIV /* This is going to be much slower */ -#define cloop(TYPE) \ -{ \ - TYPE *X= (TYPE*) in[0]; \ - TYPE *Y= (TYPE*) in[1]; \ - TYPE *Z= (TYPE*) out; \ - TYPE *Z_stop= Z + sz * 2; \ - \ - for( ; Z < Z_stop; X+= 2, Y+=2 ){ \ - double arg= atan2( X[1], X[0] ) - atan2( Y[1], Y[0] ); \ - double mod= hypot( X[1], X[0] ) / hypot( Y[1], Y[0] ); \ - *Z++= mod * cos( arg ); \ - *Z++= mod * sin( arg ); \ - } \ +#define CLOOP( TYPE ) { \ + TYPE *X = (TYPE *) in[0]; \ + TYPE *Y = (TYPE *) in[1]; \ + TYPE *Z = (TYPE *) out; \ + int i; \ + \ + for( i = 0; i < sz; i++ ) { \ + double arg = atan2( X[1], X[0] ) - atan2( Y[1], Y[0] ); \ + double mod = hypot( X[1], X[0] ) / hypot( Y[1], Y[0] ); \ + \ + Z[0] = mod * cos( arg ); \ + Z[1] = mod * sin( arg ); \ + \ + X += 2; \ + Y += 2; \ + Z += 2; \ + } \ } #else /* USE_MODARG_DIV */ -#define cloop(TYPE) \ -{ \ - TYPE *X= (TYPE*) in[0]; \ - TYPE *Y= (TYPE*) in[1]; \ - TYPE *Z= (TYPE*) out; \ - TYPE *Z_stop= Z + sz * 2; \ - \ - for( ; Z < Z_stop; X+= 2, Y+=2 ) \ - if( fabs( Y[0] ) > fabs( Y[1] )){ \ - double a= Y[1] / Y[0]; \ - double b= Y[0] + Y[1] * a; \ - *Z++= ( X[0] + X[1] * a ) / b; \ - *Z++= ( X[1] - X[0] * a ) / b; \ - } \ - else { \ - double a= Y[0] / Y[1]; \ - double b= Y[1] + Y[0] * a; \ - *Z++= ( X[0] * a + X[1] ) / b; \ - *Z++= ( X[1] * a - X[0] ) / b; \ - } \ +#define CLOOP( TYPE ) { \ + TYPE *X = (TYPE *) in[0]; \ + TYPE *Y = (TYPE *) in[1]; \ + TYPE *Z = (TYPE *) out; \ + int i; \ + \ + for( i = 0; i < sz; i++ ) { \ + if( fabs( Y[0] ) > fabs( Y[1] ) ) { \ + double a = Y[1] / Y[0]; \ + double b = Y[0] + Y[1] * a; \ + \ + Z[0] = (X[0] + X[1] * a) / b; \ + Z[1] = (X[1] - X[0] * a) / b; \ + } \ + else { \ + double a = Y[0] / Y[1]; \ + double b = Y[1] + Y[0] * a; \ + \ + Z[0] = (X[0] * a + X[1]) / b; \ + Z[1] = (X[1] * a - X[0]) / b; \ + } \ + \ + X += 2; \ + Y += 2; \ + Z += 2; \ + } \ } #endif /* USE_MODARG_DIV */ /* Real divide. */ -#define rloop(TYPE) \ -{\ - TYPE *p1 = (TYPE *) in[0];\ - TYPE *p2 = (TYPE *) in[1];\ - TYPE *q = (TYPE *) out;\ +#define RLOOP( 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];\ + for( x = 0; x < sz; x++ ) \ + q[x] = p1[x] / p2[x]; \ } static void @@ -145,31 +148,126 @@ divide_buffer( PEL **in, PEL *out, int width, IMAGE *im ) /* Divide all input types. */ switch( im->BandFmt ) { - case IM_BANDFMT_CHAR: rloop( signed char ); break; - case IM_BANDFMT_UCHAR: rloop( unsigned char ); break; - case IM_BANDFMT_SHORT: rloop( signed short ); break; - case IM_BANDFMT_USHORT: rloop( unsigned short ); break; - case IM_BANDFMT_INT: rloop( signed int ); break; - case IM_BANDFMT_UINT: rloop( unsigned int ); break; + case IM_BANDFMT_CHAR: RLOOP( signed char, float ); break; + case IM_BANDFMT_UCHAR: RLOOP( unsigned char, float ); break; + case IM_BANDFMT_SHORT: RLOOP( signed short, float ); break; + case IM_BANDFMT_USHORT: RLOOP( unsigned short, float ); break; + case IM_BANDFMT_INT: RLOOP( signed int, float ); break; + case IM_BANDFMT_UINT: RLOOP( unsigned int, float ); break; case IM_BANDFMT_FLOAT: #ifdef HAVE_LIBOIL oil_divide_f32( (float *) out, (float *) in[0], (float *) in[1], sz ); #else /*!HAVE_LIBOIL*/ - rloop( float ); + RLOOP( float, float ); #endif /*HAVE_LIBOIL*/ break; - case IM_BANDFMT_DOUBLE: rloop( double ); break; - case IM_BANDFMT_COMPLEX: cloop( float ); break; - case IM_BANDFMT_DPCOMPLEX: cloop( double ); break; + case IM_BANDFMT_DOUBLE: RLOOP( double, double ); break; + case IM_BANDFMT_COMPLEX: CLOOP( float ); break; + case IM_BANDFMT_DPCOMPLEX: CLOOP( double ); break; default: assert( 0 ); } } +/* Save a bit of typing. + */ +#define F IM_BANDFMT_FLOAT +#define X IM_BANDFMT_COMPLEX +#define D IM_BANDFMT_DOUBLE +#define DX IM_BANDFMT_DPCOMPLEX + +/* Type promotion for division. Sign and value preserving. Make sure + * these match the case statement in divide_buffer() above. + */ +static int bandfmt_divide[10] = { +/* UC C US S UI I F X D DX */ + F, F, F, F, F, F, F, X, D, DX +}; + +/** + * im_divide: + * @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_divide() type promotion + * + * + * + * input type + * output type + * + * + * + * + * uchar + * float + * + * + * char + * float + * + * + * ushort + * float + * + * + * short + * float + * + * + * uint + * float + * + * + * int + * float + * + * + * 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_multiply(), im_lintra(). + * + * Returns: 0 on success, -1 on error + */ int im_divide( IMAGE *in1, IMAGE *in2, IMAGE *out ) { @@ -187,34 +285,10 @@ im_divide( IMAGE *in1, IMAGE *in2, IMAGE *out ) */ out->Bands = IM_MAX( in1->Bands, in2->Bands ); - /* What output type will we write? float, double or complex. + /* What output type will we write? int, float or complex. */ - 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. - */ - out->BandFmt = IM_BANDFMT_DPCOMPLEX; - else - out->BandFmt = 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 ) - out->BandFmt = IM_BANDFMT_DOUBLE; - else - out->BandFmt = IM_BANDFMT_FLOAT; - } - else { - /* An int type -- output must be just float. - */ - out->BandFmt = IM_BANDFMT_FLOAT; - } + out->BandFmt = bandfmt_divide[im__format_common( in1, in2 )]; + out->Bbits = im_bits_of_fmt( out->BandFmt ); /* And process! */ diff --git a/libvips/arithmetic/im_multiply.c b/libvips/arithmetic/im_multiply.c index 1b94eaef..23f55200 100644 --- a/libvips/arithmetic/im_multiply.c +++ b/libvips/arithmetic/im_multiply.c @@ -1,13 +1,4 @@ -/* @(#) Multiply two images - * @(#) Images must have the same no of bands and can be of any type - * @(#) No check for overflow is carried out. - * @(#) - * @(#) int - * @(#) im_multiply(in1, in2, out) - * @(#) IMAGE *in1, *in2, *out; - * @(#) - * @(#) Returns 0 on success and -1 on error - * @(#) +/* im_multiply.c * * Copyright: 1990, N. Dessipris. * @@ -32,7 +23,6 @@ * 18/8/08 * - revise upcasting system * - add gtkdoc comments - * - remove separate complex case, just double size */ /* @@ -82,7 +72,32 @@ #include #endif /*WITH_DMALLOC*/ -#define LOOP( IN, OUT ) { \ +/* Complex multiply. + */ +#define CLOOP( TYPE ) { \ + TYPE *p1 = (TYPE *) in[0]; \ + TYPE *p2 = (TYPE *) in[1]; \ + TYPE *q = (TYPE *) out; \ + \ + for( x = 0; x < sz; x++ ) { \ + double x1 = p1[0]; \ + double y1 = p1[1]; \ + double x2 = p2[0]; \ + double y2 = p2[1]; \ + \ + p1 += 2; \ + p2 += 2; \ + \ + q[0] = x1 * x2 - y1 * y2; \ + q[1] = x1 * y2 + x2 * y1; \ + \ + q += 2; \ + } \ +} + +/* Real multiply. + */ +#define RLOOP( IN, OUT ) { \ IN *p1 = (IN *) in[0]; \ IN *p2 = (IN *) in[1]; \ OUT *q = (OUT *) out; \ @@ -94,9 +109,7 @@ static void multiply_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); + const int sz = width * im->Bands; int x; @@ -104,27 +117,25 @@ multiply_buffer( PEL **in, PEL *out, int width, IMAGE *im ) * bandfmt_multiply[] below. */ switch( im->BandFmt ) { - case IM_BANDFMT_CHAR: LOOP( signed char, signed short ); break; - case IM_BANDFMT_UCHAR: LOOP( unsigned char, unsigned short ); break; - case IM_BANDFMT_SHORT: LOOP( signed short, signed int ); break; - case IM_BANDFMT_USHORT: LOOP( unsigned short, unsigned int ); break; - case IM_BANDFMT_INT: LOOP( signed int, signed int ); break; - case IM_BANDFMT_UINT: LOOP( unsigned int, unsigned int ); break; + case IM_BANDFMT_CHAR: RLOOP( signed char, signed short ); break; + case IM_BANDFMT_UCHAR: RLOOP( unsigned char, unsigned short ); break; + case IM_BANDFMT_SHORT: RLOOP( signed short, signed int ); break; + case IM_BANDFMT_USHORT: RLOOP( unsigned short, unsigned int ); break; + case IM_BANDFMT_INT: RLOOP( signed int, signed int ); break; + case IM_BANDFMT_UINT: RLOOP( unsigned int, unsigned int ); break; case IM_BANDFMT_FLOAT: - case IM_BANDFMT_COMPLEX: #ifdef HAVE_LIBOIL oil_multiply_f32( (float *) out, (float *) in[0], (float *) in[1], sz ); #else /*!HAVE_LIBOIL*/ - LOOP( float, float ); + RLOOP( float, float ); #endif /*HAVE_LIBOIL*/ break; - case IM_BANDFMT_DOUBLE: - case IM_BANDFMT_DPCOMPLEX: - LOOP( double, double ); - break; + case IM_BANDFMT_COMPLEX: CLOOP( float ); break; + case IM_BANDFMT_DOUBLE: RLOOP( double, double ); break; + case IM_BANDFMT_DPCOMPLEX: CLOOP( double ); break; default: assert( 0 ); @@ -172,7 +183,7 @@ static int bandfmt_multiply[10] = { * following table is used to determine the output type: * * - * im_add() type promotion + * im_multiply() type promotion * * * diff --git a/libvips/arithmetic/im_remainder.c b/libvips/arithmetic/im_remainder.c index 8f6a0479..a4e63062 100644 --- a/libvips/arithmetic/im_remainder.c +++ b/libvips/arithmetic/im_remainder.c @@ -1,4 +1,4 @@ -/* @(#) Remainder after integer division +/* im_remainder.c * * 2/8/99 JC * - im_divide adapted to make im_remainder @@ -56,16 +56,16 @@ #include #endif /*WITH_DMALLOC*/ -#define loop(TYPE) {\ - TYPE *p1 = (TYPE *) in[0];\ - TYPE *p2 = (TYPE *) in[1];\ - TYPE *q = (TYPE *) out;\ +#define LOOP(TYPE) { \ + TYPE *p1 = (TYPE *) in[0]; \ + TYPE *p2 = (TYPE *) in[1]; \ + TYPE *q = (TYPE *) out; \ \ - for( x = 0; x < sz; x++ )\ - if( p2[x] )\ - q[x] = p1[x] % p2[x];\ - else\ - q[x] = -1;\ + for( x = 0; x < sz; x++ ) \ + if( p2[x] ) \ + q[x] = p1[x] % p2[x]; \ + else \ + q[x] = -1; \ } static void @@ -75,39 +75,38 @@ remainder_buffer( PEL **in, PEL *out, int width, IMAGE *im ) int sz = width * im->Bands; switch( im->BandFmt ) { - case IM_BANDFMT_CHAR: loop( signed char ); break; - case IM_BANDFMT_UCHAR: loop( unsigned char ); break; - case IM_BANDFMT_SHORT: loop( signed short ); break; - case IM_BANDFMT_USHORT: loop( unsigned short ); break; - case IM_BANDFMT_INT: loop( signed int ); break; - case IM_BANDFMT_UINT: loop( unsigned int ); break; + case IM_BANDFMT_CHAR: LOOP( signed char ); break; + case IM_BANDFMT_UCHAR: LOOP( unsigned char ); break; + case IM_BANDFMT_SHORT: LOOP( signed short ); break; + case IM_BANDFMT_USHORT: LOOP( unsigned short ); break; + case IM_BANDFMT_INT: LOOP( signed int ); break; + case IM_BANDFMT_UINT: LOOP( unsigned int ); break; default: assert( 0 ); } } +/* +.B im_remainder(3) +calculates the remainder after integer division of two images. The output +type is the same as the type of +.B in1 +unless +.B in1 +is float or complex, in which +case the output type is signed integer. + */ int im_remainder( IMAGE *in1, IMAGE *in2, IMAGE *out ) { - /* Basic checks. - */ - if( im_piocheck( in1, out ) || im_pincheck( in2 ) ) + if( im_piocheck( in1, out ) || + im_pincheck( in2 ) || + im_check_bands_1orn( "im_remainder", in1, in2 ) || + im_check_uncoded( "im_remainder", in1 ) || + im_check_uncoded( "im_remainder", in2 ) ) return( -1 ); - if( in1->Xsize != in2->Xsize || in1->Ysize != in2->Ysize ) { - im_error( "im_remainder", "%s", _( "not same size" ) ); - return( -1 ); - } - if( in1->Bands != in2->Bands && - (in1->Bands != 1 && in2->Bands != 1) ) { - im_error( "im_remainder", - "%s", _( "not same number of bands" ) ); - return( -1 ); - } - if( in1->Coding != IM_CODING_NONE || in2->Coding != IM_CODING_NONE ) { - im_error( "im_remainder", "%s", _( "not uncoded" ) ); - return( -1 ); - } + if( im_cp_descv( out, in1, in2, NULL ) ) return( -1 ); @@ -120,6 +119,7 @@ im_remainder( IMAGE *in1, IMAGE *in2, IMAGE *out ) */ if( im_isfloat( in1 ) || im_iscomplex( in1 ) ) out->BandFmt = IM_BANDFMT_INT; + out->Bbits = im_bits_of_fmt( out->BandFmt ); /* And process! */