redo im_deviate() as a class

This commit is contained in:
John Cupitt 2011-11-06 17:59:07 +00:00
parent 77d386cdaf
commit 1fd1559bb5
10 changed files with 352 additions and 285 deletions

View File

@ -7,7 +7,7 @@
im_lintra(), im_lintra_vec(), im_black(), im_rot90, im_rot180(), im_rot270()
im_sintra(), im_costra(), im_tantra(), im_asintra(), im_acostra(),
im_atantra(), im_exptra(), im_exp10tra(), im_logtra(), im_log10tra(),
im_abs(), im_sign(), im_max(), im_maxpos()
im_abs(), im_sign(), im_max(), im_maxpos(), im_deviate()
redone as classes
- added argument priorites to help control arg ordering
- generate has a 'stop' param to signal successful early termination

View File

@ -5,7 +5,7 @@ libarithmetic_la_SOURCES = \
abs.c \
im_bandmean.c \
im_cross_phase.c \
im_deviate.c \
deviate.c \
im_divide.c \
im_recomb.c \
im_linreg.c \

View File

@ -507,9 +507,11 @@ vips_arithmetic_operation_init( void )
extern GType vips_avg_get_type( void );
extern GType vips_min_get_type( void );
extern GType vips_max_get_type( void );
extern GType vips_deviate_get_type( void );
extern GType vips_linear_get_type( void );
extern GType vips_math_get_type( void );
extern GType vips_abs_get_type( void );
extern GType vips_sign_get_type( void );
vips_add_get_type();
vips_invert_get_type();
@ -517,8 +519,10 @@ vips_arithmetic_operation_init( void )
vips_avg_get_type();
vips_min_get_type();
vips_max_get_type();
vips_deviate_get_type();
vips_linear_get_type();
vips_math_get_type();
vips_abs_get_type();
vips_sign_get_type();
}

View File

@ -90,7 +90,7 @@ typedef struct _VipsAvg {
VipsStatistic parent_instance;
double sum;
double avg;
double out;
} VipsAvg;
typedef VipsStatisticClass VipsAvgClass;
@ -103,7 +103,7 @@ vips_avg_build( VipsObject *object )
VipsStatistic *statistic = VIPS_STATISTIC( object );
VipsAvg *avg = (VipsAvg *) object;
gint64 vals, pels;
gint64 vals;
double average;
if( VIPS_OBJECT_CLASS( vips_avg_parent_class )->build( object ) )
@ -112,9 +112,10 @@ vips_avg_build( VipsObject *object )
/* Calculate average. For complex, we accumulate re*re +
* im*im, so we need to sqrt.
*/
pels = (gint64) vips_image_get_width( statistic->in ) *
vips_image_get_height( statistic->in );
vals = pels * vips_image_get_bands( statistic->in );
vals = (gint64)
vips_image_get_width( statistic->in ) *
vips_image_get_height( statistic->in ) *
vips_image_get_bands( statistic->in );
average = avg->sum / vals;
if( vips_bandfmt_iscomplex( vips_image_get_format( statistic->in ) ) )
average = sqrt( average );
@ -234,7 +235,7 @@ vips_avg_class_init( VipsAvgClass *class )
_( "Output" ),
_( "Output value" ),
VIPS_ARGUMENT_REQUIRED_OUTPUT,
G_STRUCT_OFFSET( VipsAvg, avg ),
G_STRUCT_OFFSET( VipsAvg, out ),
-G_MAXDOUBLE, G_MAXDOUBLE, 0.0 );
}

View File

@ -0,0 +1,262 @@
/* VipsDeviate
*
* Copyright: 1990, J. Cupitt
*
* Author: J. Cupitt
* Written on: 02/08/1990
* Modified on:
* 5/5/93 JC
* - now does partial images
* - less likely to overflow
* - adapted from im_deviate
* 1/7/93 JC
* - adapted for partial v2
* - ANSIfied
* 21/2/95 JC
* - modernised again
* 11/5/95 JC
* - oops! return( NULL ) in im_deviate(), instead of return( -1 )
* 20/6/95 JC
* - now returns double, not float
* 13/1/05
* - use 64 bit arithmetic
* 8/12/06
* - add liboil support
* 2/9/09
* - gtk-doc comment
* - minor reformatting
* 4/9/09
* - use im__wrapscan()
* 31/7/10
* - remove liboil
* 6/11/11
* - rewrite as a class
*/
/*
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
*/
#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>
#include "statistic.h"
/**
* VipsDeviate:
* @in: input #IMAGE
* @out: output pixel standard deviation
*
* This operation finds the standard deviation of all pixels in @in. It
* operates on all bands of the input image: use im_stats() if you need
* to calculate an average for each band.
*
* Non-complex images only.
*
* See also: VipsDeviate.
*/
typedef struct _VipsDeviate {
VipsStatistic parent_instance;
double sum;
double sum2;
double out;
} VipsDeviate;
typedef VipsStatisticClass VipsDeviateClass;
G_DEFINE_TYPE( VipsDeviate, vips_deviate, VIPS_TYPE_STATISTIC );
static int
vips_deviate_build( VipsObject *object )
{
VipsStatistic *statistic = VIPS_STATISTIC( object );
VipsDeviate *deviate = (VipsDeviate *) object;
gint64 vals;
double s, s2;
if( statistic->in &&
vips_check_noncomplex( "VipsDeviate", statistic->in ) )
return( -1 );
if( VIPS_OBJECT_CLASS( vips_deviate_parent_class )->build( object ) )
return( -1 );
/*
NOTE: NR suggests a two-pass algorithm to minimise roundoff.
But that's too expensive for us :-( so do it the old one-pass
way.
*/
/* Calculate and return deviation. Add a fabs to stop sqrt(<=0).
*/
vals = (gint64)
vips_image_get_width( statistic->in ) *
vips_image_get_height( statistic->in ) *
vips_image_get_bands( statistic->in );
s = deviate->sum;
s2 = deviate->sum2;
g_object_set( object,
"out", sqrt( fabs( s2 - (s * s / vals) ) / (vals - 1) ),
NULL );
return( 0 );
}
/* Start function: allocate space for an array in which we can accumulate the
* sum and sum of squares for this thread.
*/
static void *
vips_deviate_start( VipsStatistic *statistic )
{
double *ss2;
ss2 = g_new( double, 2 );
ss2[0] = 0.0;
ss2[1] = 0.0;
return( (void *) ss2 );
}
/* Stop function. Add this little sum to the main sum.
*/
static int
vips_deviate_stop( VipsStatistic *statistic, void *seq )
{
VipsDeviate *deviate = (VipsDeviate *) statistic;
double *ss2 = (double *) seq;
deviate->sum += ss2[0];
deviate->sum2 += ss2[1];
g_free( ss2 );
return( 0 );
}
#define LOOP( TYPE ) { \
TYPE *p = (TYPE *) in; \
\
for( x = 0; x < sz; x++ ) { \
TYPE v = p[x]; \
\
sum += v; \
sum2 += (double) v * (double) v; \
} \
}
static int
vips_deviate_scan( VipsStatistic *statistic, void *seq,
int x, int y, void *in, int n )
{
const int sz = n * vips_image_get_bands( statistic->in );
double *ss2 = (double *) seq;
double sum;
double sum2;
sum = ss2[0];
sum2 = ss2[1];
/* Now generate code for all types.
*/
switch( vips_image_get_format( statistic->in ) ) {
case VIPS_FORMAT_UCHAR: LOOP( unsigned char ); break;
case VIPS_FORMAT_CHAR: LOOP( signed char ); break;
case VIPS_FORMAT_USHORT: LOOP( unsigned short ); break;
case VIPS_FORMAT_SHORT: LOOP( signed short ); break;
case VIPS_FORMAT_UINT: LOOP( unsigned int ); break;
case VIPS_FORMAT_INT: LOOP( signed int ); break;
case VIPS_FORMAT_FLOAT: LOOP( float ); break;
case VIPS_FORMAT_DOUBLE: LOOP( double ); break;
default:
g_assert( 0 );
}
ss2[0] = sum;
ss2[1] = sum2;
return( 0 );
}
static void
vips_deviate_class_init( VipsDeviateClass *class )
{
GObjectClass *gobject_class = (GObjectClass *) class;
VipsObjectClass *object_class = (VipsObjectClass *) class;
VipsStatisticClass *sclass = VIPS_STATISTIC_CLASS( class );
gobject_class->set_property = vips_object_set_property;
gobject_class->get_property = vips_object_get_property;
object_class->nickname = "deviate";
object_class->description = _( "find image average" );
object_class->build = vips_deviate_build;
sclass->start = vips_deviate_start;
sclass->scan = vips_deviate_scan;
sclass->stop = vips_deviate_stop;
VIPS_ARG_DOUBLE( class, "out", 2,
_( "Output" ),
_( "Output value" ),
VIPS_ARGUMENT_REQUIRED_OUTPUT,
G_STRUCT_OFFSET( VipsDeviate, out ),
-G_MAXDOUBLE, G_MAXDOUBLE, 0.0 );
}
static void
vips_deviate_init( VipsDeviate *deviate )
{
}
int
vips_deviate( VipsImage *in, double *out, ... )
{
va_list ap;
int result;
va_start( ap, out );
result = vips_call_split( "deviate", ap, in, out );
va_end( ap );
return( result );
}

View File

@ -1,276 +0,0 @@
/* im_deviate.c
*
* Copyright: 1990, J. Cupitt
*
* Author: J. Cupitt
* Written on: 02/08/1990
* Modified on:
* 5/5/93 JC
* - now does partial images
* - less likely to overflow
* - adapted from im_avg
* 1/7/93 JC
* - adapted for partial v2
* - ANSIfied
* 21/2/95 JC
* - modernised again
* 11/5/95 JC
* - oops! return( NULL ) in im_avg(), instead of return( -1 )
* 20/6/95 JC
* - now returns double, not float
* 13/1/05
* - use 64 bit arithmetic
* 8/12/06
* - add liboil support
* 2/9/09
* - gtk-doc comment
* - minor reformatting
* 4/9/09
* - use im__wrapscan()
* 31/7/10
* - remove liboil
*/
/*
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
*/
#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>
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( VipsImage *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, gboolean *stop )
{
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 vips_sink(), but the scan function works a line at a time, like
* im_wrap*(). Shared with im_min(), im_deviate() etc.
*/
int
im__wrapscan( VipsImage *in,
VipsStartFn start, im__wrapscan_fn scan, VipsStopFn 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( vips_sink( in,
wrapscan_start, wrapscan_scan, wrapscan_stop,
&wrapscan, NULL ) );
}
/* Start function: allocate space for a pair of doubles in which we can
* accumulate the sum and the sum of squares.
*/
static void *
deviate_start( IMAGE *out, void *a, void *b )
{
double *ss2;
if( !(ss2 = IM_ARRAY( NULL, 2, double )) )
return( NULL );
ss2[0] = 0.0;
ss2[1] = 0.0;
return( (void *) ss2 );
}
/* Stop function. Add this little sum to the main sum.
*/
static int
deviate_stop( void *seq, void *a, void *b )
{
double *ss2 = (double *) seq;
double *global_ss2 = (double *) b;
global_ss2[0] += ss2[0];
global_ss2[1] += ss2[1];
im_free( seq );
return( 0 );
}
/* Sum pels in this section.
*/
#define LOOP( TYPE ) { \
TYPE *p = (TYPE *) in; \
\
for( x = 0; x < sz; x++ ) { \
TYPE v = p[x]; \
\
s += v; \
s2 += (double) v * (double) v; \
} \
}
/* Loop over region, accumulating a sum in *tmp.
*/
static int
deviate_scan( void *in, int n, void *seq, void *a, void *b )
{
const IMAGE *im = (IMAGE *) a;
const int sz = n * im->Bands;
double *ss2 = (double *) seq;
int x;
double s, s2;
s = ss2[0];
s2 = ss2[1];
/* Now generate code for all types.
*/
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;
default:
g_assert( 0 );
}
ss2[0] = s;
ss2[1] = s2;
return( 0 );
}
/**
* im_deviate:
* @in: input #IMAGE
* @out: output pixel standard deviation
*
* This operation finds the standard deviation of all pixels in @in. It
* operates on all bands of the input image: use im_stats() if you need
* to calculate an average for each band.
*
* Non-complex images only.
*
* See also: im_stats(), im_bandmean(), im_avg(), im_rank()
*
* Returns: 0 on success, -1 on error
*/
int
im_deviate( IMAGE *in, double *out )
{
double global_ss2[2];
double s, s2;
gint64 vals;
/* Check our args.
*/
if( im_pincheck( in ) ||
im_check_uncoded( "im_deviate", in ) ||
im_check_noncomplex( "im_deviate", in ) )
return( -1 );
/* Loop over input, summing pixels.
*/
global_ss2[0] = 0.0;
global_ss2[1] = 0.0;
if( im__wrapscan( in,
deviate_start, deviate_scan, deviate_stop, in, global_ss2 ) )
return( -1 );
/*
NOTE: NR suggests a two-pass algorithm to minimise roundoff.
But that's too expensive for us :-( so do it the old one-pass
way.
*/
/* Calculate and return deviation. Add a fabs to stop sqrt(<=0).
*/
vals = (gint64) in->Xsize * (gint64) in->Ysize * (gint64) in->Bands;
s = global_ss2[0];
s2 = global_ss2[1];
*out = sqrt( fabs( s2 - (s * s / vals) ) / (vals - 1) );
return( 0 );
}

View File

@ -68,6 +68,74 @@
#include <vips/vips.h>
#include <vips/internal.h>
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( VipsImage *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, gboolean *stop )
{
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 vips_sink(), but the scan function works a line at a time, like
* im_wrap*(). Shared with im_min(), im_deviate() etc.
*/
int
im__wrapscan( VipsImage *in,
VipsStartFn start, im__wrapscan_fn scan, VipsStopFn 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( vips_sink( in,
wrapscan_start, wrapscan_scan, wrapscan_stop,
&wrapscan, NULL ) );
}
/* Track min/max/sum/sum-of-squares for each thread during a scan.
*/
static void *

View File

@ -934,6 +934,12 @@ im_avg( IMAGE *in, double *out )
return( vips_avg( in, out, NULL ) );
}
int
im_deviate( IMAGE *in, double *out )
{
return( vips_deviate( in, out, NULL ) );
}
int im_generate( VipsImage *im,
im_start_fn start, im_generate_fn generate, im_stop_fn stop,
void *a, void *b )

View File

@ -72,6 +72,8 @@ int vips_subtract( VipsImage *in1, VipsImage *in2, VipsImage **out, ... )
__attribute__((sentinel));
int vips_avg( VipsImage *in, double *out, ... )
__attribute__((sentinel));
int vips_deviate( VipsImage *in, double *out, ... )
__attribute__((sentinel));
int vips_min( VipsImage *in, double *out, ... )
__attribute__((sentinel));
int vips_max( VipsImage *in, double *out, ... )
@ -99,7 +101,6 @@ DOUBLEMASK *im_measure_area( VipsImage *im,
int h, int v,
int *sel, int nsel, const char *name );
DOUBLEMASK *im_stats( VipsImage *in );
int im_deviate( VipsImage *in, double *out );
int im_maxpos_avg( VipsImage *im, double *xpos, double *ypos, double *out );
int im_maxpos_vec( VipsImage *im, int *xpos, int *ypos, double *maxima, int n );
int im_minpos_vec( VipsImage *im, int *xpos, int *ypos, double *minima, int n );

View File

@ -527,6 +527,7 @@ int im_minpos( VipsImage *in, int *xpos, int *ypos, double *out );
int im_max( VipsImage *in, double *out );
int im_maxpos( VipsImage *in, int *xpos, int *ypos, double *out );
int im_avg( VipsImage *in, double *out );
int im_deviate( VipsImage *in, double *out );
int im_invert( VipsImage *in, VipsImage *out );
int im_lintra( double a, VipsImage *in, double b, VipsImage *out );
int im_lintra_vec( int n, double *a, VipsImage *in, double *b, VipsImage *out );