add binary complex ops

cross_phase is the only one atm

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

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()
 

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 \

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(); 
 }

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;
 

libvips/arithmetic/im_cross_phase.c

@@ -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 <config.h>
-#endif /*HAVE_CONFIG_H*/
-#include <vips/intl.h>
-
-#include <stdlib.h>
-#include <math.h>
-
-#include <vips/vips.h>
-
-/* 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 );
-}

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 );

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));

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

libvips/include/vips/inlines.h

@@ -27,44 +27,46 @@
 
  */
 
-#ifndef IM_INLINE_H
+#ifndef VIPS_INLINE_H
-#define IM_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 <glib.h>
 
-
 #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 ){
+static inline double 
-  double ta= fabs( a );
+vips__hypot( double a, double b ) 
-  double tb= fabs( b );
+{
+	double ta = fabs( a );
+	double tb = fabs( b );
 
-  if( ta > tb ){
+	if( ta > tb ) {
-    tb= b / a;
+		tb = b / a;
-    return ta * sqrt( 1.0 + tb * tb );
+
-  }
+		return( ta * sqrt( 1.0 + tb * tb ) );
-  else {
+	}
-    ta= a / b;
+	else {
-    return tb * sqrt( 1.0 + ta * ta );
+		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*/

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;