diff --git a/ChangeLog b/ChangeLog index 124eea22..4595bb7b 100644 --- a/ChangeLog +++ b/ChangeLog @@ -29,6 +29,8 @@ - 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 +- im_maxpos() is partial and works for complex, im_max() is now a convenience + function 25/3/09 started 7.18.0 - revised version numbers diff --git a/libvips/arithmetic/Makefile.am b/libvips/arithmetic/Makefile.am index b97abfb7..e0957222 100644 --- a/libvips/arithmetic/Makefile.am +++ b/libvips/arithmetic/Makefile.am @@ -19,7 +19,6 @@ libarithmetic_la_SOURCES = \ im_lintra.c \ im_log10tra.c \ im_logtra.c \ - im_max.c \ im_maxpos.c \ im_maxpos_avg.c \ im_maxpos_vec.c \ diff --git a/libvips/arithmetic/arith_dispatch.c b/libvips/arithmetic/arith_dispatch.c index 140f245a..1b4dcacf 100644 --- a/libvips/arithmetic/arith_dispatch.c +++ b/libvips/arithmetic/arith_dispatch.c @@ -44,7 +44,8 @@ /** * SECTION: arithmetic - * @short_description: operations which perform pixel arithmetic, trig, log + * @short_description: operations which perform pixel arithmetic, trig, log, + * stats * * @see_also: iofuncs * @stability: Stable diff --git a/libvips/arithmetic/im_max.c b/libvips/arithmetic/im_max.c deleted file mode 100644 index 8009fb7e..00000000 --- a/libvips/arithmetic/im_max.c +++ /dev/null @@ -1,277 +0,0 @@ -/* im_max.c - * - * Copyright: 1990, J. Cupitt - * - * Author: J. Cupitt - * Written on: 02/05/1990 - * Modified on : 18/03/1991, N. Dessipris - * 7/7/93 JC - * - complex case fixed - * - im_incheck() call added - * 20/6/95 JC - * - now returns double - * - modernised a little - * - now returns max square amplitude rather than amplitude for complex - * 9/5/02 JC - * - partialed - * 3/4/02 JC - * - random wrong result for >1 thread :-( (thanks Joe) - * 15/10/07 - * - oh, heh, seq->inf was not being set correctly, not that it mattered - * 4/9/09 - * - rewrite with im__value(), much simpler and fixes a race condition - * - gtkdoc 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 - - */ - -/* -#define DEBUG - */ - -#ifdef HAVE_CONFIG_H -#include -#endif /*HAVE_CONFIG_H*/ -#include - -#include -#include -#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 *global_max = (double *) b; - double *max; - - if( !(max = IM_NEW( NULL, double )) ) - return( NULL ); - *max = *global_max; - - return( (void *) max ); -} - -/* Merge the sequence value back into the per-call state. - */ -static int -max_stop( void *seq, void *a, void *b ) -{ - double *max = (double *) seq; - double *global_max = (double *) b; - - /* Merge. - */ - *global_max = IM_MAX( *global_max, *max ); - - 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 -max_scan( void *in, int n, void *seq, void *a, void *b ) -{ - const IMAGE *im = (IMAGE *) a; - const int sz = n * im->Bands; - - double *max = (double *) seq; - - int x; - double m; - - m = *max; - - 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: CLOOP( float ); break; - case IM_BANDFMT_DPCOMPLEX: CLOOP( double ); break; - - default: - g_assert( 0 ); - } - - *max = m; - - return( 0 ); -} - -/** - * im_max: - * @in: input #IMAGE - * @out: output double - * - * Finds the the maximum value of image #in and returns it at the - * location pointed by out. If input is complex, the max modulus - * is returned. im_max() finds the maximum of all bands: if you - * want to find the maximum of each band separately, use im_stats(). - * - * See also: im_maxpos(), im_min(), im_stats(). - * - * Returns: 0 on success, -1 on error - */ -int -im_max( IMAGE *in, double *out ) -{ - double global_max; - - if( im_pincheck( in ) || - im_check_uncoded( "im_max", in ) ) - return( -1 ); - - if( im__value( in, &global_max ) ) - return( -1 ); - /* We use square mod for scanning, for speed. - */ - if( im_iscomplex( in ) ) - global_max *= global_max; - - if( im__wrapscan( in, max_start, max_scan, max_stop, - in, &global_max ) ) - return( -1 ); - - /* Back to modulus. - */ - if( im_iscomplex( in ) ) - global_max = sqrt( global_max ); - - *out = global_max; - - return( 0 ); -} diff --git a/libvips/arithmetic/im_maxpos.c b/libvips/arithmetic/im_maxpos.c index dec838e8..432741a7 100644 --- a/libvips/arithmetic/im_maxpos.c +++ b/libvips/arithmetic/im_maxpos.c @@ -12,6 +12,9 @@ * - now returns double for value, like im_max() * 4/9/09 * - gtkdoc comment + * 8/9/09 + * - rewrite based on im_max() to get partial + * - move im_max() in here as a convenience function */ /* @@ -40,32 +43,143 @@ */ +/* +#define DEBUG + */ + #ifdef HAVE_CONFIG_H #include #endif /*HAVE_CONFIG_H*/ #include #include +#include +#include #include +#include #ifdef WITH_DMALLOC #include #endif /*WITH_DMALLOC*/ -/* Useful: Call a macro with the name, type pairs for all VIPS functions. +/* A position and maximum. */ -#define im_for_all_types() \ - 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: loopcmplx(float); break; \ - case IM_BANDFMT_DPCOMPLEX: loopcmplx(double); break; +typedef struct _Maxpos { + int xpos; + int ypos; + double max; +} Maxpos; + +/* New sequence value. + */ +static void * +maxpos_start( IMAGE *in, void *a, void *b ) +{ + Maxpos *global_maxpos = (Maxpos *) b; + Maxpos *maxpos; + + if( !(maxpos = IM_NEW( NULL, Maxpos )) ) + return( NULL ); + *maxpos = *global_maxpos; + + return( (void *) maxpos ); +} + +/* Merge the sequence value back into the per-call state. + */ +static int +maxpos_stop( void *seq, void *a, void *b ) +{ + Maxpos *global_maxpos = (Maxpos *) b; + Maxpos *maxpos = (Maxpos *) seq; + + /* Merge. + */ + if( maxpos->max > global_maxpos->max ) + *global_maxpos = *maxpos; + + im_free( seq ); + + return( 0 ); +} + +#define LOOP( TYPE ) { \ + TYPE *p = (TYPE *) in; \ + \ + for( x = 0; x < sz; x++ ) { \ + double v = p[x]; \ + \ + if( v > max ) { \ + max = v; \ + xpos = r->left + x / reg->im->Bands; \ + ypos = r->top + y; \ + } \ + } \ +} + +#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 > max ) { \ + max = mod; \ + xpos = r->left + x / reg->im->Bands; \ + ypos = r->top + y; \ + } \ + } \ +} + +/* Loop over region, adding to seq. + */ +static int +maxpos_scan( REGION *reg, void *seq, void *a, void *b ) +{ + const Rect *r = ®->valid; + const int sz = IM_REGION_N_ELEMENTS( reg ); + Maxpos *maxpos = (Maxpos *) seq; + + int x, y; + double max; + int xpos, ypos; + + xpos = maxpos->xpos; + ypos = maxpos->ypos; + max = maxpos->max; + + for( y = 0; y < r->height; y++ ) { + PEL *in = (PEL *) IM_REGION_ADDR( reg, r->left, r->top + y ); + + switch( reg->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: CLOOP( float ); break; + case IM_BANDFMT_DPCOMPLEX: CLOOP( double ); break; + + default: + g_assert( 0 ); + } + } + + maxpos->xpos = xpos; + maxpos->ypos = ypos; + maxpos->max = max; + + return( 0 ); +} /** * im_maxpos: @@ -75,10 +189,8 @@ * @out: returned pixel value at that position * * Function to find the maximum of an image. Works for any - * image type. Returns a double and the location of max. - * - * This is not a PIO operation! It may use a lot of memory and take a while. - * Needs a rewrite. + * image type. Returns a double and the location of max. For complex images, + * finds the pixel with the highest modulus. * * See also: im_minpos(), im_min(), im_stats(), im_maxpos_avg(). * @@ -87,69 +199,59 @@ int im_maxpos( IMAGE *in, int *xpos, int *ypos, double *out ) { - double m; - int xp=0, yp=0; - int os; + Maxpos *global_maxpos; -/* Check our args. */ - if( im_incheck( in ) ) + if( im_pincheck( in ) || + im_check_uncoded( "im_maxpos", in ) ) return( -1 ); - if( in->Coding != IM_CODING_NONE ) { - im_error( "im_maxpos", "%s", _( "not uncoded" ) ); + + if( !(global_maxpos = IM_NEW( in, Maxpos )) ) return( -1 ); - } + if( im__value( in, &global_maxpos->max ) ) + return( -1 ); + global_maxpos->xpos = 0; + global_maxpos->ypos = 0; -/* What type? First define the loop we want to perform for all types. */ -#define loop(TYPE) \ - { TYPE *p = (TYPE *) in->data; \ - int x, y; \ - m = (double) *p; \ - \ - for ( y=0; yYsize; y++ ) \ - for ( x=0; x m ) {\ - m = (double) *p; \ - xp = x; yp = y; \ - }\ - p++ ;\ - }\ - } - -#define loopcmplx(TYPE) \ - { TYPE *p = (TYPE *) in->data; \ - double re=(double)*p;\ - double im=(double)*(p+1);\ - double mod = re * re + im * im;\ - int x, y; \ - m = mod; \ - \ - for ( y=0; yYsize; y++ ) \ - for ( x=0; x m ) {\ - m = mod; \ - xp = x; yp = y; \ - }\ - }\ - } - -/* Now generate code for all types. */ - os = in->Xsize * in->Bands; - switch( in->BandFmt ) { - im_for_all_types(); - default: { - g_assert( 0 ); - return( -1 ); - } - } - - /* Return maxima and position of maxima. Nasty: we divide the xpos by - * the number of bands to get the position in pixels. + /* We use square mod for scanning, for speed. */ - *out = m; - *xpos = xp / in->Bands; - *ypos = yp; + if( im_iscomplex( in ) ) + global_maxpos->max *= global_maxpos->max; + + if( im_iterate( in, maxpos_start, maxpos_scan, maxpos_stop, + in, global_maxpos ) ) + return( -1 ); + + /* Back to modulus. + */ + if( im_iscomplex( in ) ) + global_maxpos->max = sqrt( global_maxpos->max ); + + if( xpos ) + *xpos = global_maxpos->xpos; + if( ypos ) + *ypos = global_maxpos->ypos; + if( out ) + *out = global_maxpos->max; return( 0 ); } + +/** + * im_max: + * @in: input #IMAGE + * @out: output double + * + * Finds the the maximum value of image #in and returns it at the + * location pointed by out. If input is complex, the max modulus + * is returned. im_max() finds the maximum of all bands: if you + * want to find the maximum of each band separately, use im_stats(). + * + * See also: im_maxpos(), im_min(), im_stats(). + * + * Returns: 0 on success, -1 on error + */ +int +im_max( IMAGE *in, double *out ) +{ + return( im_maxpos( in, NULL, NULL, out ) ); +} diff --git a/libvips/arithmetic/im_stats.c b/libvips/arithmetic/im_stats.c index ab7172ff..c7aa6254 100644 --- a/libvips/arithmetic/im_stats.c +++ b/libvips/arithmetic/im_stats.c @@ -1,4 +1,5 @@ -/* @(#) im_stats: find general image statistics for all bands separately +/* im_stats.c + * (C) Kirk Martinez 1993 23/4/93 J.Cupitt - adapted to partial images @@ -22,6 +23,7 @@ * some architectures * 7/9/09 * - rework based on new im__wrapscan() / im_max() ideas for a proper fix + * - gtkdoc comment */ /* @@ -180,11 +182,19 @@ stats_scan( void *in, int n, void *seq, void *a, void *b ) return( 0 ); } -/* Find the statistics of an image. Take any non-complex format. Write the - * stats to a DOUBLEMASK of size 6 by (in->Bands+1). We hold a row for each - * band, plus one row for all bands. Row n has 6 elements, which are, in - * order, (minimum, maximum, sum, sum^2, mean, deviation) for band n. Row 0 has - * the figures for all bands together. +/** + * im_stats: + * @in: image to scan + * + * Find many image statistics in a single pass through the data. Returns a + * #DOUBLEMASK of 6 columns by n + 1 (where n is number of bands in image @in) + * rows. Columns are statistics, and are, in order: minimum, maximum, sum, + * sum of squares, mean, standard deviation. Row 0 has statistics for all + * bands together, row 1 has stats for band 1, and so on. + * + * See also: im_maxpos(), im_min(), im_deviate(). + * + * Returns: 0 on success, -1 on error */ DOUBLEMASK * im_stats( IMAGE *im )