From 1de6c6bcc6b52d6418df4b742ed2d882a4431ae9 Mon Sep 17 00:00:00 2001 From: John Cupitt Date: Thu, 6 Dec 2012 14:52:51 +0000 Subject: [PATCH] add binary complex ops cross_phase is the only one atm --- ChangeLog | 1 + TODO | 2 +- libvips/arithmetic/Makefile.am | 1 - libvips/arithmetic/arithmetic.c | 2 + libvips/arithmetic/complex.c | 239 ++++++++++++++++++++++++++++ libvips/arithmetic/im_cross_phase.c | 150 ----------------- libvips/arithmetic/math2.c | 2 +- libvips/include/vips/arithmetic.h | 17 ++ libvips/include/vips/enumtypes.h | 2 + libvips/include/vips/inlines.h | 40 ++--- libvips/iofuncs/enumtypes.c | 17 ++ 11 files changed, 301 insertions(+), 172 deletions(-) delete mode 100644 libvips/arithmetic/im_cross_phase.c diff --git a/ChangeLog b/ChangeLog index 7ed54f1c..b6534538 100644 --- a/ChangeLog +++ b/ChangeLog @@ -28,6 +28,7 @@ - deprecate im_maxpos_avg(): too specialised to be worth maintaining - deprecate im_linreg(): easily done by combining other operators - deprecate im_point(): easily done by combining other operators +- add binary complex operations, with cross_phase as the only one so far 14/11/12 started 7.30.6 - capture tiff warnings earlier diff --git a/TODO b/TODO index 592735d3..790879cb 100644 --- a/TODO +++ b/TODO @@ -1,4 +1,4 @@ -- im_cross_phase.c needs a rewrite? +- im_cross_phase.c needs a vips7 compat wrapper - add vips_band()/vips_bor() diff --git a/libvips/arithmetic/Makefile.am b/libvips/arithmetic/Makefile.am index 84ac1c05..c780f262 100644 --- a/libvips/arithmetic/Makefile.am +++ b/libvips/arithmetic/Makefile.am @@ -3,7 +3,6 @@ noinst_LTLIBRARIES = libarithmetic.la libarithmetic_la_SOURCES = \ abs.c \ complex.c \ - im_cross_phase.c \ deviate.c \ divide.c \ measure.c \ diff --git a/libvips/arithmetic/arithmetic.c b/libvips/arithmetic/arithmetic.c index 9ebf9d84..350d2eb7 100644 --- a/libvips/arithmetic/arithmetic.c +++ b/libvips/arithmetic/arithmetic.c @@ -700,6 +700,7 @@ vips_arithmetic_operation_init( void ) extern GType vips_math2_get_type( void ); extern GType vips_math2_const_get_type( void ); extern GType vips_complex_get_type( void ); + extern GType vips_complex2_get_type( void ); extern GType vips_complexget_get_type( void ); extern GType vips_complexform_get_type( void ); @@ -728,6 +729,7 @@ vips_arithmetic_operation_init( void ) vips_math2_get_type(); vips_math2_const_get_type(); vips_complex_get_type(); + vips_complex2_get_type(); vips_complexget_get_type(); vips_complexform_get_type(); } diff --git a/libvips/arithmetic/complex.c b/libvips/arithmetic/complex.c index 89c202f3..797ce1ef 100644 --- a/libvips/arithmetic/complex.c +++ b/libvips/arithmetic/complex.c @@ -334,6 +334,245 @@ vips_conj( VipsImage *in, VipsImage **out, ... ) return( result ); } +typedef struct _VipsComplex2 { + VipsUnary parent_instance; + + VipsOperationComplex2 cmplx; + +} VipsComplex2; + +typedef VipsUnaryClass VipsComplex2Class; + +G_DEFINE_TYPE( VipsComplex2, vips_complex2, VIPS_TYPE_BINARY ); + +#define LOOP2( IN, OUT, OP ) { \ + IN *p1 = (IN *) in[0]; \ + IN *p2 = (IN *) in[1]; \ + OUT *q = (OUT *) out; \ + \ + for( x = 0; x < sz; x++ ) { \ + OP( q, p1[x], 0.0, p2[x], 0.0 ); \ + \ + q += 2; \ + } \ +} + +#define CLOOP2( IN, OUT, OP ) { \ + IN *p1 = (IN *) in[0]; \ + IN *p2 = (IN *) in[1]; \ + OUT *q = (OUT *) out; \ + \ + for( x = 0; x < sz; x++ ) { \ + OP( q, p1[0], p1[1], p2[0], p2[1] ); \ + \ + p1 += 2; \ + p2 += 2; \ + q += 2; \ + } \ +} + +#define SWITCH2( OP ) \ + switch( vips_image_get_format( im ) ) { \ + case VIPS_FORMAT_UCHAR: \ + LOOP2( unsigned char, float, OP ); break; \ + case VIPS_FORMAT_CHAR: \ + LOOP2( signed char, float, OP ); break; \ + case VIPS_FORMAT_USHORT: \ + LOOP2( unsigned short, float, OP ); break; \ + case VIPS_FORMAT_SHORT: \ + LOOP2( signed short, float, OP ); break; \ + case VIPS_FORMAT_UINT: \ + LOOP2( unsigned int, float, OP ); break; \ + case VIPS_FORMAT_INT: \ + LOOP2( signed int, float, OP ); break; \ + case VIPS_FORMAT_FLOAT: \ + LOOP2( float, float, OP ); break; \ + case VIPS_FORMAT_DOUBLE: \ + LOOP2( double, double, OP ); break;\ + case VIPS_FORMAT_COMPLEX: \ + CLOOP2( float, float, OP ); break; \ + case VIPS_FORMAT_DPCOMPLEX: \ + CLOOP2( double, double, OP ); break;\ + \ + default: \ + g_assert( 0 ); \ + } + +/* There doesn't seem to be much difference in speed between these two methods + * (on an Athlon64), so I use the modulus argument version, since atan2() is + * in c89 but hypot() is c99. + * + * If you think that it might be faster on your platform, uncomment the + * following: + */ +#define USE_MODARG_DIV + +#ifdef USE_MODARG_DIV + +#define CROSS( Q, X1, Y1, X2, Y2 ) { \ + double arg = atan2( X2, X1 ) - atan2( Y2, Y1 ); \ + \ + Q[0] = cos( arg ); \ + Q[1] = sin( arg ); \ +} + +#else /* USE_MODARG_DIV */ + +#define CROSS( Q, X1, Y1, X2, Y2 ) { \ + if( ABS( Y1 ) > ABS( Y2 ) ) { \ + double a = Y2 / Y1; \ + double b = Y1 + Y2 * a; \ + double re = (X1 + X2 * a) / b; \ + double im = (X2 - X1 * a) / b; \ + double mod = vips__hypot( re, im ); \ + \ + Q[0] = re / mod; \ + Q[1] = im / mod; \ + } \ + else { \ + double a = Y1 / Y2; \ + double b = Y2 + Y1 * a; \ + double re = (X1 * a + X2) / b; \ + double im = (X2 * a - X1) / b; \ + double mod = vips__hypot( re, im ); \ + \ + Q[0] = re / mod; \ + Q[1] = im / mod; \ + } \ +} + +#endif /* USE_MODARG_DIV */ + +static void +vips_complex2_buffer( VipsArithmetic *arithmetic, + VipsPel *out, VipsPel **in, int width ) +{ + VipsComplex2 *cmplx = (VipsComplex2 *) arithmetic; + VipsImage *im = arithmetic->ready[0]; + const int sz = width * vips_image_get_bands( im ); + + int x; + + switch( cmplx->cmplx ) { + case VIPS_OPERATION_COMPLEX2_CROSS_PHASE: SWITCH2( CROSS ); break; + + default: + g_assert( 0 ); + } +} + +/* Save a bit of typing. + */ +#define UC VIPS_FORMAT_UCHAR +#define C VIPS_FORMAT_CHAR +#define US VIPS_FORMAT_USHORT +#define S VIPS_FORMAT_SHORT +#define UI VIPS_FORMAT_UINT +#define I VIPS_FORMAT_INT +#define F VIPS_FORMAT_FLOAT +#define X VIPS_FORMAT_COMPLEX +#define D VIPS_FORMAT_DOUBLE +#define DX VIPS_FORMAT_DPCOMPLEX + +static const VipsBandFormat vips_bandfmt_complex2[10] = { +/* UC C US S UI I F X D DX */ + X, X, X, X, X, X, X, X, DX, DX +}; + +static void +vips_complex2_class_init( VipsComplex2Class *class ) +{ + GObjectClass *gobject_class = G_OBJECT_CLASS( class ); + VipsObjectClass *object_class = (VipsObjectClass *) class; + VipsArithmeticClass *aclass = VIPS_ARITHMETIC_CLASS( class ); + + gobject_class->set_property = vips_object_set_property; + gobject_class->get_property = vips_object_get_property; + + object_class->nickname = "complex2"; + object_class->description = + _( "perform a binary complex operation on two images" ); + + vips_arithmetic_set_format_table( aclass, vips_bandfmt_complex2 ); + + aclass->process_line = vips_complex2_buffer; + + VIPS_ARG_ENUM( class, "cmplx", 200, + _( "Operation" ), + _( "binary complex operation to perform" ), + VIPS_ARGUMENT_REQUIRED_INPUT, + G_STRUCT_OFFSET( VipsComplex2, cmplx ), + VIPS_TYPE_OPERATION_COMPLEX, + VIPS_OPERATION_COMPLEX2_CROSS_PHASE ); +} + +static void +vips_complex2_init( VipsComplex2 *cmplx ) +{ +} + +static int +vips_complex2v( VipsImage *left, VipsImage *right, VipsImage **out, + VipsOperationComplex2 cmplx, va_list ap ) +{ + return( vips_call_split( "complex2", ap, left, right, out, cmplx ) ); +} + +/** + * vips_complex2: + * @left: input #VipsImage + * @right: input #VipsImage + * @out: output #VipsImage + * @cmplx: complex2 operation to perform + * @...: %NULL-terminated list of optional named arguments + * + * Perform various binary operations on complex images. + * + * Angles are expressed in degrees. The output type is complex unless the + * input is double or dpcomplex, in which case the output is dpcomplex. + * + * Returns: 0 on success, -1 on error + */ +int +vips_complex2( VipsImage *left, VipsImage *right, VipsImage **out, + VipsOperationComplex2 cmplx, ... ) +{ + va_list ap; + int result; + + va_start( ap, cmplx ); + result = vips_complex2v( left, right, out, cmplx, ap ); + va_end( ap ); + + return( result ); +} + +/** + * vips_cross_phase: + * @left: input #VipsImage + * @right: input #VipsImage + * @out: output #VipsImage + * @...: %NULL-terminated list of optional named arguments + * + * Perform #VIPS_OPERATION_COMPLEX2_CROSS_PHASE on an image. + * See vips_complex2(). + * + * Returns: 0 on success, -1 on error + */ +int +vips_cross_phase( VipsImage *left, VipsImage *right, VipsImage **out, ... ) +{ + va_list ap; + int result; + + va_start( ap, out ); + result = vips_complex2v( left, right, out, + VIPS_OPERATION_COMPLEX2_CROSS_PHASE, ap ); + va_end( ap ); + + return( result ); +} + typedef struct _VipsComplexget { VipsUnary parent_instance; diff --git a/libvips/arithmetic/im_cross_phase.c b/libvips/arithmetic/im_cross_phase.c deleted file mode 100644 index c16fba1b..00000000 --- a/libvips/arithmetic/im_cross_phase.c +++ /dev/null @@ -1,150 +0,0 @@ -/* im_cross_phase.c - * - * Copyright: 2008, Nottingham Trent University - * - * Author: Tom Vajzovic - * Written on: 2008-01-09 - * - * 2008-02-04 tcv: - * - exp( i.th ) == cos(th)+i.sin(th) NOT sin(th)+i.cos(th) - * - add quadratic version (ifdef'd out ATM - still using trigonometric one) - * 2/9/09 - * - gtk-doc comment - */ - -/* - - 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 - -/* There doesn't seem to be much difference in speed between these two methods (on an Athlon64), - * so I use the modulus argument version, since atan2() is in c89 but hypot() is c99. - * - * If you think that it might be faster on your platform, uncomment the following: - */ -#define USE_MODARG_DIV - -#ifdef USE_MODARG_DIV - -#define COMPLEX_PHASE_FN( TYPE, ABS ) \ -static void \ -complex_phase_ ## TYPE ( void *in1, void *in2, void *out, int n, void *im, void *unrequired ){ \ - \ - TYPE *X= (TYPE*) in1; \ - TYPE *Y= (TYPE*) in2; \ - TYPE *Z= (TYPE*) out; \ - TYPE *Z_stop= Z + 2 * n * ((IMAGE*)im)-> Bands; \ - \ - for( ; Z < Z_stop; X+= 2, Y+= 2 ){ \ - double arg= atan2( X[1], X[0] ) - atan2( Y[1], Y[0] ); \ - *Z++= cos( arg ); \ - *Z++= sin( arg ); \ - } \ -} - -#else /* USE_MODARG_DIV */ - -#define COMPLEX_PHASE_FN( TYPE, ABS ) \ -static void \ -complex_phase_ ## TYPE ( void *in1, void *in2, void *out, int n, void *im, void *unrequired ){ \ - \ - TYPE *X= (TYPE*) in1; \ - TYPE *Y= (TYPE*) in2; \ - TYPE *Z= (TYPE*) out; \ - TYPE *Z_stop= Z + 2 * n * ((IMAGE*)im)-> Bands; \ - \ - for( ; Z < Z_stop; X+= 2, Y+= 2 ) \ - \ - if( ABS( Y[0] ) > ABS( Y[1] )){ \ - double a= Y[1] / Y[0]; \ - double b= Y[0] + Y[1] * a; \ - double re= ( X[0] + X[1] * a ) / b; \ - double im= ( X[1] - X[0] * a ) / b; \ - double mod= im__hypot( re, im ); \ - *Z++= re / mod; \ - *Z++= im / mod; \ - } \ - else { \ - double a= Y[0] / Y[1]; \ - double b= Y[1] + Y[0] * a; \ - double re= ( X[0] * a + X[1] ) / b; \ - double im= ( X[1] * a - X[0] ) / b; \ - double mod= im__hypot( re, im ); \ - *Z++= re / mod; \ - *Z++= im / mod; \ - } \ -} - -#endif /* USE_MODARG_DIV */ - -COMPLEX_PHASE_FN( float, fabsf ) -COMPLEX_PHASE_FN( double, fabs ) - -/** - * im_cross_phase: - * @a: input #IMAGE 1 - * @b: input #IMAGE 2 - * @out: output #IMAGE - * - * Find the phase of the cross power spectrum of two complex images, - * expressed as a complex image where the modulus of each pixel is - * one. - * - * I.E. find (a.b*)/|a.b*| where - * . represents complex multiplication - * * represents the complex conjugate - * || represents the complex modulus - * - * See also: im_multiply(), im_sign(). - * - * Returns: 0 on success, -1 on error - */ -int im_cross_phase( IMAGE *a, IMAGE *b, IMAGE *out ){ -#define FUNCTION_NAME "im_phase" - - if( im_pincheck( a ) || im_pincheck( b ) || im_poutcheck( out )) - return -1; - - if( im_check_size_same( FUNCTION_NAME, a, b ) || - im_check_bands_same( FUNCTION_NAME, a, b ) || - im_check_format_same( FUNCTION_NAME, a, b ) || - im_check_uncoded( FUNCTION_NAME, a ) || - im_check_uncoded( FUNCTION_NAME, b ) || - im_check_complex( FUNCTION_NAME, a ) || - im_check_complex( FUNCTION_NAME, b ) ) - return -1; - - return im_cp_descv( out, a, b, NULL ) || im_wraptwo( a, b, out, - IM_BANDFMT_COMPLEX == a-> BandFmt ? complex_phase_float : complex_phase_double, a, NULL ); -} diff --git a/libvips/arithmetic/math2.c b/libvips/arithmetic/math2.c index d926fc3c..f44b96ef 100644 --- a/libvips/arithmetic/math2.c +++ b/libvips/arithmetic/math2.c @@ -196,7 +196,7 @@ vips_math2_class_init( VipsMath2Class *class ) gobject_class->get_property = vips_object_get_property; object_class->nickname = "math2"; - object_class->description = _( "pow( left, right)" ); + object_class->description = _( "binary math operations" ); object_class->build = vips_math2_build; vips_arithmetic_set_format_table( aclass, vips_bandfmt_math2 ); diff --git a/libvips/include/vips/arithmetic.h b/libvips/include/vips/arithmetic.h index 4ec79cb4..5dbad4f3 100644 --- a/libvips/include/vips/arithmetic.h +++ b/libvips/include/vips/arithmetic.h @@ -147,6 +147,17 @@ typedef enum { VIPS_OPERATION_COMPLEX_LAST } VipsOperationComplex; +/** + * VipsOperationComplex2: + * @VIPS_OPERATION_COMPLEX2_CROSS_PHASE: convert to polar coordinates + * + * See also: vips_complex2(). + */ +typedef enum { + VIPS_OPERATION_COMPLEX2_CROSS_PHASE, + VIPS_OPERATION_COMPLEX2_LAST +} VipsOperationComplex2; + /** * VipsOperationComplexget: * @VIPS_OPERATION_COMPLEXGET_REAL: get real component @@ -230,6 +241,12 @@ int vips_rect( VipsImage *in, VipsImage **out, ... ) int vips_conj( VipsImage *in, VipsImage **out, ... ) __attribute__((sentinel)); +int vips_complex2( VipsImage *left, VipsImage *right, VipsImage **out, + VipsOperationComplex2 cmplx, ... ) + __attribute__((sentinel)); +int vips_cross_phase( VipsImage *left, VipsImage *right, VipsImage **out, ... ) + __attribute__((sentinel)); + int vips_complexget( VipsImage *in, VipsImage **out, VipsOperationComplexget get, ... ) __attribute__((sentinel)); diff --git a/libvips/include/vips/enumtypes.h b/libvips/include/vips/enumtypes.h index 54c9c742..5eb5a758 100644 --- a/libvips/include/vips/enumtypes.h +++ b/libvips/include/vips/enumtypes.h @@ -34,6 +34,8 @@ GType vips_operation_boolean_get_type (void) G_GNUC_CONST; #define VIPS_TYPE_OPERATION_BOOLEAN (vips_operation_boolean_get_type()) GType vips_operation_complex_get_type (void) G_GNUC_CONST; #define VIPS_TYPE_OPERATION_COMPLEX (vips_operation_complex_get_type()) +GType vips_operation_complex2_get_type (void) G_GNUC_CONST; +#define VIPS_TYPE_OPERATION_COMPLEX2 (vips_operation_complex2_get_type()) GType vips_operation_complexget_get_type (void) G_GNUC_CONST; #define VIPS_TYPE_OPERATION_COMPLEXGET (vips_operation_complexget_get_type()) /* enumerations from "../../../libvips/include/vips/conversion.h" */ diff --git a/libvips/include/vips/inlines.h b/libvips/include/vips/inlines.h index 9d1a60e6..d86da108 100644 --- a/libvips/include/vips/inlines.h +++ b/libvips/include/vips/inlines.h @@ -27,44 +27,46 @@ */ -#ifndef IM_INLINE_H -#define IM_INLINE_H +#ifndef VIPS_INLINE_H +#define VIPS_INLINE_H -/* glib promises to define inline in a portable way */ +/* glib promises to define inline in a portable way + */ #include - #ifdef __cplusplus extern "C" { #endif /*__cplusplus*/ - #ifdef HAVE_HYPOT -#define im__hypot hypot +#define vips__hypot hypot #else /* HAVE_HYPOT */ -static inline double im__hypot( double a, double b ){ - double ta= fabs( a ); - double tb= fabs( b ); +static inline double +vips__hypot( double a, double b ) +{ + double ta = fabs( a ); + double tb = fabs( b ); - if( ta > tb ){ - tb= b / a; - return ta * sqrt( 1.0 + tb * tb ); - } - else { - ta= a / b; - return tb * sqrt( 1.0 + ta * ta ); - } + if( ta > tb ) { + tb = b / a; + + return( ta * sqrt( 1.0 + tb * tb ) ); + } + else { + ta = a / b; + + return( tb * sqrt( 1.0 + ta * ta ) ); + } } #endif /* HAVE_HYPOT */ - #ifdef __cplusplus } #endif /*__cplusplus*/ -#endif /*IM_INLINE_H*/ +#endif /*VIPS_INLINE_H*/ diff --git a/libvips/iofuncs/enumtypes.c b/libvips/iofuncs/enumtypes.c index b8ee6f33..1b5a7e38 100644 --- a/libvips/iofuncs/enumtypes.c +++ b/libvips/iofuncs/enumtypes.c @@ -368,6 +368,23 @@ vips_operation_complex_get_type( void ) return( etype ); } GType +vips_operation_complex2_get_type( void ) +{ + static GType etype = 0; + + if( etype == 0 ) { + static const GEnumValue values[] = { + {VIPS_OPERATION_COMPLEX2_CROSS_PHASE, "VIPS_OPERATION_COMPLEX2_CROSS_PHASE", "cross-phase"}, + {VIPS_OPERATION_COMPLEX2_LAST, "VIPS_OPERATION_COMPLEX2_LAST", "last"}, + {0, NULL, NULL} + }; + + etype = g_enum_register_static( "VipsOperationComplex2", values ); + } + + return( etype ); +} +GType vips_operation_complexget_get_type( void ) { static GType etype = 0;