partial im_maxpos

This commit is contained in:
John Cupitt 2009-09-08 16:32:55 +00:00
parent 5901a7ca87
commit 155cf95543
6 changed files with 196 additions and 359 deletions

View File

@ -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

View File

@ -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 \

View File

@ -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: <link linkend="VIPS-iofuncs">iofuncs</link>
* @stability: Stable

View File

@ -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 <config.h>
#endif /*HAVE_CONFIG_H*/
#include <vips/intl.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <vips/vips.h>
#include <vips/internal.h>
#ifdef WITH_DMALLOC
#include <dmalloc.h>
#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 = &reg->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 );
}

View File

@ -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 <config.h>
#endif /*HAVE_CONFIG_H*/
#include <vips/intl.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <vips/vips.h>
#include <vips/internal.h>
#ifdef WITH_DMALLOC
#include <dmalloc.h>
#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 = &reg->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; y<in->Ysize; y++ ) \
for ( x=0; x<os; x++ ) {\
if( (double) *p > 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; y<in->Ysize; y++ ) \
for ( x=0; x<os; x++ ) {\
re = (double)*p++; im = (double)*p++; \
mod = re * re + im * im; \
if( mod > 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 ) );
}

View File

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