From 94eaae68ef7da27faa1e16815a05b06e1ca05f05 Mon Sep 17 00:00:00 2001 From: John Cupitt Date: Tue, 8 Sep 2009 20:32:24 +0000 Subject: [PATCH] new im_maxpos_avg() --- ChangeLog | 1 + TODO | 2 - libvips/arithmetic/im_maxpos_avg.c | 295 +++++++++++++++++------------ 3 files changed, 179 insertions(+), 119 deletions(-) diff --git a/ChangeLog b/ChangeLog index 609e4002..195c7de9 100644 --- a/ChangeLog +++ b/ChangeLog @@ -31,6 +31,7 @@ - im_avg() works for complex, returning average modulus - im_maxpos()/im_minpos() are partial and work for complex - im_max()/im_min() are now convenience functions +- im_maxpos_avg() handles complex and multi-band images 25/3/09 started 7.18.0 - revised version numbers diff --git a/TODO b/TODO index 31ac536b..a39b41cb 100644 --- a/TODO +++ b/TODO @@ -1,6 +1,4 @@ -- maybe rewrite im_maxpos_avg() - - im_maxpos_vec needs gtkdoc diff --git a/libvips/arithmetic/im_maxpos_avg.c b/libvips/arithmetic/im_maxpos_avg.c index cdb799fb..9206544c 100644 --- a/libvips/arithmetic/im_maxpos_avg.c +++ b/libvips/arithmetic/im_maxpos_avg.c @@ -12,12 +12,15 @@ * - renamed avg as sum, a bit clearer * 2/9/09 * - gtkdoc comment + * 8/9/08 + * - rewrite from im_maxpos() + * - now handles many bands, complex, faster */ /* 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 @@ -40,40 +43,163 @@ */ - -/** HEADERS **/ +/* +#define DEBUG + */ #ifdef HAVE_CONFIG_H #include -#endif /* HAVE_CONFIG_H */ +#endif /*HAVE_CONFIG_H*/ #include +#include +#include +#include + #include +#include #ifdef WITH_DMALLOC #include -#endif /*WITH_DMALLOC */ +#endif /*WITH_DMALLOC*/ +/* A position and maximum. + */ +typedef struct _Maxposavg { + int xpos; + int ypos; + double max; + int occurences; +} Maxposavg; -/** LOCAL TYPES **/ +/* New sequence value. + */ +static void * +maxposavg_start( IMAGE *in, void *a, void *b ) +{ + Maxposavg *global_maxposavg = (Maxposavg *) b; + Maxposavg *maxposavg; -typedef struct { - double x_sum; - double y_sum; - double val; - unsigned int occurrences; + if( !(maxposavg = IM_NEW( NULL, Maxposavg )) ) + return( NULL ); + *maxposavg = *global_maxposavg; -} pos_avg_t; + return( (void *) maxposavg ); +} +/* Merge the sequence value back into the per-call state. + */ +static int +maxposavg_stop( void *seq, void *a, void *b ) +{ + Maxposavg *global_maxposavg = (Maxposavg *) b; + Maxposavg *maxposavg = (Maxposavg *) seq; -/** LOCAL FUNCTIONS DECLARATIONS **/ + /* Merge. + */ + if( maxposavg->max > global_maxposavg->max ) + *global_maxposavg = *maxposavg; -static void *maxpos_avg_start( IMAGE *im , void *, void * ); -static int maxpos_avg_scan( REGION *reg, void *seq, void *, void * ); -static int maxpos_avg_stop( void *seq, void *, void * ); + im_free( seq ); + return( 0 ); +} -/** EXPORTED FUNCTION **/ +#define LOOP( TYPE ) { \ + TYPE *p = (TYPE *) in; \ + TYPE m; \ + \ + m = max; \ + \ + for( x = 0; x < sz; x++ ) { \ + TYPE v = p[x]; \ + \ + if( v == m ) { \ + xpos += r->left + x / reg->im->Bands; \ + ypos += r->top + y; \ + occurences += 1; \ + } \ + else if( v > m ) { \ + m = v; \ + xpos = r->left + x / reg->im->Bands; \ + ypos = r->top + y; \ + occurences = 1; \ + } \ + } \ + \ + max = m; \ +} + +#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 ) { \ + xpos += r->left + x / reg->im->Bands; \ + ypos += r->top + y; \ + occurences += 1; \ + } \ + else if( mod > max ) { \ + max = mod; \ + xpos = r->left + x / reg->im->Bands; \ + ypos = r->top + y; \ + occurences = 1; \ + } \ + } \ +} + +/* Loop over region, adding to seq. + */ +static int +maxposavg_scan( REGION *reg, void *seq, void *a, void *b ) +{ + const Rect *r = ®->valid; + const int sz = IM_REGION_N_ELEMENTS( reg ); + Maxposavg *maxposavg = (Maxposavg *) seq; + + int x, y; + double max; + int xpos, ypos, occurences; + + xpos = maxposavg->xpos; + ypos = maxposavg->ypos; + max = maxposavg->max; + occurences = maxposavg->occurences; + + 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 ); + } + } + + maxposavg->xpos = xpos; + maxposavg->ypos = ypos; + maxposavg->max = max; + maxposavg->occurences = occurences; + + return( 0 ); +} /** * im_maxpos_avg: @@ -84,117 +210,52 @@ static int maxpos_avg_stop( void *seq, void *, void * ); * * Function to find the maximum of an image. Returns coords and value at * double precision. In the event of a draw, returns average of all - * drawing coords, and interpolated value at that position. + * drawing coords. * * See also: im_maxpos(), im_min(), im_stats(). * * Returns: 0 on success, -1 on error */ -int im_maxpos_avg( IMAGE *im, double *xpos, double *ypos, double *out ){ -#define FUNCTION_NAME "im_maxpos_avg" +int +im_maxpos_avg( IMAGE *in, double *xpos, double *ypos, double *out ) +{ + Maxposavg *global_maxposavg; - pos_avg_t master= { 0.0, 0.0, 0.0, 0 }; + if( im_pincheck( in ) || + im_check_uncoded( "im_maxpos_avg", in ) ) + return( -1 ); - if( im_pincheck( im ) || - im_check_uncoded( FUNCTION_NAME, im ) || - im_check_noncomplex( FUNCTION_NAME, im ) || - im_check_mono( FUNCTION_NAME, im ) ) - return -1; + if( !(global_maxposavg = IM_NEW( in, Maxposavg )) ) + return( -1 ); + if( im__value( in, &global_maxposavg->max ) ) + return( -1 ); + global_maxposavg->xpos = 0; + global_maxposavg->ypos = 0; + global_maxposavg->occurences = 1; - if( ! xpos || ! ypos || ! out ){ - im_error( FUNCTION_NAME, "%s", _("invalid argument") ); - return -1; - } + /* We use square mod for scanning, for speed. + */ + if( im_iscomplex( in ) ) + global_maxposavg->max *= global_maxposavg->max; - if( im_iterate( im, maxpos_avg_start, maxpos_avg_scan, maxpos_avg_stop, &master, NULL ) ) - return -1; + if( im_iterate( in, maxposavg_start, maxposavg_scan, maxposavg_stop, + in, global_maxposavg ) ) + return( -1 ); - *xpos= master. x_sum / master. occurrences; - *ypos= master. y_sum / master. occurrences; + /* Back to modulus. + */ + if( im_iscomplex( in ) ) + global_maxposavg->max = sqrt( global_maxposavg->max ); - return im_point_bilinear( im, *xpos, *ypos, 0, out ); + if( xpos ) + *xpos = (double) global_maxposavg->xpos / + global_maxposavg->occurences; + if( ypos ) + *ypos = (double) global_maxposavg->ypos / + global_maxposavg->occurences; + if( out ) + *out = global_maxposavg->max; -#undef FUNCTION_NAME -} - -static void *maxpos_avg_start( IMAGE *im, void *a, void *b ){ - pos_avg_t *seq; - - seq= IM_NEW( NULL, pos_avg_t ); - if( ! seq ) - return NULL; - - seq-> x_sum= 0.0; - seq-> y_sum= 0.0; - seq-> val= 0.0; - seq-> occurrences= 0; - - return (void *) seq; -} - -/* should be void (always returns 0) */ -static int maxpos_avg_scan( REGION *reg, void *vseq, void *a, void *b ) { - - pos_avg_t *seq= (pos_avg_t *) vseq; - const int right= reg-> valid. left + reg-> valid. width; - const int bottom= reg-> valid. top + reg-> valid. height; - int x; - int y; - -#define LOOPS(type){ \ - type *read= (type*) IM_REGION_ADDR( reg, reg-> valid. left, reg-> valid. top ) - reg-> valid. left; \ - size_t skip= IM_REGION_LSKIP( reg ) / sizeof( type ); \ - \ - for( y= reg-> valid. top; y < bottom; ++y, read+= skip ) \ - for( x= reg-> valid. left; x < right; ++x ) \ - if( !seq-> occurrences || \ - read[x] > seq-> val ){ \ - seq-> val= read[x]; \ - seq-> x_sum= x; \ - seq-> y_sum= y; \ - seq-> occurrences= 1; \ - } \ - else if( read[x] == seq-> val ){ \ - seq-> x_sum+= x; \ - seq-> y_sum+= y; \ - ++ (seq-> occurrences); \ - } \ -} - - switch( reg-> im-> BandFmt ){ - case IM_BANDFMT_CHAR: LOOPS( gint8 ) break; - case IM_BANDFMT_UCHAR: LOOPS( guint8 ) break; - case IM_BANDFMT_SHORT: LOOPS( gint16 ) break; - case IM_BANDFMT_USHORT: LOOPS( guint16 ) break; - case IM_BANDFMT_INT: LOOPS( gint32 ) break; - case IM_BANDFMT_UINT: LOOPS( guint32 ) break; - case IM_BANDFMT_FLOAT: LOOPS( float ) break; - case IM_BANDFMT_DOUBLE: LOOPS( double ) break; - } - - return 0; -#undef LOOPS -} - -/* should be void (always returns 0) */ -static int maxpos_avg_stop( void *vseq, void *a, void *b ) { - - pos_avg_t *seq = (pos_avg_t *) vseq; - pos_avg_t *master = (pos_avg_t *) a; - - if( !master->occurrences || - seq-> val > master-> val ){ - master-> val= seq-> val; - master-> x_sum= seq-> x_sum; - master-> y_sum= seq-> y_sum; - master-> occurrences= seq-> occurrences; - } - else if( seq-> val == master-> val ){ - master-> x_sum+= seq-> x_sum; - master-> y_sum+= seq-> y_sum; - master-> occurrences+= seq-> occurrences; - } - im_free( seq ); - return 0; + return( 0 ); }