This commit is contained in:
John Cupitt 2012-12-04 14:57:41 +00:00
parent 7a202cf95e
commit 713ecf8252
4 changed files with 62 additions and 560 deletions

2
TODO
View File

@ -1,3 +1,5 @@
- im_maxpos_avg() should be NaN-avoiding?
- now we've removed round-to-nearest from NN, we need something extra in the
affine transform to displace the input cods

View File

@ -9,7 +9,6 @@ libarithmetic_la_SOURCES = \
divide.c \
im_linreg.c \
im_maxpos_avg.c \
im_maxpos_vec.c \
measure.c \
multiply.c \
im_point_bilinear.c \

View File

@ -1,492 +0,0 @@
/* im_maxpos_vec.c
*
* Copyright: 2006, The Nottingham Trent University
* Copyright: 2006, Tom Vajzovic
*
* Author: Tom Vajzovic
*
* Written on: 2006-09-01
*
* 9/9/09
* - gtkdoc comments
*/
/*
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
*/
/** HEADERS **/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif /* HAVE_CONFIG_H */
#include <vips/intl.h>
#include <stdlib.h>
#include <float.h>
#include <vips/vips.h>
/** TYPE DEFINITIONS **/
typedef struct {
int *xs;
int *ys;
double *vals;
int *ptrs;
int start;
} maxpos_list;
/** LOCAL FUNCTIONS DECLARATIONS **/
static maxpos_list *maxpos_list_alloc( int n );
static void maxpos_list_free( maxpos_list *list );
static void maxpos_list_init( maxpos_list *list, int n );
static void *maxpos_vec_start( IMAGE *unrequired, void *, void * );
static int maxpos_vec_scan( REGION *reg, void *seq, void *, void *, gboolean * );
static void add_to_maxpos_list( maxpos_list *list, int x, int y, double val );
static int maxpos_vec_stop( void *seq, void *, void * );
static void minpos_list_init( maxpos_list *list, int n );
static void *minpos_vec_start( IMAGE *unrequired, void *, void * );
static int minpos_vec_scan( REGION *reg, void *seq, void *, void *, gboolean * );
static void add_to_minpos_list( maxpos_list *list, int x, int y, double val );
static int minpos_vec_stop( void *seq, void *, void * );
/** EXPORTED FUNCTIONS **/
/**
* im_maxpos_vec:
* @im: image to scan
* @xpos: array to return x positions
* @ypos: array to return y positions
* @maxima: array to return values
* @n: number of maxima to search for
*
* Find the coordinates and values of the @n @maxima of an image.
*
* For 8 and 16-bit images, it's much faster to find the histogram and then
* calculate a threshold from that. See im_mpercent().
*
* See also: im_minpos(), im_min(), im_stats(), im_maxpos_avg().
*
* Returns: 0 on success, -1 on error
*/
int im_maxpos_vec( IMAGE *im, int *xpos, int *ypos, double *maxima, int n ){
#define FUNCTION_NAME "im_maxpos_vec"
/* number of sequences used is beyond my control at this level, but I note that */
/* efficiency decreases as more sequences are used - speed may still increase */
int result;
int *pointers= im_malloc( NULL, n * sizeof( int* ) );
maxpos_list master_list= { xpos, ypos, maxima, pointers, 0 };
if( im_pincheck( im ) )
return -1;
if( !pointers )
return -1;
if( ! ( vips_bandfmt_isint( im->BandFmt ) ||
vips_bandfmt_isfloat( im->BandFmt ) ) ){
im_error( FUNCTION_NAME, "%s", _( "scalar images only" ) );
return -1;
}
if( 1 != im-> Bands ){
im_error( FUNCTION_NAME, "%s", _( "single band images only" ) );
return -1;
}
if( IM_CODING_NONE != im-> Coding ){
im_error( FUNCTION_NAME, "%s", _( "uncoded images only" ) );
return -1;
}
if( ! xpos || ! ypos || ! maxima || n < 1 ){
im_error( FUNCTION_NAME, "%s", _( "invalid argument" ) );
return -1;
}
maxpos_list_init( &master_list, n );
result= vips_sink( im, maxpos_vec_start, maxpos_vec_scan, maxpos_vec_stop, &n, &master_list );
im_free( pointers );
return result;
#undef FUNCTION_NAME
}
/**
* im_minpos_vec:
* @im: image to scan
* @xpos: array to return x positions
* @ypos: array to return y positions
* @minima: array to return values
* @n: number of minima to search for
*
* Find the coordinates and values of the @n @minima of an image.
*
* For 8 and 16-bit images, it's much faster to find the histogram and then
* calculate a threshold from that. See im_mpercent().
*
* See also: im_minpos(), im_min(), im_stats(), im_maxpos_avg().
*
* Returns: 0 on success, -1 on error
*/
int im_minpos_vec( IMAGE *im, int *xpos, int *ypos, double *minima, int n ){
#define FUNCTION_NAME "im_minpos_vec"
/* number of sequences used is beyond my control at this level, but I note that */
/* effeciency decreases as more sequences are used - speed may still increase */
int result;
int *pointers= im_malloc( NULL, n * sizeof( int* ) );
maxpos_list master_list= { xpos, ypos, minima, pointers, 0 };
if( im_pincheck( im ) )
return -1;
if( !pointers )
return -1;
if( ! ( vips_bandfmt_isint( im->BandFmt ) ||
vips_bandfmt_isfloat( im->BandFmt ) ) ){
im_error( FUNCTION_NAME, "%s", _( "scalar images only" ) );
return -1;
}
if( 1 != im-> Bands ){
im_error( FUNCTION_NAME, "%s", _( "single band images only" ) );
return -1;
}
if( IM_CODING_NONE != im-> Coding ){
im_error( FUNCTION_NAME, "%s", _( "uncoded images only" ) );
return -1;
}
if( ! xpos || ! ypos || ! minima || n < 1 ){
im_error( FUNCTION_NAME, "%s", _( "invalid argument" ) );
return -1;
}
minpos_list_init( &master_list, n );
result= vips_sink( im, minpos_vec_start, minpos_vec_scan, minpos_vec_stop, &n, &master_list );
im_free( pointers );
return result;
#undef FUNCTION_NAME
}
/** LOCAL FUNCTION DEFINITIONS **/
static maxpos_list *maxpos_list_alloc( int n ){
maxpos_list *list= im_malloc( NULL, sizeof( maxpos_list ) );
if( ! list )
return NULL;
list-> xs= im_malloc( NULL, 3 * n * sizeof( int ) );
list-> vals= im_malloc( NULL, n * sizeof( double ) );
if( ! list-> xs || ! list-> vals ){
im_free( list-> xs );
im_free( list-> vals );
im_free( list );
return NULL;
}
list-> ys= list-> xs + n;
list-> ptrs= list-> ys + n;
return list;
}
static void maxpos_list_free( maxpos_list *list ){
im_free( list-> xs );
im_free( list-> vals );
im_free( list );
}
static void maxpos_list_init( maxpos_list *list, int n ){
int i;
for( i= 0; i < n; ++i ){
list-> xs[ i ]= 0;
list-> ys[ i ]= 0;
list-> vals[ i ]= 0;
list-> ptrs[ i ]= i + 1;
}
list-> ptrs[ n - 1 ]= -1;
list-> start= 0;
}
static void *maxpos_vec_start( IMAGE *unrequired, void *a, void *b ){
int *n = (int *) a;
maxpos_list *list= maxpos_list_alloc( *n );
if( ! list )
return NULL;
maxpos_list_init( list, *n );
return list;
}
static int maxpos_vec_scan( REGION *reg, void *seq, void *a, void *b, gboolean *stop ){
maxpos_list *list = (maxpos_list *) seq;
#define MAXPOS_VEC_SCAN( type ){ \
\
int y= reg-> valid. top; \
int x; \
int ymax= y + reg-> valid. height; \
int xmax= reg-> valid. left + reg-> valid. width; \
\
type *row= (type*)IM_REGION_ADDR( reg, reg-> valid. left, y ) - reg-> valid. left; \
size_t skip= IM_REGION_LSKIP( reg ) / sizeof( type ); \
\
for( ; y < ymax; ++y, row+= skip ) \
for( x= reg-> valid. left; x < xmax; ++x ) \
if( row[ x ] > list-> vals[ list-> start ] ) \
add_to_maxpos_list( list, x, y, row[ x ] ); \
}
switch( reg-> im-> BandFmt ){
case IM_BANDFMT_UCHAR: MAXPOS_VEC_SCAN( guint8 ) break;
case IM_BANDFMT_CHAR: MAXPOS_VEC_SCAN( gint8 ) break;
case IM_BANDFMT_USHORT: MAXPOS_VEC_SCAN( guint16 ) break;
case IM_BANDFMT_SHORT: MAXPOS_VEC_SCAN( gint16 ) break;
case IM_BANDFMT_UINT: MAXPOS_VEC_SCAN( guint32 ) break;
case IM_BANDFMT_INT: MAXPOS_VEC_SCAN( gint32 ) break;
case IM_BANDFMT_FLOAT: MAXPOS_VEC_SCAN( float ) break;
case IM_BANDFMT_DOUBLE: MAXPOS_VEC_SCAN( double ) break;
default:
g_assert( 0 );
}
#undef MAXPOS_VEC_SCAN
return 0;
}
static void add_to_maxpos_list( maxpos_list *list, int x, int y, double val ){
int pointer= list-> start;
while( -1 != list-> ptrs[ pointer ] && val > list-> vals[ list-> ptrs[ pointer ] ] )
pointer= list-> ptrs[ pointer ];
list-> xs[ list-> start ]= x;
list-> ys[ list-> start ]= y;
list-> vals[ list-> start ]= val;
if( list-> start != pointer ){
/* we are adding mid-chain not at the very bottom */
int second= list-> ptrs[ list-> start ];
list-> ptrs[ list-> start ]= list-> ptrs[ pointer ];
list-> ptrs[ pointer ]= list-> start;
list-> start= second;
}
}
static int maxpos_vec_stop( void *seq, void *a, void *b ){
/* reverse list */
maxpos_list *list = (maxpos_list *) seq;
maxpos_list *master_list = (maxpos_list *) b;
int prev= -1;
int pointer= list-> start;
while( -1 != list-> ptrs[ pointer ] ){
int next= list-> ptrs[ pointer ];
list-> ptrs[ pointer ]= prev;
prev= pointer;
pointer= next;
}
list-> ptrs[ pointer ]= prev;
list-> start= pointer;
/* add to main list */
for( ; -1 != pointer; pointer= list-> ptrs[ pointer ] )
/* loop over all the ones found in this sequence */
if( list-> vals[ pointer ] > master_list-> vals[ master_list-> start ] )
add_to_maxpos_list( master_list, list-> xs[ pointer ], list-> ys[ pointer ], list-> vals[ pointer ] );
else
break;
/* since we are now high->low, if this one isn't big enough, none of the rest are */
maxpos_list_free( list );
return 0;
}
static void minpos_list_init( maxpos_list *list, int n ){
int i;
for( i= 0; i < n; ++i ){
list-> xs[ i ]= 0;
list-> ys[ i ]= 0;
list-> vals[ i ]= DBL_MAX;
list-> ptrs[ i ]= i + 1;
}
list-> ptrs[ n - 1 ]= -1;
list-> start= 0;
}
static void *minpos_vec_start( IMAGE *unrequired, void *a, void *b ){
int *n = (int *) a;
maxpos_list *list= maxpos_list_alloc( *n );
if( ! list )
return NULL;
minpos_list_init( list, *n );
return list;
}
static int minpos_vec_scan( REGION *reg, void *seq, void *a, void *b, gboolean *stop ){
maxpos_list *list = (maxpos_list *) seq;
#define MINPOS_VEC_SCAN( type ){ \
\
int y= reg-> valid. top; \
int x; \
int ymax= y + reg-> valid. height; \
int xmax= reg-> valid. left + reg-> valid. width; \
\
type *row= (type*)IM_REGION_ADDR( reg, reg-> valid. left, y ) - reg-> valid. left; \
size_t skip= IM_REGION_LSKIP( reg ) / sizeof( type ); \
\
for( ; y < ymax; ++y, row+= skip ) \
for( x= reg-> valid. left; x < xmax; ++x ) \
if( row[ x ] < list-> vals[ list-> start ] ) \
add_to_minpos_list( list, x, y, row[ x ] ); \
}
switch( reg-> im-> BandFmt ){
case IM_BANDFMT_UCHAR: MINPOS_VEC_SCAN( guint8 ) break;
case IM_BANDFMT_CHAR: MINPOS_VEC_SCAN( gint8 ) break;
case IM_BANDFMT_USHORT: MINPOS_VEC_SCAN( guint16 ) break;
case IM_BANDFMT_SHORT: MINPOS_VEC_SCAN( gint16 ) break;
case IM_BANDFMT_UINT: MINPOS_VEC_SCAN( guint32 ) break;
case IM_BANDFMT_INT: MINPOS_VEC_SCAN( gint32 ) break;
case IM_BANDFMT_FLOAT: MINPOS_VEC_SCAN( float ) break;
case IM_BANDFMT_DOUBLE: MINPOS_VEC_SCAN( double ) break;
default:
g_assert( 0 );
}
#undef MINPOS_VEC_SCAN
return 0;
}
static void add_to_minpos_list( maxpos_list *list, int x, int y, double val ){
int pointer= list-> start;
while( -1 != list-> ptrs[ pointer ] && val < list-> vals[ list-> ptrs[ pointer ] ] )
pointer= list-> ptrs[ pointer ];
list-> xs[ list-> start ]= x;
list-> ys[ list-> start ]= y;
list-> vals[ list-> start ]= val;
if( list-> start != pointer ){
/* we are adding mid-chain not at the very bottom */
int second= list-> ptrs[ list-> start ];
list-> ptrs[ list-> start ]= list-> ptrs[ pointer ];
list-> ptrs[ pointer ]= list-> start;
list-> start= second;
}
}
static int minpos_vec_stop( void *seq, void *a, void *b ){
/* reverse list */
maxpos_list *list = (maxpos_list *) seq;
maxpos_list *master_list = (maxpos_list *) b;
int prev= -1;
int pointer= list-> start;
while( -1 != list-> ptrs[ pointer ] ){
int next= list-> ptrs[ pointer ];
list-> ptrs[ pointer ]= prev;
prev= pointer;
pointer= next;
}
list-> ptrs[ pointer ]= prev;
list-> start= pointer;
/* add to main list */
for( ; -1 != pointer; pointer= list-> ptrs[ pointer ] )
/* loop over all the ones found in this sequence */
if( list-> vals[ pointer ] < master_list-> vals[ master_list-> start ] )
add_to_minpos_list( master_list, list-> xs[ pointer ], list-> ys[ pointer ], list-> vals[ pointer ] );
else
break;
/* since we are now high->low, if this one isn't big enough, none of the rest are */
maxpos_list_free( list );
return 0;
}

View File

@ -21,6 +21,8 @@
* 24/2/12
* - avoid NaN in float/double/complex images
* - allow +/- INFINITY as a result
* 4/12/12
* - track top n values
*/
/*
@ -98,7 +100,14 @@ typedef struct _VipsMax {
/* Max number of values we track.
*/
int max_values;
int size;
/* The single max. Can be unset if, for example, the whole image is
* NaN.
*/
double max;
int x;
int y;
/* Global state here.
*/
@ -110,7 +119,7 @@ vips_values_init( VipsValues *values, VipsMax *max )
{
values->max = max;
values->size = max->max_values;
values->size = max->size;
values->n = 0;
values->value = VIPS_ARRAY( max, values->size, double );
values->x_pos = VIPS_ARRAY( max, values->size, int );
@ -171,32 +180,31 @@ vips_max_build( VipsObject *object )
{
VipsStatistic *statistic = VIPS_STATISTIC( object );
VipsMax *max = (VipsMax *) object;
VipsValues *values = &max->values;
int i;
vips_values_init( *max->values, max );
vips_values_init( values, max );
if( VIPS_OBJECT_CLASS( vips_max_parent_class )->build( object ) )
return( -1 );
/* For speed we accumulate max ** 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
* will trigger an error later.
*/
if( max->set ) {
double m;
/* For speed we accumulate max^2 for complex.
*/
m = max->max;
if( vips_bandfmt_iscomplex(
vips_image_get_format( statistic->in ) ) )
m = sqrt( m );
if( values->n > 0 ) {
/* We have to set the props via g_object_set() to stop vips
* complaining they are unset.
*/
g_object_set( max,
"out", m,
"x", max->x,
"y", max->y,
"out", values->value[values->n - 1],
"x", values->x_pos[values->n - 1],
"y", values->y_pos[values->n - 1],
NULL );
}
@ -235,34 +243,29 @@ vips_max_stop( VipsStatistic *statistic, void *seq )
return( 0 );
}
/* real max with an upper bound.
/* Real max with an upper bound. Stop if our array fills with maxval.
*/
#define LOOPU( TYPE, UPPER ) { \
TYPE *p = (TYPE *) in; \
TYPE m; \
\
if( max->set ) \
m = max->max; \
if( values->n ) \
m = values->value[0]; \
else { \
m = p[0]; \
max->x = x; \
max->y = y; \
} \
\
for( i = 0; i < sz; i++ ) { \
if( p[i] > m ) { \
m = p[i]; \
max->x = x + i / bands; \
max->y = y; \
vips_values_add( values, p[i], x + i / bands, y ); \
m = values->value[0]; \
\
if( m >= UPPER ) { \
statistic->stop = TRUE; \
break; \
} \
} \
} \
\
max->max = m; \
max->set = TRUE; \
}
/* float/double max ... no limits, and we have to avoid NaN.
@ -276,29 +279,18 @@ vips_max_stop( VipsStatistic *statistic, void *seq )
TYPE m; \
gboolean set; \
\
set = max->set; \
m = max->max; \
set = values->n > 0; \
if( set ) \
m = values->value[0]; \
\
for( i = 0; i < sz; i++ ) { \
if( set ) { \
if( p[i] > m ) { \
m = p[i]; \
max->x = x + i / bands; \
max->y = y; \
} \
} \
else if( !isnan( p[i] ) ) { \
m = p[i]; \
max->x = x + i / bands; \
max->y = y; \
if( (set && p[i] > m) || \
(!set && !isnan( p[i] )) ) { \
vips_values_add( values, p[i], x + i / bands, y ); \
m = values->value[0]; \
set = TRUE; \
} \
} \
\
if( set ) { \
max->max = m; \
max->set = TRUE; \
} \
}
/* As LOOPF, but complex. Track max(mod) to avoid sqrt().
@ -308,8 +300,9 @@ vips_max_stop( VipsStatistic *statistic, void *seq )
TYPE m; \
gboolean set; \
\
set = max->set; \
m = max->max; \
set = values->n > 0; \
if( set ) \
m = values->value[0]; \
\
for( i = 0; i < sz; i++ ) { \
TYPE mod; \
@ -317,25 +310,13 @@ vips_max_stop( VipsStatistic *statistic, void *seq )
mod = p[0] * p[0] + p[1] * p[1]; \
p += 2; \
\
if( set ) { \
if( mod > m ) { \
m = mod; \
max->x = x + i / bands; \
max->y = y; \
} \
} \
else if( !isnan( mod ) ) { \
m = mod; \
max->x = x + i / bands; \
max->y = y; \
if( (set && mod > m) || \
(!set && !isnan( mod )) ) { \
vips_values_add( values, mod, x + i / bands, y ); \
m = values->value[0]; \
set = TRUE; \
} \
} \
\
if( set ) { \
max->max = m; \
max->set = TRUE; \
} \
}
/* Loop over region, adding to seq.
@ -344,7 +325,7 @@ static int
vips_max_scan( VipsStatistic *statistic, void *seq,
int x, int y, void *in, int n )
{
VipsMax *max = (VipsMax *) seq;
VipsValues *values = (VipsValues *) seq;
const int bands = vips_image_get_bands( statistic->in );
const int sz = n * bands;
@ -413,17 +394,26 @@ vips_max_class_init( VipsMaxClass *class )
G_STRUCT_OFFSET( VipsMax, x ),
0, 1000000, 0 );
VIPS_ARG_INT( class, "y", 2,
VIPS_ARG_INT( class, "y", 3,
_( "y" ),
_( "Vertical position of maximum" ),
VIPS_ARGUMENT_OPTIONAL_OUTPUT,
G_STRUCT_OFFSET( VipsMax, y ),
0, 1000000, 0 );
VIPS_ARG_INT( class, "size", 4,
_( "Size" ),
_( "Number of maximum values to find" ),
VIPS_ARGUMENT_OPTIONAL_INPUT,
G_STRUCT_OFFSET( VipsMax, size ),
0, 1000000, 10 );
}
static void
vips_max_init( VipsMax *max )
{
max->size = 10;
}
/**
@ -436,11 +426,12 @@ vips_max_init( VipsMax *max )
*
* @x: horizontal position of maximum
* @y: vertical position of maximum
* @size: number of maxima to find
*
* This operation finds the maximum value in an image.
*
* If the image contains several maximum values, only the first one found is
* returned.
* If the image contains several maximum values, only the first @size
* found are returned.
*
* It operates on all
* bands of the input image: use vips_stats() if you need to find an
@ -448,6 +439,8 @@ vips_max_init( VipsMax *max )
*
* For complex images, this operation finds the maximum modulus.
*
* You can read out the position of the maximum with @x and @y.
*
* See also: vips_min(), vips_stats().
*
* Returns: 0 on success, -1 on error