vips_min() done too

This commit is contained in:
John Cupitt 2012-12-05 14:40:01 +00:00
parent 1e332d2f29
commit a4433f1b9f
3 changed files with 261 additions and 129 deletions

View File

@ -24,7 +24,7 @@
- preserve jpeg app13 (photoshop ipct) - preserve jpeg app13 (photoshop ipct)
- nearest neighbour goes back to round down ... round to nearest caused a - nearest neighbour goes back to round down ... round to nearest caused a
range of annoying problems, such as strange half-pixels along edges range of annoying problems, such as strange half-pixels along edges
- vips_max() tracks the top n maxima - vips_max() / _min() track the top n maxima / minima
14/11/12 started 7.30.6 14/11/12 started 7.30.6
- capture tiff warnings earlier - capture tiff warnings earlier

View File

@ -79,7 +79,7 @@
typedef struct _VipsValues { typedef struct _VipsValues {
struct _VipsMax *max; struct _VipsMax *max;
/* The max number of values we track. /* Number of values we track.
*/ */
int size; int size;
@ -88,7 +88,7 @@ typedef struct _VipsValues {
int n; int n;
/* Position and values. We track mod**2 for complex and do a sqrt() at /* Position and values. We track mod**2 for complex and do a sqrt() at
* the end. The three arrays are sorted by values, smallest first. * the end. The three arrays are sorted by @value, smallest first.
*/ */
double *value; double *value;
int *x_pos; int *x_pos;
@ -109,8 +109,8 @@ typedef struct _VipsMax {
int x; int x;
int y; int y;
/* And the postions and values we found as VipsArrays for returning to /* And the positions and values we found as VipsArrays for returning
* our caller. * to our caller.
*/ */
VipsArrayDouble *max_array; VipsArrayDouble *max_array;
VipsArrayInt *x_array; VipsArrayInt *x_array;
@ -143,7 +143,7 @@ vips_values_add( VipsValues *values, double v, int x, int y )
/* Find insertion point. /* Find insertion point.
*/ */
for( i = 0; i < values->n; i++ ) for( i = 0; i < values->n; i++ )
if( values->value[i] > v ) if( v <= values->value[i] )
break; break;
/* Array full? /* Array full?
@ -307,9 +307,8 @@ vips_max_stop( VipsStatistic *statistic, void *seq )
/* float/double max ... no limits, and we have to avoid NaN. /* float/double max ... no limits, and we have to avoid NaN.
* *
* NaN compares false to every float value, so if we were to take the first * NaN compares false to every float value, so we don't need to test for NaN
* point in this buffer as our start max (as we do above) and it was NaN, we'd * in the second loop.
* never replace it with a true value.
*/ */
#define LOOPF( TYPE ) { \ #define LOOPF( TYPE ) { \
TYPE *p = (TYPE *) in; \ TYPE *p = (TYPE *) in; \

View File

@ -10,11 +10,11 @@
* 23/7/93 JC * 23/7/93 JC
* - im_incheck() added * - im_incheck() added
* 20/6/95 JC * 20/6/95 JC
* - now returns double for value, like im_max() * - now returns double for value, like im_min()
* 4/9/09 * 4/9/09
* - gtkdoc comment * - gtkdoc comment
* 8/9/09 * 8/9/09
* - rewrite, from im_maxpos() * - rewrite, from im_minpos()
* 30/8/11 * 30/8/11
* - rewrite as a class * - rewrite as a class
* 5/9/11 * 5/9/11
@ -22,6 +22,9 @@
* 24/2/12 * 24/2/12
* - avoid NaN in float/double/complex images * - avoid NaN in float/double/complex images
* - allow +/- INFINITY as a result * - allow +/- INFINITY as a result
* 4/12/12
* - from min.c
* - track and return bottom n values
*/ */
/* /*
@ -69,21 +72,114 @@
#include "statistic.h" #include "statistic.h"
/* Track min values and position here. We need one of these for each thread,
* and one for the main value.
*
* We will generally only be tracking a small (<10?) number of values, so
* simple arrays will be fastest.
*/
typedef struct _VipsValues {
struct _VipsMin *min;
/* The min number of values we track.
*/
int size;
/* How many values we have in the arrays.
*/
int n;
/* Position and values. We track mod**2 for complex and do a sqrt() at
* the end. The three arrays are sorted by @value, largest first.
*/
double *value;
int *x_pos;
int *y_pos;
} VipsValues;
typedef struct _VipsMin { typedef struct _VipsMin {
VipsStatistic parent_instance; VipsStatistic parent_instance;
gboolean set; /* FALSE means no value yet */ /* Number of values we track.
*/
int size;
/* The current miniumum. When scanning complex images, we keep the /* The single min. Can be unset if, for example, the whole image is
* square of the modulus here and do a single sqrt() right at the end. * NaN.
*/ */
double min; double min;
int x;
int y;
/* And its position. /* And the positions and values we found as VipsArrays for returning
* to our caller.
*/ */
int x, y; VipsArrayDouble *min_array;
VipsArrayInt *x_array;
VipsArrayInt *y_array;
/* Global state here.
*/
VipsValues values;
} VipsMin; } VipsMin;
static void
vips_values_init( VipsValues *values, VipsMin *min )
{
values->min = min;
values->size = min->size;
values->n = 0;
values->value = VIPS_ARRAY( min, values->size, double );
values->x_pos = VIPS_ARRAY( min, values->size, int );
values->y_pos = VIPS_ARRAY( min, values->size, int );
}
/* Add a value. Do nothing if the value is too large.
*/
static void
vips_values_add( VipsValues *values, double v, int x, int y )
{
int i, j;
/* Find insertion point.
*/
for( i = 0; i < values->n; i++ )
if( v >= values->value[i] )
break;
/* Array full?
*/
if( values->n == values->size ) {
if( i > 0 ) {
/* We need to move stuff to the left to make space,
* shunting the largest out.
*/
for( j = 0; j < i - 1; j++ ) {
values->value[j] = values->value[j + 1];
values->x_pos[j] = values->x_pos[j + 1];
values->y_pos[j] = values->y_pos[j + 1];
}
values->value[i - 1] = v;
values->x_pos[i - 1] = x;
values->y_pos[i - 1] = y;
}
}
else {
/* Not full, move stuff to the right into empty space.
*/
for( j = values->n; j > i; j-- ) {
values->value[j] = values->value[j - 1];
values->x_pos[j] = values->x_pos[j - 1];
values->y_pos[j] = values->y_pos[j - 1];
}
values->value[i] = v;
values->x_pos[i] = x;
values->y_pos[i] = y;
values->n += 1;
}
}
typedef VipsStatisticClass VipsMinClass; typedef VipsStatisticClass VipsMinClass;
G_DEFINE_TYPE( VipsMin, vips_min, VIPS_TYPE_STATISTIC ); G_DEFINE_TYPE( VipsMin, vips_min, VIPS_TYPE_STATISTIC );
@ -93,47 +189,72 @@ vips_min_build( VipsObject *object )
{ {
VipsStatistic *statistic = VIPS_STATISTIC( object ); VipsStatistic *statistic = VIPS_STATISTIC( object );
VipsMin *min = (VipsMin *) object; VipsMin *min = (VipsMin *) object;
VipsValues *values = &min->values;
int i;
vips_values_init( values, min );
if( VIPS_OBJECT_CLASS( vips_min_parent_class )->build( object ) ) if( VIPS_OBJECT_CLASS( vips_min_parent_class )->build( object ) )
return( -1 ); return( -1 );
/* For speed we accumulate min ** 2 for complex.
*/
if( vips_bandfmt_iscomplex( vips_image_get_format( statistic->in ) ) )
for( i = 0; i < values->n; i++ )
values->value[i] = sqrt( values->value[i] );
/* Don't set if there's no value (eg. if every pixel is NaN). This /* Don't set if there's no value (eg. if every pixel is NaN). This
* will trigger an error later. * will trigger an error later.
*/ */
if( min->set ) { if( values->n > 0 ) {
double m; VipsArrayDouble *out_array;
VipsArrayInt *x_array;
VipsArrayInt *y_array;
/* For speed we accumulate min^2 for complex. out_array = vips_array_double_new( values->value, values->n );
*/ x_array = vips_array_int_new( values->x_pos, values->n );
m = min->min; y_array = vips_array_int_new( values->y_pos, values->n );
if( vips_bandfmt_iscomplex(
vips_image_get_format( statistic->in ) ) )
m = sqrt( m );
/* We have to set the props via g_object_set() to stop vips /* We have to set the props via g_object_set() to stop vips
* complaining they are unset. * complaining they are unset.
*/ */
g_object_set( min, g_object_set( min,
"out", m, "out", values->value[values->n - 1],
"x", min->x, "x", values->x_pos[values->n - 1],
"y", min->y, "y", values->y_pos[values->n - 1],
"out_array", out_array,
"x_array", x_array,
"y_array", y_array,
NULL ); NULL );
vips_area_unref( (VipsArea *) out_array );
vips_area_unref( (VipsArea *) x_array );
vips_area_unref( (VipsArea *) y_array );
} }
#ifdef DEBUG
printf( "vips_min_build: %d values found\n", values->n );
for( i = 0; i < values->n; i++ )
printf( "%d) %g\t%d\t%d\n",
i,
values->value[i],
values->x_pos[i], values->y_pos[i] );
#endif /*DEBUG*/
return( 0 ); return( 0 );
} }
/* New sequence value. Make a private VipsMin for this thread. /* New sequence value. Make a private VipsValues for this thread.
*/ */
static void * static void *
vips_min_start( VipsStatistic *statistic ) vips_min_start( VipsStatistic *statistic )
{ {
VipsMin *min; VipsValues *values;
min = g_new( VipsMin, 1 ); values = g_new( VipsValues, 1 );
min->set = FALSE; vips_values_init( values, (VipsMin *) statistic );
return( (void *) min ); return( (void *) values );
} }
/* Merge the sequence value back into the per-call state. /* Merge the sequence value back into the per-call state.
@ -141,122 +262,97 @@ vips_min_start( VipsStatistic *statistic )
static int static int
vips_min_stop( VipsStatistic *statistic, void *seq ) vips_min_stop( VipsStatistic *statistic, void *seq )
{ {
VipsMin *global = (VipsMin *) statistic; VipsMin *min = (VipsMin *) statistic;
VipsMin *min = (VipsMin *) seq; VipsValues *values = (VipsValues *) seq;
if( min->set && int i;
(!global->set || min->min < global->min) ) {
global->min = min->min;
global->x = min->x;
global->y = min->y;
global->set = TRUE;
}
g_free( min ); for( i = 0; i < values->n; i++ )
vips_values_add( &min->values,
values->value[i], values->x_pos[i], values->y_pos[i] );
g_free( values );
return( 0 ); return( 0 );
} }
/* real min with a lower bound. /* Real min with a lower bound.
*
* Add values to the buffer if they are less than the buffer maximum. If
* the buffer isn't full, there is no maximum.
*
* Avoid a double test by splitting the loop into two phases: before and after
* the buffer fills.
*
* Stop if our array fills with minval.
*/ */
#define LOOPL( TYPE, LOWER ) { \ #define LOOPU( TYPE, LOWER ) { \
TYPE *p = (TYPE *) in; \ TYPE *p = (TYPE *) in; \
TYPE m; \ TYPE m; \
\ \
if( min->set ) \ for( i = 0; i < sz && values->n < values->size; i++ ) \
m = min->min; \ vips_values_add( values, p[i], x + i / bands, y ); \
else { \ m = values->value[0]; \
m = p[0]; \
min->x = x; \
min->y = y; \
} \
\ \
for( i = 0; i < sz; i++ ) { \ for( ; i < sz; i++ ) { \
if( p[i] < m ) { \ if( p[i] < m ) { \
m = p[i]; \ vips_values_add( values, p[i], x + i / bands, y ); \
min->x = x + i / bands; \ m = values->value[0]; \
min->y = y; \ \
if( m <= LOWER ) { \ if( m <= LOWER ) { \
statistic->stop = TRUE; \ statistic->stop = TRUE; \
break; \ break; \
} \ } \
} \ } \
} \ } \
\
min->min = m; \
min->set = TRUE; \
} }
/* float/double min ... no limits, and we have to avoid NaN. /* float/double min ... no limits, and we have to avoid NaN.
* *
* NaN compares false to every float value, so if we were to take the first * NaN compares false to every float value, so we don't need to test for NaN
* point in this buffer as our start min (as we do above) and it was NaN, we'd * in the second loop.
* never replace it with a true value.
*/ */
#define LOOPF( TYPE ) { \ #define LOOPF( TYPE ) { \
TYPE *p = (TYPE *) in; \ TYPE *p = (TYPE *) in; \
TYPE m; \ TYPE m; \
gboolean set; \
\ \
set = min->set; \ for( i = 0; i < sz && values->n < values->size; i++ ) \
m = min->min; \ if( !isnan( p[i] ) ) \
vips_values_add( values, p[i], x + i / bands, y ); \
m = values->value[0]; \
\ \
for( i = 0; i < sz; i++ ) { \ for( ; i < sz; i++ ) \
if( set ) { \ if( p[i] < m ) { \
if( p[i] < m ) { \ vips_values_add( values, p[i], x + i / bands, y ); \
m = p[i]; \ m = values->value[0]; \
min->x = x + i / bands; \
min->y = y; \
} \
} \ } \
else if( !isnan( p[i] ) ) { \
m = p[i]; \
min->x = x + i / bands; \
min->y = y; \
set = TRUE; \
} \
} \
\
if( set ) { \
min->min = m; \
min->set = TRUE; \
} \
} }
/* As LOOPF, but complex. Track min(mod) to avoid sqrt(). /* As LOOPF, but complex. Track min(mod ** 2) to avoid sqrt().
*/ */
#define LOOPC( TYPE ) { \ #define LOOPC( TYPE ) { \
TYPE *p = (TYPE *) in; \ TYPE *p = (TYPE *) in; \
TYPE m; \ TYPE m; \
gboolean set; \
\ \
set = min->set; \ for( i = 0; i < sz && values->n < values->size; i++ ) { \
m = min->min; \ TYPE mod2 = p[0] * p[0] + p[1] * p[1]; \
\ \
for( i = 0; i < sz; i++ ) { \ if( !isnan( mod2 ) ) \
TYPE mod; \ vips_values_add( values, p[i], x + i / bands, y ); \
\ \
mod = p[0] * p[0] + p[1] * p[1]; \
p += 2; \ p += 2; \
\
if( set ) { \
if( mod > m ) { \
m = mod; \
min->x = x + i / bands; \
min->y = y; \
} \
} \
else if( !isnan( mod ) ) { \
m = mod; \
min->x = x + i / bands; \
min->y = y; \
set = TRUE; \
} \
} \ } \
m = values->value[0]; \
\ \
if( set ) { \ for( ; i < sz; i++ ) { \
min->min = m; \ TYPE mod2 = p[0] * p[0] + p[1] * p[1]; \
min->set = TRUE; \ \
if( mod2 < m ) { \
vips_values_add( values, mod2, x + i / bands, y ); \
m = values->value[0]; \
} \
\
p += 2; \
} \ } \
} }
@ -266,7 +362,7 @@ static int
vips_min_scan( VipsStatistic *statistic, void *seq, vips_min_scan( VipsStatistic *statistic, void *seq,
int x, int y, void *in, int n ) int x, int y, void *in, int n )
{ {
VipsMin *min = (VipsMin *) seq; VipsValues *values = (VipsValues *) seq;
const int bands = vips_image_get_bands( statistic->in ); const int bands = vips_image_get_bands( statistic->in );
const int sz = n * bands; const int sz = n * bands;
@ -274,26 +370,26 @@ vips_min_scan( VipsStatistic *statistic, void *seq,
switch( vips_image_get_format( statistic->in ) ) { switch( vips_image_get_format( statistic->in ) ) {
case VIPS_FORMAT_UCHAR: case VIPS_FORMAT_UCHAR:
LOOPL( unsigned char, 0 ); break; LOOPU( unsigned char, 0 ); break;
case VIPS_FORMAT_CHAR: case VIPS_FORMAT_CHAR:
LOOPL( signed char, SCHAR_MIN ); break; LOOPU( signed char, SCHAR_MIN ); break;
case VIPS_FORMAT_USHORT: case VIPS_FORMAT_USHORT:
LOOPL( unsigned short, 0 ); break; LOOPU( unsigned short, 0 ); break;
case VIPS_FORMAT_SHORT: case VIPS_FORMAT_SHORT:
LOOPL( signed short, SHRT_MIN ); break; LOOPU( signed short, SHRT_MIN ); break;
case VIPS_FORMAT_UINT: case VIPS_FORMAT_UINT:
LOOPL( unsigned int, 0 ); break; LOOPU( unsigned int, 0 ); break;
case VIPS_FORMAT_INT: case VIPS_FORMAT_INT:
LOOPL( signed int, INT_MIN ); break; LOOPU( signed int, INT_MIN ); break;
case VIPS_FORMAT_FLOAT: case VIPS_FORMAT_FLOAT:
LOOPF( float ); break; LOOPF( float ); break;
case VIPS_FORMAT_DOUBLE: case VIPS_FORMAT_DOUBLE:
LOOPF( double ); break; LOOPF( double ); break;
case VIPS_FORMAT_COMPLEX: case VIPS_FORMAT_COMPLEX:
LOOPC( float ); break; LOOPC( float ); break;
case VIPS_FORMAT_DPCOMPLEX: case VIPS_FORMAT_DPCOMPLEX:
LOOPC( double ); break; LOOPC( double ); break;
default: default:
@ -335,34 +431,67 @@ vips_min_class_init( VipsMinClass *class )
G_STRUCT_OFFSET( VipsMin, x ), G_STRUCT_OFFSET( VipsMin, x ),
0, 1000000, 0 ); 0, 1000000, 0 );
VIPS_ARG_INT( class, "y", 2, VIPS_ARG_INT( class, "y", 3,
_( "y" ), _( "y" ),
_( "Vertical position of minimum" ), _( "Vertical position of minimum" ),
VIPS_ARGUMENT_OPTIONAL_OUTPUT, VIPS_ARGUMENT_OPTIONAL_OUTPUT,
G_STRUCT_OFFSET( VipsMin, y ), G_STRUCT_OFFSET( VipsMin, y ),
0, 1000000, 0 ); 0, 1000000, 0 );
VIPS_ARG_INT( class, "size", 4,
_( "Size" ),
_( "Number of minimum values to find" ),
VIPS_ARGUMENT_OPTIONAL_INPUT,
G_STRUCT_OFFSET( VipsMin, size ),
0, 1000000, 10 );
VIPS_ARG_BOXED( class, "out_array", 6,
_( "Output array" ),
_( "Array of output values" ),
VIPS_ARGUMENT_OPTIONAL_OUTPUT,
G_STRUCT_OFFSET( VipsMin, min_array ),
VIPS_TYPE_ARRAY_DOUBLE );
VIPS_ARG_BOXED( class, "x_array", 7,
_( "x array" ),
_( "Array of horizontal positions" ),
VIPS_ARGUMENT_OPTIONAL_OUTPUT,
G_STRUCT_OFFSET( VipsMin, x_array ),
VIPS_TYPE_ARRAY_INT );
VIPS_ARG_BOXED( class, "y_array", 8,
_( "y array" ),
_( "Array of vertical positions" ),
VIPS_ARGUMENT_OPTIONAL_OUTPUT,
G_STRUCT_OFFSET( VipsMin, y_array ),
VIPS_TYPE_ARRAY_INT );
} }
static void static void
vips_min_init( VipsMin *min ) vips_min_init( VipsMin *min )
{ {
min->size = 1;
} }
/** /**
* vips_min: * vips_min:
* @in: input #VipsImage * @in: input #VipsImage
* @out: output pixel maximum * @out: output pixel minimum
* @...: %NULL-terminated list of optional named arguments * @...: %NULL-terminated list of optional named arguments
* *
* Optional arguments: * Optional arguments:
* *
* @x: horizontal position of minimum * @x: horizontal position of minimum
* @y: vertical position of minimum * @y: vertical position of minimum
* @size: number of minima to find
* @out_array: return array of minimum values
* @x_array: corresponding horizontal positions
* @y_array: corresponding vertical positions
* *
* This operation finds the minimum value in an image. * This operation finds the minimum value in an image.
* *
* If the image contains several minimum values, only the first one found is * If the image contains several minimum values, only the first @size
* returned. * found are returned.
* *
* It operates on all * It operates on all
* bands of the input image: use vips_stats() if you need to find an * bands of the input image: use vips_stats() if you need to find an
@ -370,7 +499,11 @@ vips_min_init( VipsMin *min )
* *
* For complex images, this operation finds the minimum modulus. * For complex images, this operation finds the minimum modulus.
* *
* See also: vips_max(), vips_stats(). * You can read out the position of the minimum with @x and @y. You can read
* out arrays of the values and positions of the top @size minima with
* @out_array, @x_array and @y_array.
*
* See also: vips_min(), vips_stats().
* *
* Returns: 0 on success, -1 on error * Returns: 0 on success, -1 on error
*/ */