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