im_minpos() and im_min() rewrite

This commit is contained in:
John Cupitt 2009-09-08 16:57:34 +00:00
parent 155cf95543
commit ffea411943
7 changed files with 265 additions and 286 deletions

View File

@ -29,8 +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
- im_maxpos()/im_minpos() are partial and work for complex
- im_max()/im_min() are now a convenience functions
25/3/09 started 7.18.0
- revised version numbers

4
TODO
View File

@ -1,6 +1,5 @@
- im_stats() needs debugging
- im_minpos(), im_maxpos() need partialing
- maybe rewrite im_maxpos_avg()
- im_maxpos_vec needs gtkdoc
@ -8,7 +7,6 @@
- 1-bit PNG read is broken?
> The bug is that 1bit depth PNG addresses are incorrectly interpreted. At

View File

@ -23,7 +23,6 @@ libarithmetic_la_SOURCES = \
im_maxpos_avg.c \
im_maxpos_vec.c \
im_measure.c \
im_min.c \
im_minpos.c \
im_multiply.c \
im_powtra.c \

View File

@ -26,6 +26,8 @@
* 7/9/09
* - rewrite for im__wrapiter()
* - add complex case (needed for im_max())
* 8/9/09
* - wrapscan stuff moved here
*/
/*
@ -74,6 +76,74 @@
#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 ) );
}
/* Start function: allocate space for a double in which we can accumulate the
* sum.
*/

View File

@ -106,16 +106,21 @@ maxpos_stop( void *seq, void *a, void *b )
#define LOOP( TYPE ) { \
TYPE *p = (TYPE *) in; \
TYPE m; \
\
m = max; \
\
for( x = 0; x < sz; x++ ) { \
double v = p[x]; \
TYPE v = p[x]; \
\
if( v > max ) { \
max = v; \
if( v > m ) { \
m = v; \
xpos = r->left + x / reg->im->Bands; \
ypos = r->top + y; \
} \
} \
\
max = m; \
}
#define CLOOP( TYPE ) { \

View File

@ -1,206 +0,0 @@
/* im_min.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 min square amplitude rather than amplitude for complex
* 9/5/02 JC
* - partialed, based in im_max()
* 3/4/02 JC
* - random wrong result for >1 thread :-( (thanks Joe)
* 4/9/09
* - rewrite from im_max()
*/
/*
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*/
/* New sequence value.
*/
static void *
min_start( IMAGE *in, void *a, void *b )
{
double *global_min = (double *) b;
double *min;
if( !(min = IM_NEW( NULL, double )) )
return( NULL );
*min = *global_min;
return( (void *) min );
}
/* Merge the sequence value back into the per-call state.
*/
static int
min_stop( void *seq, void *a, void *b )
{
double *min = (double *) seq;
double *global_min = (double *) b;
/* 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
min_scan( void *in, int n, void *seq, void *a, void *b )
{
const IMAGE *im = (IMAGE *) a;
const int sz = n * im->Bands;
double *min = (double *) seq;
int x;
double m;
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: CLOOP( float ); break;
case IM_BANDFMT_DPCOMPLEX: CLOOP( double ); break;
default:
g_assert( 0 );
}
*min = m;
return( 0 );
}
/**
* im_min:
* @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 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().
*
* Returns: 0 on success, -1 on error
*/
int
im_min( IMAGE *in, double *out )
{
double global_min;
if( im_pincheck( in ) ||
im_check_uncoded( "im_min", in ) )
return( -1 );
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 );
/* Back to modulus.
*/
if( im_iscomplex( in ) )
global_min = sqrt( global_min );
*out = global_min;
return( 0 );
}

View File

@ -13,6 +13,8 @@
* - now returns double for value, like im_max()
* 4/9/09
* - gtkdoc comment
* 8/9/09
* - rewrite, from im_maxpos()
*/
/*
@ -41,19 +43,149 @@
*/
/*
#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*/
/* A position and minimum.
*/
typedef struct _Minpos {
int xpos;
int ypos;
double min;
} Minpos;
/* New sequence value.
*/
static void *
minpos_start( IMAGE *in, void *a, void *b )
{
Minpos *global_minpos = (Minpos *) b;
Minpos *minpos;
if( !(minpos = IM_NEW( NULL, Minpos )) )
return( NULL );
*minpos = *global_minpos;
return( (void *) minpos );
}
/* Merge the sequence value back into the per-call state.
*/
static int
minpos_stop( void *seq, void *a, void *b )
{
Minpos *global_minpos = (Minpos *) b;
Minpos *minpos = (Minpos *) seq;
/* Merge.
*/
if( minpos->min > global_minpos->min )
*global_minpos = *minpos;
im_free( seq );
return( 0 );
}
#define LOOP( TYPE ) { \
TYPE *p = (TYPE *) in; \
TYPE m; \
\
m = min; \
\
for( x = 0; x < sz; x++ ) { \
TYPE v = p[x]; \
\
if( v < m ) { \
m = v; \
xpos = r->left + x / reg->im->Bands; \
ypos = r->top + y; \
} \
} \
\
min = 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 < min ) { \
min = mod; \
xpos = r->left + x / reg->im->Bands; \
ypos = r->top + y; \
} \
} \
}
/* Loop over region, adding to seq.
*/
static int
minpos_scan( REGION *reg, void *seq, void *a, void *b )
{
const Rect *r = &reg->valid;
const int sz = IM_REGION_N_ELEMENTS( reg );
Minpos *minpos = (Minpos *) seq;
int x, y;
double min;
int xpos, ypos;
xpos = minpos->xpos;
ypos = minpos->ypos;
min = minpos->min;
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 );
}
}
minpos->xpos = xpos;
minpos->ypos = ypos;
minpos->min = min;
return( 0 );
}
/**
* im_minpos:
* @in: image to search
@ -62,88 +194,69 @@
* @out: returned pixel value at that position
*
* Function to find the minimum of an image. Works for any
* image type. Returns a double and the location of min.
* image type. Returns a double and the location of min. For complex images,
* finds the pixel with the smallest modulus.
*
* This is not a PIO operation! It may use a lot of memory and take a while.
* Needs a rewrite.
*
* See also: im_maxpos(), im_max(), im_stats(), im_maxpos_avg().
* See also: im_maxpos(), im_min(), im_stats(), im_maxpos_avg().
*
* Returns: 0 on success, -1 on error
*/
int
im_minpos( IMAGE *in, int *xpos, int *ypos, double *out )
{
double m;
int xp=0, yp=0;
int os;
Minpos *global_minpos;
/* Check our args. */
if( im_incheck( in ) )
if( im_pincheck( in ) ||
im_check_uncoded( "im_minpos", in ) )
return( -1 );
if( in->Coding != IM_CODING_NONE )
{
im_error("im_minpos", "%s", _("input must be uncoded"));
if( !(global_minpos = IM_NEW( in, Minpos )) )
return( -1 );
}
if( im__value( in, &global_minpos->min ) )
return( -1 );
global_minpos->xpos = 0;
global_minpos->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 ) {
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;
default:
g_assert( 0 );
}
/* Take out bands on x.
/* We use square mod for scanning, for speed.
*/
*out = m;
*xpos = xp / in->Bands;
*ypos = yp;
if( im_iscomplex( in ) )
global_minpos->min *= global_minpos->min;
if( im_iterate( in, minpos_start, minpos_scan, minpos_stop,
in, global_minpos ) )
return( -1 );
/* Back to modulus.
*/
if( im_iscomplex( in ) )
global_minpos->min = sqrt( global_minpos->min );
if( xpos )
*xpos = global_minpos->xpos;
if( ypos )
*ypos = global_minpos->ypos;
if( out )
*out = global_minpos->min;
return( 0 );
}
/**
* im_min:
* @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 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_min(), im_stats().
*
* Returns: 0 on success, -1 on error
*/
int
im_min( IMAGE *in, double *out )
{
return( im_minpos( in, NULL, NULL, out ) );
}