diff --git a/ChangeLog b/ChangeLog index f6c3e50e..124eea22 100644 --- a/ChangeLog +++ b/ChangeLog @@ -26,8 +26,9 @@ - im_remainder() has a float result for float types - im_measure() allows sel == NULL, meaning all patches - added "deprecated" package -- faster, simpler, better im_max() +- faster, simpler, better im_max(), im_min, im_avg(), im_deviate() - im_max() returns true modulus, not square of modulus, for complex images +- im_avg() works for complex, returning average modulus 25/3/09 started 7.18.0 - revised version numbers diff --git a/TODO b/TODO index f5a27593..382a62cf 100644 --- a/TODO +++ b/TODO @@ -1,8 +1,4 @@ -- new im_max() needs testing - - im_avg() needs to work for complex: return modulus - - use im_max() model to rewrite im_min(), im_stats() +- im_stats() needs debugging - im_minpos(), im_maxpos() need partialing diff --git a/libvips/arithmetic/im_avg.c b/libvips/arithmetic/im_avg.c index 1b735645..b17a7305 100644 --- a/libvips/arithmetic/im_avg.c +++ b/libvips/arithmetic/im_avg.c @@ -23,6 +23,9 @@ * - add liboil support * 18/8/09 * - gtkdoc, minor reformatting + * 7/9/09 + * - rewrite for im__wrapiter() + * - add complex case (needed for im_max()) */ /* @@ -59,7 +62,6 @@ #include #include #include -#include #include #include @@ -72,30 +74,32 @@ #include #endif /*WITH_DMALLOC*/ -/* Start function: allocate space for a double in which we can accululate the +/* Start function: allocate space for a double in which we can accumulate the * sum. */ static void * -start_fn( IMAGE *out, void *a, void *b ) +avg_start( IMAGE *out, void *a, void *b ) { - double *tmp; + double *sum; - if( !(tmp = IM_ARRAY( out, 1, double )) ) + if( !(sum = IM_NEW( NULL, double )) ) return( NULL ); - *tmp = 0.0; + *sum = 0.0; - return( (void *) tmp ); + return( (void *) sum ); } /* Stop function. Add this little sum to the main sum. */ static int -stop_fn( void *seq, void *a, void *b ) +avg_stop( void *seq, void *a, void *b ) { - double *tmp = (double *) seq; - double *sum = (double *) a; + double *sum = (double *) seq; + double *global_sum = (double *) b; - *sum += *tmp; + *global_sum += *sum; + + im_free( seq ); return( 0 ); } @@ -103,64 +107,76 @@ stop_fn( void *seq, void *a, void *b ) /* Sum pels in this section. */ #define LOOP( TYPE ) { \ - TYPE *p; \ + TYPE *p = (TYPE *) in; \ \ - for( y = to; y < bo; y++ ) { \ - p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \ - \ - for( x = 0; x < sz; x++ ) \ - sum += *p++; \ - } \ + for( x = 0; x < sz; x++ ) \ + m += p[x]; \ } +#define CLOOP( TYPE ) { \ + TYPE *p = (TYPE *) in; \ + \ + for( x = 0; x < sz; x++ ) { \ + double mod, re, im; \ + \ + re = p[0]; \ + im = p[1]; \ + p += 2; \ + mod = re * re + im * im; \ + \ + m += mod; \ + } \ +} + /* Loop over region, accumulating a sum in *tmp. */ static int -scan_fn( REGION *reg, void *seq, void *a, void *b ) -{ - double *tmp = (double *) seq; - Rect *r = ®->valid; - IMAGE *im = reg->im; - int le = r->left; - int to = r->top; - int bo = IM_RECT_BOTTOM(r); - int sz = IM_REGION_N_ELEMENTS( reg ); - double sum = 0.0; - int x, y; +avg_scan( void *in, int n, void *seq, void *a, void *b ) +{ + const IMAGE *im = (IMAGE *) a; + const int sz = n * im->Bands; + + double *sum = (double *) seq; + + int x; + double m; + + m = *sum; /* Now generate code for all types. */ switch( im->BandFmt ) { - case IM_BANDFMT_UCHAR: LOOP( unsigned char ); break; - case IM_BANDFMT_CHAR: LOOP( signed char ); break; - case IM_BANDFMT_USHORT: LOOP( unsigned short ); break; - case IM_BANDFMT_SHORT: LOOP( signed short ); break; - case IM_BANDFMT_UINT: LOOP( unsigned int ); break; - case IM_BANDFMT_INT: LOOP( signed int ); break; - case IM_BANDFMT_FLOAT: LOOP( float ); break; + case IM_BANDFMT_UCHAR: LOOP( unsigned char ); break; + case IM_BANDFMT_CHAR: LOOP( signed char ); break; + case IM_BANDFMT_USHORT: LOOP( unsigned short ); break; + case IM_BANDFMT_SHORT: LOOP( signed short ); break; + case IM_BANDFMT_UINT: LOOP( unsigned int ); break; + case IM_BANDFMT_INT: LOOP( signed int ); break; + case IM_BANDFMT_FLOAT: LOOP( float ); break; case IM_BANDFMT_DOUBLE: #ifdef HAVE_LIBOIL - for( y = to; y < bo; y++ ) { - double *p = (double *) IM_REGION_ADDR( reg, le, y ); - double t; +{ + double *p = (double *) in; + double t; - oil_sum_f64( &t, p, sizeof( double ), sz ); + oil_sum_f64( &t, p, sizeof( double ), sz ); - sum += t; - } + m += t; +} #else /*!HAVE_LIBOIL*/ LOOP( double ); #endif /*HAVE_LIBOIL*/ break; + case IM_BANDFMT_COMPLEX: CLOOP( float ); break; + case IM_BANDFMT_DPCOMPLEX: CLOOP( double ); break; + default: - assert( 0 ); + g_assert( 0 ); } - /* Add to sum for this sequence. - */ - *tmp += sum; + *sum = m; return( 0 ); } @@ -172,9 +188,7 @@ scan_fn( REGION *reg, void *seq, void *a, void *b ) * * This operation finds the average value in an image. It operates on all * bands of the input image: use im_stats() if you need to calculate an - * average for each band. - * - * Non-complex images only. + * average for each band. For complex images, return the average modulus. * * See also: im_stats(), im_bandmean(), im_deviate(), im_rank() * @@ -183,7 +197,7 @@ scan_fn( REGION *reg, void *seq, void *a, void *b ) int im_avg( IMAGE *in, double *out ) { - double sum = 0.0; + double global_sum; gint64 vals, pels; /* Check our args. @@ -195,14 +209,35 @@ im_avg( IMAGE *in, double *out ) /* Loop over input, summing pixels. */ - if( im_iterate( in, start_fn, scan_fn, stop_fn, &sum, NULL ) ) + global_sum = 0.0; + if( im__wrapscan( in, + avg_start, avg_scan, avg_stop, in, &global_sum ) ) return( -1 ); - /* Calculate and return average. + /* Calculate and return average. For complex, we accumulate re*re + + * im*im, so we need to sqrt. */ pels = (gint64) in->Xsize * in->Ysize; vals = pels * in->Bands; - *out = sum / vals; + *out = global_sum / vals; + if( im_iscomplex( in ) ) + *out = sqrt( *out ); + + return( 0 ); +} + +/* Get the value of pixel (0, 0). Use this to init the min/max value for + * im_max()/im_stats()/etc. + */ +int +im__value( IMAGE *im, double *value ) +{ + IMAGE *t; + + if( !(t = im_open_local( im, "im__value", "p" )) || + im_extract_area( im, t, 0, 0, 1, 1 ) || + im_avg( t, value ) ) + return( -1 ); return( 0 ); } diff --git a/libvips/arithmetic/im_deviate.c b/libvips/arithmetic/im_deviate.c index a6a63a79..dcb14da3 100644 --- a/libvips/arithmetic/im_deviate.c +++ b/libvips/arithmetic/im_deviate.c @@ -25,6 +25,8 @@ * 2/9/09 * - gtk-doc comment * - minor reformatting + * 4/9/09 + * - use im__wrapscan() */ /* @@ -77,28 +79,30 @@ * accumulate the sum and the sum of squares. */ static void * -start_fn( IMAGE *out, void *a, void *b ) +deviate_start( IMAGE *out, void *a, void *b ) { - double *tmp; + double *ss2; - if( !(tmp = IM_ARRAY( out, 2, double )) ) + if( !(ss2 = IM_ARRAY( NULL, 2, double )) ) return( NULL ); - tmp[0] = 0.0; - tmp[1] = 0.0; + ss2[0] = 0.0; + ss2[1] = 0.0; - return( (void *) tmp ); + return( (void *) ss2 ); } /* Stop function. Add this little sum to the main sum. */ static int -stop_fn( void *seq, void *a, void *b ) +deviate_stop( void *seq, void *a, void *b ) { - double *tmp = (double *) seq; - double *sum = (double *) a; + double *ss2 = (double *) seq; + double *global_ss2 = (double *) b; - sum[0] += tmp[0]; - sum[1] += tmp[1]; + global_ss2[0] += ss2[0]; + global_ss2[1] += ss2[1]; + + im_free( seq ); return( 0 ); } @@ -106,60 +110,56 @@ stop_fn( void *seq, void *a, void *b ) /* Sum pels in this section. */ #define LOOP( TYPE ) { \ - TYPE *p; \ + TYPE *p = (TYPE *) in; \ \ - for( y = to; y < bo; y++ ) { \ - p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \ + for( x = 0; x < sz; x++ ) { \ + TYPE v = p[x]; \ \ - for( x = 0; x < sz; x++ ) { \ - TYPE v = p[x]; \ - \ - s += v; \ - s2 += (double) v * (double) v; \ - } \ + s += v; \ + s2 += (double) v * (double) v; \ } \ } -/* Loop over region, adding information to the appropriate fields of tmp. +/* Loop over region, accumulating a sum in *tmp. */ static int -scan_fn( REGION *reg, void *seq, void *a, void *b ) -{ - double *tmp = (double *) seq; - Rect *r = ®->valid; - IMAGE *im = reg->im; - int le = r->left; - int to = r->top; - int bo = IM_RECT_BOTTOM(r); - int sz = IM_REGION_N_ELEMENTS( reg ); - double s = 0.0; - double s2 = 0.0; - int x, y; +deviate_scan( void *in, int n, void *seq, void *a, void *b ) +{ + const IMAGE *im = (IMAGE *) a; + const int sz = n * im->Bands; + + double *ss2 = (double *) seq; + + int x; + double s, s2; + + s = ss2[0]; + s2 = ss2[1]; /* Now generate code for all types. */ switch( im->BandFmt ) { - case IM_BANDFMT_UCHAR: LOOP( unsigned char ); break; - case IM_BANDFMT_CHAR: LOOP( signed char ); break; - case IM_BANDFMT_USHORT: LOOP( unsigned short ); break; - case IM_BANDFMT_SHORT: LOOP( signed short ); break; - case IM_BANDFMT_UINT: LOOP( unsigned int ); break; - case IM_BANDFMT_INT: LOOP( signed int ); break; - case IM_BANDFMT_FLOAT: LOOP( float ); break; + case IM_BANDFMT_UCHAR: LOOP( unsigned char ); break; + case IM_BANDFMT_CHAR: LOOP( signed char ); break; + case IM_BANDFMT_USHORT: LOOP( unsigned short ); break; + case IM_BANDFMT_SHORT: LOOP( signed short ); break; + case IM_BANDFMT_UINT: LOOP( unsigned int ); break; + case IM_BANDFMT_INT: LOOP( signed int ); break; + case IM_BANDFMT_FLOAT: LOOP( float ); break; case IM_BANDFMT_DOUBLE: #ifdef HAVE_LIBOIL - for( y = to; y < bo; y++ ) { - double *p = (double *) IM_REGION_ADDR( reg, le, y ); - double t; - double t2; +{ + double *p = (double *) in; + double t; + double t2; - oil_sum_f64( &t, p, sizeof( double ), sz ); - oil_squaresum_f64( &t2, p, sz ); + oil_sum_f64( &t, p, sizeof( double ), sz ); + oil_squaresum_f64( &t2, p, sz ); - s += t; - s2 += t2; - } + s += t; + s2 += t2; +} #else /*!HAVE_LIBOIL*/ LOOP( double ); #endif /*HAVE_LIBOIL*/ @@ -169,10 +169,8 @@ scan_fn( REGION *reg, void *seq, void *a, void *b ) g_assert( 0 ); } - /* Add to sum for this sequence. - */ - tmp[0] += s; - tmp[1] += s2; + ss2[0] += s; + ss2[1] += s2; return( 0 ); } @@ -195,7 +193,7 @@ scan_fn( REGION *reg, void *seq, void *a, void *b ) int im_deviate( IMAGE *in, double *out ) { - double sum[2] = { 0.0, 0.0 }; + double global_ss2[2]; gint64 N; /* Check our args. @@ -207,8 +205,10 @@ im_deviate( IMAGE *in, double *out ) /* Loop over input, summing pixels. */ - if( im_iterate( in, start_fn, scan_fn, stop_fn, &sum, NULL ) ) - return( -1 ); + global_ss2[0] = 0.0; + global_ss2[1] = 0.0; + if( im__wrapscan( in, + deviate_start, deviate_scan, deviate_stop, in, &global_ss2 ) ) /* @@ -221,7 +221,8 @@ im_deviate( IMAGE *in, double *out ) /* Calculate and return deviation. Add a fabs to stop sqrt(<=0). */ N = (gint64) in->Xsize * in->Ysize * in->Bands; - *out = sqrt( fabs( sum[1] - (sum[0] * sum[0] / N) ) / (N - 1) ); + *out = sqrt( fabs( global_ss2[1] - + (global_ss2[0] * global_ss2[0] / N) ) / (N - 1) ); return( 0 ); } diff --git a/libvips/arithmetic/im_max.c b/libvips/arithmetic/im_max.c index 81e8b8bd..8009fb7e 100644 --- a/libvips/arithmetic/im_max.c +++ b/libvips/arithmetic/im_max.c @@ -63,35 +63,106 @@ #include #include +#include #ifdef WITH_DMALLOC #include #endif /*WITH_DMALLOC*/ +typedef struct _Wrapscan { + IMAGE *in; + im_start_fn start; + im__wrapscan_fn scan; + im_stop_fn stop; + void *a; + void *b; +} Wrapscan; + +static void * +wrapscan_start( IMAGE *in, void *a, void *b ) +{ + Wrapscan *wrapscan = (Wrapscan *) a; + + return( wrapscan->start( in, wrapscan->a, wrapscan->b ) ); +} + +static int +wrapscan_stop( void *seq, void *a, void *b ) +{ + Wrapscan *wrapscan = (Wrapscan *) a; + + return( wrapscan->stop( seq, wrapscan->a, wrapscan->b ) ); +} + +static int +wrapscan_scan( REGION *reg, void *seq, void *a, void *b ) +{ + Wrapscan *wrapscan = (Wrapscan *) a; + Rect *r = ®->valid; + int lsk = IM_REGION_LSKIP( reg ); + + int y; + PEL *p; + + p = (PEL *) IM_REGION_ADDR( reg, r->left, r->top ); + for( y = 0; y < r->height; y++ ) { + if( wrapscan->scan( p, r->width, seq, + wrapscan->a, wrapscan->b ) ) + return( -1 ); + p += lsk; + } + + return( 0 ); +} + +/* Like im_iterate(), but the scan function works a line at a time, like + * im_wrap*(). Shared with im_min(), im_deviate() etc. + */ +int +im__wrapscan( IMAGE *in, + im_start_fn start, im__wrapscan_fn scan, im_stop_fn stop, + void *a, void *b ) +{ + Wrapscan wrapscan; + + wrapscan.in = in; + wrapscan.start = start; + wrapscan.scan = scan; + wrapscan.stop = stop; + wrapscan.a = a; + wrapscan.b = b; + + return( im_iterate( in, + wrapscan_start, wrapscan_scan, wrapscan_stop, + &wrapscan, NULL ) ); +} + /* New sequence value. */ static void * max_start( IMAGE *in, void *a, void *b ) { - double *value = (double *) a; - double *seq = IM_NEW( NULL, double ); + double *global_max = (double *) b; + double *max; - *seq = *value; + if( !(max = IM_NEW( NULL, double )) ) + return( NULL ); + *max = *global_max; - return( (void *) seq ); + return( (void *) max ); } /* Merge the sequence value back into the per-call state. */ static int -max_stop( void *vseq, void *a, void *b ) +max_stop( void *seq, void *a, void *b ) { - double *seq = (double *) vseq; - double *value = (double *) a; + double *max = (double *) seq; + double *global_max = (double *) b; /* Merge. */ - *value = IM_MAX( *value, *seq ); + *global_max = IM_MAX( *global_max, *max ); im_free( seq ); @@ -99,53 +170,46 @@ max_stop( void *vseq, void *a, void *b ) } #define LOOP( TYPE ) { \ - for( y = to; y < bo; y++ ) { \ - TYPE *p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \ + TYPE *p = (TYPE *) in; \ + \ + for( x = 0; x < sz; x++ ) { \ + double v = p[x]; \ \ - for( x = 0; x < nel; x++ ) { \ - double v = p[x]; \ - \ - if( v > m ) \ - m = v; \ - } \ + if( v > m ) \ + m = v; \ } \ } #define CLOOP( TYPE ) { \ - for( y = to; y < bo; y++ ) { \ - TYPE *p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \ + TYPE *p = (TYPE *) in; \ + \ + for( x = 0; x < sz; x++ ) { \ + double mod, re, im; \ \ - for( x = 0; x < nel; x++ ) { \ - double mod, re, im; \ - \ - re = p[0]; \ - im = p[1]; \ - p += 2; \ - mod = re * re + im * im; \ - \ - if( mod > m ) \ - m = mod; \ - } \ + re = p[0]; \ + im = p[1]; \ + p += 2; \ + mod = re * re + im * im; \ + \ + if( mod > m ) \ + m = mod; \ } \ } /* Loop over region, adding to seq. */ static int -max_scan( REGION *reg, void *vseq, void *a, void *b ) +max_scan( void *in, int n, void *seq, void *a, void *b ) { - double *seq = (double *) vseq; - Rect *r = ®->valid; - IMAGE *im = reg->im; - int le = r->left; - int to = r->top; - int bo = IM_RECT_BOTTOM(r); - int nel = IM_REGION_N_ELEMENTS( reg ); + const IMAGE *im = (IMAGE *) a; + const int sz = n * im->Bands; - int x, y; + double *max = (double *) seq; + + int x; double m; - m = *seq; + m = *max; switch( im->BandFmt ) { case IM_BANDFMT_UCHAR: LOOP( unsigned char ); break; @@ -163,30 +227,7 @@ max_scan( REGION *reg, void *vseq, void *a, void *b ) g_assert( 0 ); } - *seq = m; - -#ifdef DEBUG - printf( "im_max: left = %d, top = %d, width = %d, height = %d\n", - r->left, r->top, r->width, r->height ); - printf( " (max = %g)\n", *seq ); -#endif /*DEBUG*/ - - return( 0 ); -} - -/* Get the value of pixel (0, 0). Use this to init the min/max value for - * threads. Shared with im_min(), im_stats() etc. This will return mod for - * complex. - */ -int -im__value( IMAGE *im, double *value ) -{ - IMAGE *t; - - if( !(t = im_open_local( im, "im__value", "p" )) || - im_extract_area( im, t, 0, 0, 1, 1 ) || - im_avg( t, value ) ) - return( -1 ); + *max = m; return( 0 ); } @@ -208,28 +249,29 @@ im__value( IMAGE *im, double *value ) int im_max( IMAGE *in, double *out ) { - double value; + double global_max; if( im_pincheck( in ) || im_check_uncoded( "im_max", in ) ) return( -1 ); - if( im__value( in, &value ) ) + if( im__value( in, &global_max ) ) return( -1 ); /* We use square mod for scanning, for speed. */ if( im_iscomplex( in ) ) - value *= value; + global_max *= global_max; - if( im_iterate( in, max_start, max_scan, max_stop, &value, NULL ) ) + if( im__wrapscan( in, max_start, max_scan, max_stop, + in, &global_max ) ) return( -1 ); /* Back to modulus. */ if( im_iscomplex( in ) ) - value = sqrt( value ); + global_max = sqrt( global_max ); - *out = value; + *out = global_max; return( 0 ); } diff --git a/libvips/arithmetic/im_min.c b/libvips/arithmetic/im_min.c index 75cb095f..8f77849c 100644 --- a/libvips/arithmetic/im_min.c +++ b/libvips/arithmetic/im_min.c @@ -16,6 +16,8 @@ * - partialed, based in im_max() * 3/4/02 JC * - random wrong result for >1 thread :-( (thanks Joe) + * 4/9/09 + * - rewrite from im_max() */ /* @@ -56,160 +58,105 @@ #include #include #include -#include #include +#include #ifdef WITH_DMALLOC #include #endif /*WITH_DMALLOC*/ -/* Per-call state. - */ -typedef struct _MinInfo { - /* Parameters. - */ - IMAGE *in; - double *out; - - /* Global min so far. - */ - double value; - int valid; /* zero means value is unset */ -} MinInfo; - -/* Per thread state. - */ -typedef struct _Seq { - MinInfo *inf; - - double value; - int valid; /* zero means value is unset */ -} Seq; - /* New sequence value. */ static void * -start_fn( IMAGE *im, void *a, void *b ) +min_start( IMAGE *in, void *a, void *b ) { - MinInfo *inf = (MinInfo *) a; - Seq *seq = IM_NEW( NULL, Seq ); + double *global_min = (double *) b; + double *min; - seq->inf = inf; - seq->valid = 0; + if( !(min = IM_NEW( NULL, double )) ) + return( NULL ); + *min = *global_min; - return( seq ); + return( (void *) min ); } /* Merge the sequence value back into the per-call state. */ static int -stop_fn( void *vseq, void *a, void *b ) +min_stop( void *seq, void *a, void *b ) { - Seq *seq = (Seq *) vseq; - MinInfo *inf = (MinInfo *) a; + double *min = (double *) seq; + double *global_min = (double *) b; - if( seq->valid ) { - if( !inf->valid ) - /* Just copy. - */ - inf->value = seq->value; - else - /* Merge. - */ - inf->value = IM_MIN( inf->value, seq->value ); - - inf->valid = 1; - } + /* Merge. + */ + *global_min = IM_MIN( *global_min, *min ); im_free( seq ); return( 0 ); } +#define LOOP( TYPE ) { \ + TYPE *p = (TYPE *) in; \ + \ + for( x = 0; x < sz; x++ ) { \ + double v = p[x]; \ + \ + if( v < m ) \ + m = v; \ + } \ +} + +#define CLOOP( TYPE ) { \ + TYPE *p = (TYPE *) in; \ + \ + for( x = 0; x < sz; x++ ) { \ + double mod, re, im; \ + \ + re = p[0]; \ + im = p[1]; \ + p += 2; \ + mod = re * re + im * im; \ + \ + if( mod < m ) \ + m = mod; \ + } \ +} + /* Loop over region, adding to seq. */ static int -scan_fn( REGION *reg, void *vseq, void *a, void *b ) +min_scan( void *in, int n, void *seq, void *a, void *b ) { - Seq *seq = (Seq *) vseq; - Rect *r = ®->valid; - IMAGE *im = reg->im; - int le = r->left; - int to = r->top; - int bo = IM_RECT_BOTTOM(r); - int nel = IM_REGION_N_ELEMENTS( reg ); + const IMAGE *im = (IMAGE *) a; + const int sz = n * im->Bands; - int x, y; + double *min = (double *) seq; + int x; double m; -#ifdef DEBUG - printf( "im_min: left = %d, top = %d, width = %d, height = %d\n", - r->left, r->top, r->width, r->height ); -#endif /*DEBUG*/ - -#define loop(TYPE) { \ - m = *((TYPE *) IM_REGION_ADDR( reg, le, to )); \ - \ - for( y = to; y < bo; y++ ) { \ - TYPE *p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \ - \ - for( x = 0; x < nel; x++ ) { \ - double v = p[x]; \ - \ - if( v < m ) \ - m = v; \ - } \ - } \ -} - -#define complex_loop(TYPE) { \ - TYPE *p = (TYPE *) IM_REGION_ADDR( reg, le, to ); \ - double real = p[0]; \ - double imag = p[1]; \ - \ - m = real * real + imag * imag; \ - \ - for( y = to; y < bo; y++ ) { \ - TYPE *p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \ - \ - for( x = 0; x < nel * 2; x += 2 ) { \ - double mod; \ - \ - real = p[x]; \ - imag = p[x + 1]; \ - mod = real * real + imag * imag; \ - \ - if( mod < m ) \ - m = mod; \ - } \ - } \ -} + m = *min; switch( im->BandFmt ) { - case IM_BANDFMT_UCHAR: loop( unsigned char ); break; - case IM_BANDFMT_CHAR: loop( signed char ); break; - case IM_BANDFMT_USHORT: loop( unsigned short ); break; - case IM_BANDFMT_SHORT: loop( signed short ); break; - case IM_BANDFMT_UINT: loop( unsigned int ); break; - case IM_BANDFMT_INT: loop( signed int ); break; - case IM_BANDFMT_FLOAT: loop( float ); break; - case IM_BANDFMT_DOUBLE: loop( double ); break; - case IM_BANDFMT_COMPLEX: complex_loop( float ); break; - case IM_BANDFMT_DPCOMPLEX: complex_loop( double ); break; + case IM_BANDFMT_UCHAR: LOOP( unsigned char ); break; + case IM_BANDFMT_CHAR: LOOP( signed char ); break; + case IM_BANDFMT_USHORT: LOOP( unsigned short ); break; + case IM_BANDFMT_SHORT: LOOP( signed short ); break; + case IM_BANDFMT_UINT: LOOP( unsigned int ); break; + case IM_BANDFMT_INT: LOOP( signed int ); break; + case IM_BANDFMT_FLOAT: LOOP( float ); break; + case IM_BANDFMT_DOUBLE: LOOP( double ); break; + case IM_BANDFMT_COMPLEX: CLOOP( float ); break; + case IM_BANDFMT_DPCOMPLEX: CLOOP( double ); break; default: - assert( 0 ); + g_assert( 0 ); } - if( seq->valid ) { - seq->value = IM_MIN( seq->value, m ); - } - else { - seq->value = m; - seq->valid = 1; - } + *min = m; return( 0 ); } @@ -219,10 +166,10 @@ scan_fn( REGION *reg, void *vseq, void *a, void *b ) * @in: input #IMAGE * @out: output double * - * Finds the the minimum value of image @in and returns it at the - * location pointed by @out. If input is complex, the min square amplitude - * (re*re+im*im) is returned. im_min() finds the maximum of all bands: if you - * want to find the maximum of each band separately, use im_stats(). + * Finds the the minimum value of image #in and returns it at the + * location pointed by out. If input is complex, the min modulus + * is returned. im_min() finds the minimum of all bands: if you + * want to find the minimum of each band separately, use im_stats(). * * See also: im_minpos(), im_max(), im_stats(). * @@ -230,21 +177,30 @@ scan_fn( REGION *reg, void *vseq, void *a, void *b ) */ int im_min( IMAGE *in, double *out ) -{ - MinInfo inf; - - inf.in = in; - inf.out = out; - inf.valid = 0; +{ + double global_min; if( im_pincheck( in ) || im_check_uncoded( "im_min", in ) ) return( -1 ); - if( im_iterate( in, start_fn, scan_fn, stop_fn, &inf, NULL ) ) + if( im__value( in, &global_min ) ) + return( -1 ); + /* We use square mod for scanning, for speed. + */ + if( im_iscomplex( in ) ) + global_min *= global_min; + + if( im__wrapscan( in, min_start, min_scan, min_stop, + in, &global_min ) ) return( -1 ); - *out = inf.value; + /* Back to modulus. + */ + if( im_iscomplex( in ) ) + global_min = sqrt( global_min ); + + *out = global_min; return( 0 ); } diff --git a/libvips/arithmetic/im_stats.c b/libvips/arithmetic/im_stats.c index 6b86f2c9..b9e8518d 100644 --- a/libvips/arithmetic/im_stats.c +++ b/libvips/arithmetic/im_stats.c @@ -17,12 +17,12 @@ * 13/1/05 * - use 64 bit arithmetic - -1/9/09 - - argh nope min/max was broken again for >1CPU in short pipelines on - some architectures - -*/ + * 1/9/09 + * - argh nope min/max was broken again for >1CPU in short pipelines on + * some architectures + * 7/9/09 + * - rework based on new im__wrapscan() / im_max() ideas for a proper fix + */ /* @@ -62,183 +62,119 @@ #include #include #include -#include #include +#include #ifdef WITH_DMALLOC #include #endif /*WITH_DMALLOC*/ -/* Make and initialise a DOUBLEMASK suitable for grabbing statistics. +/* Track min/max/sum/sum-of-squares for each thread during a scan. */ static void * -make_mask( IMAGE *im, void *a, void *b ) +stats_start( IMAGE *im, void *a, void *b ) { - DOUBLEMASK *out; + double *global_stats = (double *) b; - /* Make temp output. - */ - if( !(out = im_create_dmask( "stats", 6, im->Bands + 1 )) ) + double *stats; + int i; + + if( !(stats = IM_ARRAY( NULL, 4 * im->Bands, double )) ) return( NULL ); - - /* Set offset to magic value: 42 indicates we have not yet initialised - * max and min for this mask. - */ - out->offset = 42; -#ifdef DEBUG - printf( "make_mask: created %p\n", out ); -#endif /*DEBUG*/ + for( i = 0; i < 4 * im->Bands; i++ ) + stats[i] = global_stats[i]; - return( out ); + return( stats ); } -/* Merge a temp DOUBLEMASK into the real DOUBLEMASK. Row 0 is unused, row 1 - * has the stats for band 1. These are: (minimum, maximum, sum, sum^2). If the - * offset of out is 42, then it has not been inited yet and we just copy. +/* Merge a thread's output back into the global stats. */ static int -merge_mask( void *seq, void *a, void *b ) +stats_stop( void *seq, void *a, void *b ) { - DOUBLEMASK *tmp = (DOUBLEMASK *) seq; - DOUBLEMASK *out = (DOUBLEMASK *) a; - double *rowi, *rowo; - int z; + const IMAGE *im = (IMAGE *) a; + double *global_stats = (double *) b; + double *stats = (double *) seq; -#ifdef DEBUG - printf( "merge_mask: tmp = %p, out = %p\n", tmp, out ); - im_print_dmask( tmp ); - im_print_dmask( out ); -#endif /*DEBUG*/ + int i; - /* Merge, or just copy? Also, tmp might be uninited, so allow for that - * too. - */ - if( out->offset == 42 && tmp->offset != 42 ) { -#ifdef DEBUG - printf( "merge_mask: copying\n", tmp ); -#endif /*DEBUG*/ + for( i = 0; i < 4 * im->Bands; i += 4 ) { + global_stats[0] = IM_MIN( global_stats[0], stats[0] ); + global_stats[1] = IM_MAX( global_stats[1], stats[1] ); + global_stats[2] += stats[2]; + global_stats[3] += stats[3]; - /* Copy info from tmp. - */ - for( z = 1; z < tmp->ysize; z++ ) { - rowi = tmp->coeff + z * 6; - rowo = out->coeff + z * 6; - - rowo[0] = rowi[0]; - rowo[1] = rowi[1]; - rowo[2] = rowi[2]; - rowo[3] = rowi[3]; - } - - out->offset = 0; - } - else if( out->offset != 42 && tmp->offset != 42 ) { -#ifdef DEBUG - printf( "merge_mask: merging\n" ); -#endif /*DEBUG*/ - - /* Add info from tmp. - */ - for( z = 1; z < tmp->ysize; z++ ) { - rowi = tmp->coeff + z * 6; - rowo = out->coeff + z * 6; - - rowo[0] = IM_MIN( rowi[0], rowo[0] ); - rowo[1] = IM_MAX( rowi[1], rowo[1] ); - rowo[2] += rowi[2]; - rowo[3] += rowi[3]; - } + global_stats += 4; + stats += 4; } - /* Can now free tmp. - */ - im_free_dmask( tmp ); + im_free( seq ); return( 0 ); } -/* Loop over region, adding information to the appropriate fields of tmp. - * We set max, min, sum, sum of squares. Our caller fills in the rest. +/* We scan lines bands times to avoid repeating band loops. + * Use temp variables of same type for min/max for faster comparisons. + */ +#define LOOP( TYPE ) { \ + for( z = 0; z < im->Bands; z++ ) { \ + TYPE *q = (TYPE *) in + z; \ + double *row = stats + z * 4; \ + TYPE small, big; \ + double sum, sum2; \ + \ + small = row[0]; \ + big = row[1]; \ + sum = row[2]; \ + sum2 = row[3]; \ + \ + for( x = 0; x < n; x++ ) { \ + TYPE value = *q; \ + \ + sum += value;\ + sum2 += (double) value * (double) value;\ + if( value > big ) \ + big = value; \ + else if( value < small ) \ + small = value;\ + \ + q += im->Bands; \ + }\ + \ + row[0] = small; \ + row[1] = big; \ + row[2] = sum; \ + row[3] = sum2; \ + }\ +} + +/* Loop over region, adding to seq. */ static int -scan_fn( REGION *reg, void *seq, void *a, void *b ) +stats_scan( void *in, int n, void *seq, void *a, void *b ) { - DOUBLEMASK *tmp = (DOUBLEMASK *) seq; - Rect *r = ®->valid; - IMAGE *im = reg->im; - int bands = im->Bands; - int le = r->left; - int ri = IM_RECT_RIGHT(r); - int to = r->top; - int bo = IM_RECT_BOTTOM(r); - int x, y, z; + const IMAGE *im = (IMAGE *) a; -/* What type? First define the loop we want to perform for all types. - * We scan lines bands times to avoid repeating band loops. - * Use temp variables of same type for min/max for faster comparisons. - * Use double to sum bands. - */ -#define non_complex_loop(TYPE) \ - { TYPE *p, *q; \ - TYPE value, small, big; \ - double *row; \ - \ - /* Have min and max been initialised? \ - */ \ - if( tmp->offset == 42 ) { \ - /* Init min and max for each band. \ - */ \ - p = (TYPE *) IM_REGION_ADDR( reg, le, to ); \ - for( z = 0; z < bands; z++ ) { \ - row = tmp->coeff + (z + 1) * 6; \ - row[0] = p[z]; \ - row[1] = p[z]; \ - } \ - tmp->offset = 0; \ - } \ - \ - for( y = to; y < bo; y++ ) { \ - p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \ - \ - for( z = 0; z < bands; z++ ) { \ - q = p + z; \ - row = tmp->coeff + (z + 1)*6; \ - small = row[0]; \ - big = row[1]; \ - \ - for( x = le; x < ri; x++ ) { \ - value = *q; \ - q += bands; \ - row[2] += value;\ - row[3] += (double)value*(double)value;\ - if( value > big ) \ - big = value; \ - else if( value < small ) \ - small = value;\ - }\ - \ - row[0] = small; \ - row[1] = big; \ - }\ - }\ - } + double *stats = (double *) seq; + + int x, z; /* Now generate code for all types. */ switch( im->BandFmt ) { - case IM_BANDFMT_UCHAR: non_complex_loop(unsigned char); break; - case IM_BANDFMT_CHAR: non_complex_loop(signed char); break; - case IM_BANDFMT_USHORT: non_complex_loop(unsigned short); break; - case IM_BANDFMT_SHORT: non_complex_loop(signed short); break; - case IM_BANDFMT_UINT: non_complex_loop(unsigned int); break; - case IM_BANDFMT_INT: non_complex_loop(signed int); break; - case IM_BANDFMT_DOUBLE: non_complex_loop(double); break; - case IM_BANDFMT_FLOAT: non_complex_loop(float); break; + case IM_BANDFMT_UCHAR: LOOP( unsigned char ); break; + case IM_BANDFMT_CHAR: LOOP( signed char ); break; + case IM_BANDFMT_USHORT: LOOP( unsigned short ); break; + case IM_BANDFMT_SHORT: LOOP( signed short ); break; + case IM_BANDFMT_UINT: LOOP( unsigned int ); break; + case IM_BANDFMT_INT: LOOP( signed int ); break; + case IM_BANDFMT_DOUBLE: LOOP( double ); break; + case IM_BANDFMT_FLOAT: LOOP( float ); break; default: - assert( 0 ); + g_assert( 0 ); } return( 0 ); @@ -251,58 +187,69 @@ scan_fn( REGION *reg, void *seq, void *a, void *b ) * the figures for all bands together. */ DOUBLEMASK * -im_stats( IMAGE *in ) +im_stats( IMAGE *im ) { DOUBLEMASK *out; - double *row, *base; + double *row; gint64 pels, vals, z; + double *global_stats; + int i, j; + double value; - /* Check our args. - */ - if( im_pincheck( in ) ) + if( im_pincheck( im ) || + im_check_noncomplex( "im_stats", im ) || + im_check_uncoded( "im_stats", im ) ) return( NULL ); - if( im_iscomplex( in ) ) { - im_error( "im_stats", "%s", _( "bad input type" ) ); - return( NULL ); - } - if( in->Coding != IM_CODING_NONE ) { - im_error( "im_stats", "%s", _( "not uncoded" ) ); - return( NULL ); - } - /* Make output. - */ - pels = (gint64) in->Xsize * in->Ysize; - vals = pels * in->Bands; - if( !(out = make_mask( in, NULL, NULL )) ) + if( !(global_stats = IM_ARRAY( im, 4 * im->Bands, double )) ) return( NULL ); + if( im__value( im, &value ) ) + return( NULL ); + for( i = 0; i < 4 * im->Bands; i += 4 ) { + global_stats[i + 0] = value; + global_stats[i + 1] = value; + global_stats[i + 2] = 0.0; + global_stats[i + 3] = 0.0; + } /* Loop over input, calculating min, max, sum, sum^2 for each band * separately. */ - if( im_iterate( in, make_mask, scan_fn, merge_mask, out, NULL ) ) { - im_free_dmask( out ); + if( im__wrapscan( im, stats_start, stats_scan, stats_stop, + im, &global_stats ) ) return( NULL ); - } /* Calculate mean, deviation, plus overall stats. */ - base = out->coeff; - base[0] = base[6]; /* Init global max/min */ - base[1] = base[7]; - for( z = 0; z < in->Bands; z++ ) { - row = base + (z + 1) * 6; - base[0] = IM_MIN( base[0], row[0] ); - base[1] = IM_MAX( base[1], row[1] ); - base[2] += row[2]; - base[3] += row[3]; + if( !(out = im_create_dmask( "stats", 6, im->Bands + 1 )) ) + return( NULL ); + + /* Init global max/min/sum/sum2. + */ + out->coeff[0] = value; + out->coeff[1] = value; + out->coeff[2] = 0.0; + out->coeff[3] = 0.0; + + pels = (gint64) im->Xsize * im->Ysize; + vals = pels * im->Bands; + + for( i = 0; i < im->Bands; i++ ) { + row = out->coeff + (z + 1) * 6; + for( j = 0; j < 4; j++ ) + row[j] = global_stats[i * 4 + j]; + + out->coeff[0] = IM_MIN( out->coeff[0], row[0] ); + out->coeff[1] = IM_MAX( out->coeff[1], row[1] ); + out->coeff[2] += row[2]; + out->coeff[3] += row[3]; row[4] = row[2] / pels; row[5] = sqrt( fabs( row[3] - (row[2] * row[2] / pels) ) / (pels - 1) ); } - base[4] = base[2] / vals; - base[5] = sqrt( fabs( base[3] - (base[2] * base[2] / vals) ) / - (vals - 1) ); + out->coeff[4] = out->coeff[2] / vals; + out->coeff[5] = sqrt( fabs( out->coeff[3] - + (out->coeff[2] * out->coeff[2] / vals) ) / (vals - 1) ); #ifdef DEBUG printf( "im_stats:\n" ); diff --git a/libvips/include/vips/internal.h b/libvips/include/vips/internal.h index 44041e34..4864b239 100644 --- a/libvips/include/vips/internal.h +++ b/libvips/include/vips/internal.h @@ -136,6 +136,10 @@ int im__cast_and_call( IMAGE *in1, IMAGE *in2, IMAGE *out, VipsBandFmt im__format_common( IMAGE *in1, IMAGE *in2 ); int im__math( const char *name, IMAGE *in, IMAGE *out, im_wrapone_fn gen ); int im__value( IMAGE *im, double *value ); +typedef int (*im__wrapscan_fn)( void *p, int n, void *seq, void *a, void *b ); +int im__wrapscan( IMAGE *in, + im_start_fn start, im__wrapscan_fn scan, im_stop_fn stop, + void *a, void *b ); int im__test_kill( IMAGE *im ); void *im__mmap( int fd, int writeable, size_t length, gint64 offset );