sergiotarxz/libvips
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

stuff

Browse Source
This commit is contained in:
John Cupitt 2009-09-07 17:05:19 +00:00
parent 616bb9ab7f
commit 42dbc81827
8 changed files with 471 additions and 489 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
ChangeLog
View File

@ -26,8 +26,9 @@
- im_remainder() has a float result for float types - im_remainder() has a float result for float types
- im_measure() allows sel == NULL, meaning all patches - im_measure() allows sel == NULL, meaning all patches
- added "deprecated" package - added "deprecated" package
- faster, simpler, better im_max() - faster, simpler, better im_max(), im_min, im_avg(), im_deviate()
- im_max() returns true modulus, not square of modulus, for complex images - im_max() returns true modulus, not square of modulus, for complex images
- im_avg() works for complex, returning average modulus
25/3/09 started 7.18.0 25/3/09 started 7.18.0
- revised version numbers - revised version numbers

6
TODO
View File

@ -1,8 +1,4 @@
- new im_max() needs testing - im_stats() needs debugging
im_avg() needs to work for complex: return modulus
use im_max() model to rewrite im_min(), im_stats()
- im_minpos(), im_maxpos() need partialing - im_minpos(), im_maxpos() need partialing

117
libvips/arithmetic/im_avg.c
View File

@ -23,6 +23,9 @@
* - add liboil support * - add liboil support
* 18/8/09 * 18/8/09
* - gtkdoc, minor reformatting * - gtkdoc, minor reformatting
* 7/9/09
* - rewrite for im__wrapiter()
* - add complex case (needed for im_max())
*/ */
/* /*
@ -59,7 +62,6 @@
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <math.h> #include <math.h>
#include <assert.h>
#include <vips/vips.h> #include <vips/vips.h>
#include <vips/internal.h> #include <vips/internal.h>
@ -72,30 +74,32 @@
#include <dmalloc.h> #include <dmalloc.h>
#endif /*WITH_DMALLOC*/ #endif /*WITH_DMALLOC*/
/* Start function: allocate space for a double in which we can accululate the /* Start function: allocate space for a double in which we can accumulate the
* sum. * sum.
*/ */
static void * static void *
start_fn( IMAGE *out, void *a, void *b ) avg_start( IMAGE *out, void *a, void *b )
{ {
double *tmp; double *sum;
if( !(tmp = IM_ARRAY( out, 1, double )) ) if( !(sum = IM_NEW( NULL, double )) )
return( NULL ); return( NULL );
*tmp = 0.0; *sum = 0.0;
return( (void *) tmp ); return( (void *) sum );
} }
/* Stop function. Add this little sum to the main sum. /* Stop function. Add this little sum to the main sum.
*/ */
static int static int
stop_fn( void *seq, void *a, void *b ) avg_stop( void *seq, void *a, void *b )
{ {
double *tmp = (double *) seq; double *sum = (double *) seq;
double *sum = (double *) a; double *global_sum = (double *) b;
*sum += *tmp; *global_sum += *sum;
im_free( seq );
return( 0 ); return( 0 );
} }
@ -103,30 +107,41 @@ stop_fn( void *seq, void *a, void *b )
/* Sum pels in this section. /* Sum pels in this section.
*/ */
#define LOOP( TYPE ) { \ #define LOOP( TYPE ) { \
TYPE *p; \ TYPE *p = (TYPE *) in; \
\
for( y = to; y < bo; y++ ) { \
p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \
\ \
for( x = 0; x < sz; x++ ) \ for( x = 0; x < sz; x++ ) \
sum += *p++; \ m += p[x]; \
}
#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; \
\
m += mod; \
} \ } \
} }
/* Loop over region, accumulating a sum in *tmp. /* Loop over region, accumulating a sum in *tmp.
*/ */
static int static int
scan_fn( REGION *reg, void *seq, void *a, void *b ) avg_scan( void *in, int n, void *seq, void *a, void *b )
{ {
double *tmp = (double *) seq; const IMAGE *im = (IMAGE *) a;
Rect *r = &reg->valid; const int sz = n * im->Bands;
IMAGE *im = reg->im;
int le = r->left; double *sum = (double *) seq;
int to = r->top;
int bo = IM_RECT_BOTTOM(r); int x;
int sz = IM_REGION_N_ELEMENTS( reg ); double m;
double sum = 0.0;
int x, y; m = *sum;
/* Now generate code for all types. /* Now generate code for all types.
*/ */
@ -141,26 +156,27 @@ scan_fn( REGION *reg, void *seq, void *a, void *b )
case IM_BANDFMT_DOUBLE: case IM_BANDFMT_DOUBLE:
#ifdef HAVE_LIBOIL #ifdef HAVE_LIBOIL
for( y = to; y < bo; y++ ) { {
double *p = (double *) IM_REGION_ADDR( reg, le, y ); double *p = (double *) in;
double t; double t;
oil_sum_f64( &t, p, sizeof( double ), sz ); oil_sum_f64( &t, p, sizeof( double ), sz );
sum += t; m += t;
} }
#else /*!HAVE_LIBOIL*/ #else /*!HAVE_LIBOIL*/
LOOP( double ); LOOP( double );
#endif /*HAVE_LIBOIL*/ #endif /*HAVE_LIBOIL*/
break; break;
case IM_BANDFMT_COMPLEX: CLOOP( float ); break;
case IM_BANDFMT_DPCOMPLEX: CLOOP( double ); break;
default: default:
assert( 0 ); g_assert( 0 );
} }
/* Add to sum for this sequence. *sum = m;
*/
*tmp += sum;
return( 0 ); return( 0 );
} }
@ -172,9 +188,7 @@ scan_fn( REGION *reg, void *seq, void *a, void *b )
* *
* This operation finds the average value in an image. It operates on all * This operation finds the average value in an image. It operates on all
* bands of the input image: use im_stats() if you need to calculate an * bands of the input image: use im_stats() if you need to calculate an
* average for each band. * average for each band. For complex images, return the average modulus.
*
* Non-complex images only.
* *
* See also: im_stats(), im_bandmean(), im_deviate(), im_rank() * See also: im_stats(), im_bandmean(), im_deviate(), im_rank()
* *
@ -183,7 +197,7 @@ scan_fn( REGION *reg, void *seq, void *a, void *b )
int int
im_avg( IMAGE *in, double *out ) im_avg( IMAGE *in, double *out )
{ {
double sum = 0.0; double global_sum;
gint64 vals, pels; gint64 vals, pels;
/* Check our args. /* Check our args.
@ -195,14 +209,35 @@ im_avg( IMAGE *in, double *out )
/* Loop over input, summing pixels. /* Loop over input, summing pixels.
*/ */
if( im_iterate( in, start_fn, scan_fn, stop_fn, &sum, NULL ) ) global_sum = 0.0;
if( im__wrapscan( in,
avg_start, avg_scan, avg_stop, in, &global_sum ) )
return( -1 ); return( -1 );
/* Calculate and return average. /* Calculate and return average. For complex, we accumulate re*re +
* im*im, so we need to sqrt.
*/ */
pels = (gint64) in->Xsize * in->Ysize; pels = (gint64) in->Xsize * in->Ysize;
vals = pels * in->Bands; vals = pels * in->Bands;
*out = sum / vals; *out = global_sum / vals;
if( im_iscomplex( in ) )
*out = sqrt( *out );
return( 0 );
}
/* Get the value of pixel (0, 0). Use this to init the min/max value for
* im_max()/im_stats()/etc.
*/
int
im__value( IMAGE *im, double *value )
{
IMAGE *t;
if( !(t = im_open_local( im, "im__value", "p" )) ||
im_extract_area( im, t, 0, 0, 1, 1 ) ||
im_avg( t, value ) )
return( -1 );
return( 0 ); return( 0 );
} }

79
libvips/arithmetic/im_deviate.c
View File

@ -25,6 +25,8 @@
* 2/9/09 * 2/9/09
* - gtk-doc comment * - gtk-doc comment
* - minor reformatting * - minor reformatting
* 4/9/09
* - use im__wrapscan()
*/ */
/* /*
@ -77,28 +79,30 @@
* accumulate the sum and the sum of squares. * accumulate the sum and the sum of squares.
*/ */
static void * static void *
start_fn( IMAGE *out, void *a, void *b ) deviate_start( IMAGE *out, void *a, void *b )
{ {
double *tmp; double *ss2;
if( !(tmp = IM_ARRAY( out, 2, double )) ) if( !(ss2 = IM_ARRAY( NULL, 2, double )) )
return( NULL ); return( NULL );
tmp[0] = 0.0; ss2[0] = 0.0;
tmp[1] = 0.0; ss2[1] = 0.0;
return( (void *) tmp ); return( (void *) ss2 );
} }
/* Stop function. Add this little sum to the main sum. /* Stop function. Add this little sum to the main sum.
*/ */
static int static int
stop_fn( void *seq, void *a, void *b ) deviate_stop( void *seq, void *a, void *b )
{ {
double *tmp = (double *) seq; double *ss2 = (double *) seq;
double *sum = (double *) a; double *global_ss2 = (double *) b;
sum[0] += tmp[0]; global_ss2[0] += ss2[0];
sum[1] += tmp[1]; global_ss2[1] += ss2[1];
im_free( seq );
return( 0 ); return( 0 );
} }
@ -106,10 +110,7 @@ stop_fn( void *seq, void *a, void *b )
/* Sum pels in this section. /* Sum pels in this section.
*/ */
#define LOOP( TYPE ) { \ #define LOOP( TYPE ) { \
TYPE *p; \ TYPE *p = (TYPE *) in; \
\
for( y = to; y < bo; y++ ) { \
p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \
\ \
for( x = 0; x < sz; x++ ) { \ for( x = 0; x < sz; x++ ) { \
TYPE v = p[x]; \ TYPE v = p[x]; \
@ -117,24 +118,23 @@ stop_fn( void *seq, void *a, void *b )
s += v; \ s += v; \
s2 += (double) v * (double) v; \ s2 += (double) v * (double) v; \
} \ } \
} \
} }
/* Loop over region, adding information to the appropriate fields of tmp. /* Loop over region, accumulating a sum in *tmp.
*/ */
static int static int
scan_fn( REGION *reg, void *seq, void *a, void *b ) deviate_scan( void *in, int n, void *seq, void *a, void *b )
{ {
double *tmp = (double *) seq; const IMAGE *im = (IMAGE *) a;
Rect *r = &reg->valid; const int sz = n * im->Bands;
IMAGE *im = reg->im;
int le = r->left; double *ss2 = (double *) seq;
int to = r->top;
int bo = IM_RECT_BOTTOM(r); int x;
int sz = IM_REGION_N_ELEMENTS( reg ); double s, s2;
double s = 0.0;
double s2 = 0.0; s = ss2[0];
int x, y; s2 = ss2[1];
/* Now generate code for all types. /* Now generate code for all types.
*/ */
@ -149,8 +149,8 @@ scan_fn( REGION *reg, void *seq, void *a, void *b )
case IM_BANDFMT_DOUBLE: case IM_BANDFMT_DOUBLE:
#ifdef HAVE_LIBOIL #ifdef HAVE_LIBOIL
for( y = to; y < bo; y++ ) { {
double *p = (double *) IM_REGION_ADDR( reg, le, y ); double *p = (double *) in;
double t; double t;
double t2; double t2;
@ -159,7 +159,7 @@ scan_fn( REGION *reg, void *seq, void *a, void *b )
s += t; s += t;
s2 += t2; s2 += t2;
} }
#else /*!HAVE_LIBOIL*/ #else /*!HAVE_LIBOIL*/
LOOP( double ); LOOP( double );
#endif /*HAVE_LIBOIL*/ #endif /*HAVE_LIBOIL*/
@ -169,10 +169,8 @@ scan_fn( REGION *reg, void *seq, void *a, void *b )
g_assert( 0 ); g_assert( 0 );
} }
/* Add to sum for this sequence. ss2[0] += s;
*/ ss2[1] += s2;
tmp[0] += s;
tmp[1] += s2;
return( 0 ); return( 0 );
} }
@ -195,7 +193,7 @@ scan_fn( REGION *reg, void *seq, void *a, void *b )
int int
im_deviate( IMAGE *in, double *out ) im_deviate( IMAGE *in, double *out )
{ {
double sum[2] = { 0.0, 0.0 }; double global_ss2[2];
gint64 N; gint64 N;
/* Check our args. /* Check our args.
@ -207,8 +205,10 @@ im_deviate( IMAGE *in, double *out )
/* Loop over input, summing pixels. /* Loop over input, summing pixels.
*/ */
if( im_iterate( in, start_fn, scan_fn, stop_fn, &sum, NULL ) ) global_ss2[0] = 0.0;
return( -1 ); global_ss2[1] = 0.0;
if( im__wrapscan( in,
deviate_start, deviate_scan, deviate_stop, in, &global_ss2 ) )
/* /*
@ -221,7 +221,8 @@ im_deviate( IMAGE *in, double *out )
/* Calculate and return deviation. Add a fabs to stop sqrt(<=0). /* Calculate and return deviation. Add a fabs to stop sqrt(<=0).
*/ */
N = (gint64) in->Xsize * in->Ysize * in->Bands; N = (gint64) in->Xsize * in->Ysize * in->Bands;
*out = sqrt( fabs( sum[1] - (sum[0] * sum[0] / N) ) / (N - 1) ); *out = sqrt( fabs( global_ss2[1] -
(global_ss2[0] * global_ss2[0] / N) ) / (N - 1) );
return( 0 ); return( 0 );
} }

154
libvips/arithmetic/im_max.c
View File

@ -63,35 +63,106 @@
#include <math.h> #include <math.h>
#include <vips/vips.h> #include <vips/vips.h>
#include <vips/internal.h>
#ifdef WITH_DMALLOC #ifdef WITH_DMALLOC
#include <dmalloc.h> #include <dmalloc.h>
#endif /*WITH_DMALLOC*/ #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. /* New sequence value.
*/ */
static void * static void *
max_start( IMAGE *in, void *a, void *b ) max_start( IMAGE *in, void *a, void *b )
{ {
double *value = (double *) a; double *global_max = (double *) b;
double *seq = IM_NEW( NULL, double ); double *max;
*seq = *value; if( !(max = IM_NEW( NULL, double )) )
return( NULL );
*max = *global_max;
return( (void *) seq ); return( (void *) max );
} }
/* Merge the sequence value back into the per-call state. /* Merge the sequence value back into the per-call state.
*/ */
static int static int
max_stop( void *vseq, void *a, void *b ) max_stop( void *seq, void *a, void *b )
{ {
double *seq = (double *) vseq; double *max = (double *) seq;
double *value = (double *) a; double *global_max = (double *) b;
/* Merge. /* Merge.
*/ */
*value = IM_MAX( *value, *seq ); *global_max = IM_MAX( *global_max, *max );
im_free( seq ); im_free( seq );
@ -99,23 +170,20 @@ max_stop( void *vseq, void *a, void *b )
} }
#define LOOP( TYPE ) { \ #define LOOP( TYPE ) { \
for( y = to; y < bo; y++ ) { \ TYPE *p = (TYPE *) in; \
TYPE *p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \
\ \
for( x = 0; x < nel; x++ ) { \ for( x = 0; x < sz; x++ ) { \
double v = p[x]; \ double v = p[x]; \
\ \
if( v > m ) \ if( v > m ) \
m = v; \ m = v; \
} \ } \
} \
} }
#define CLOOP( TYPE ) { \ #define CLOOP( TYPE ) { \
for( y = to; y < bo; y++ ) { \ TYPE *p = (TYPE *) in; \
TYPE *p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \
\ \
for( x = 0; x < nel; x++ ) { \ for( x = 0; x < sz; x++ ) { \
double mod, re, im; \ double mod, re, im; \
\ \
re = p[0]; \ re = p[0]; \
@ -126,26 +194,22 @@ max_stop( void *vseq, void *a, void *b )
if( mod > m ) \ if( mod > m ) \
m = mod; \ m = mod; \
} \ } \
} \
} }
/* Loop over region, adding to seq. /* Loop over region, adding to seq.
*/ */
static int static int
max_scan( REGION *reg, void *vseq, void *a, void *b ) max_scan( void *in, int n, void *seq, void *a, void *b )
{ {
double *seq = (double *) vseq; const IMAGE *im = (IMAGE *) a;
Rect *r = &reg->valid; const int sz = n * im->Bands;
IMAGE *im = reg->im;
int le = r->left;
int to = r->top;
int bo = IM_RECT_BOTTOM(r);
int nel = IM_REGION_N_ELEMENTS( reg );
int x, y; double *max = (double *) seq;
int x;
double m; double m;
m = *seq; m = *max;
switch( im->BandFmt ) { switch( im->BandFmt ) {
case IM_BANDFMT_UCHAR: LOOP( unsigned char ); break; case IM_BANDFMT_UCHAR: LOOP( unsigned char ); break;
@ -163,30 +227,7 @@ max_scan( REGION *reg, void *vseq, void *a, void *b )
g_assert( 0 ); g_assert( 0 );
} }
*seq = m; *max = m;
#ifdef DEBUG
printf( "im_max: left = %d, top = %d, width = %d, height = %d\n",
r->left, r->top, r->width, r->height );
printf( " (max = %g)\n", *seq );
#endif /*DEBUG*/
return( 0 );
}
/* Get the value of pixel (0, 0). Use this to init the min/max value for
* threads. Shared with im_min(), im_stats() etc. This will return mod for
* complex.
*/
int
im__value( IMAGE *im, double *value )
{
IMAGE *t;
if( !(t = im_open_local( im, "im__value", "p" )) ||
im_extract_area( im, t, 0, 0, 1, 1 ) ||
im_avg( t, value ) )
return( -1 );
return( 0 ); return( 0 );
} }
@ -208,28 +249,29 @@ im__value( IMAGE *im, double *value )
int int
im_max( IMAGE *in, double *out ) im_max( IMAGE *in, double *out )
{ {
double value; double global_max;
if( im_pincheck( in ) || if( im_pincheck( in ) ||
im_check_uncoded( "im_max", in ) ) im_check_uncoded( "im_max", in ) )
return( -1 ); return( -1 );
if( im__value( in, &value ) ) if( im__value( in, &global_max ) )
return( -1 ); return( -1 );
/* We use square mod for scanning, for speed. /* We use square mod for scanning, for speed.
*/ */
if( im_iscomplex( in ) ) if( im_iscomplex( in ) )
value *= value; global_max *= global_max;
if( im_iterate( in, max_start, max_scan, max_stop, &value, NULL ) ) if( im__wrapscan( in, max_start, max_scan, max_stop,
in, &global_max ) )
return( -1 ); return( -1 );
/* Back to modulus. /* Back to modulus.
*/ */
if( im_iscomplex( in ) ) if( im_iscomplex( in ) )
value = sqrt( value ); global_max = sqrt( global_max );
*out = value; *out = global_max;
return( 0 ); return( 0 );
} }

188
libvips/arithmetic/im_min.c
View File

@ -16,6 +16,8 @@
* - partialed, based in im_max() * - partialed, based in im_max()
* 3/4/02 JC * 3/4/02 JC
* - random wrong result for >1 thread :-( (thanks Joe) * - random wrong result for >1 thread :-( (thanks Joe)
* 4/9/09
* - rewrite from im_max()
*/ */
/* /*
@ -56,160 +58,105 @@
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <math.h> #include <math.h>
#include <assert.h>
#include <vips/vips.h> #include <vips/vips.h>
#include <vips/internal.h>
#ifdef WITH_DMALLOC #ifdef WITH_DMALLOC
#include <dmalloc.h> #include <dmalloc.h>
#endif /*WITH_DMALLOC*/ #endif /*WITH_DMALLOC*/
/* Per-call state.
*/
typedef struct _MinInfo {
/* Parameters.
*/
IMAGE *in;
double *out;
/* Global min so far.
*/
double value;
int valid; /* zero means value is unset */
} MinInfo;
/* Per thread state.
*/
typedef struct _Seq {
MinInfo *inf;
double value;
int valid; /* zero means value is unset */
} Seq;
/* New sequence value. /* New sequence value.
*/ */
static void * static void *
start_fn( IMAGE *im, void *a, void *b ) min_start( IMAGE *in, void *a, void *b )
{ {
MinInfo *inf = (MinInfo *) a; double *global_min = (double *) b;
Seq *seq = IM_NEW( NULL, Seq ); double *min;
seq->inf = inf; if( !(min = IM_NEW( NULL, double )) )
seq->valid = 0; return( NULL );
*min = *global_min;
return( seq ); return( (void *) min );
} }
/* Merge the sequence value back into the per-call state. /* Merge the sequence value back into the per-call state.
*/ */
static int static int
stop_fn( void *vseq, void *a, void *b ) min_stop( void *seq, void *a, void *b )
{ {
Seq *seq = (Seq *) vseq; double *min = (double *) seq;
MinInfo *inf = (MinInfo *) a; double *global_min = (double *) b;
if( seq->valid ) {
if( !inf->valid )
/* Just copy.
*/
inf->value = seq->value;
else
/* Merge. /* Merge.
*/ */
inf->value = IM_MIN( inf->value, seq->value ); *global_min = IM_MIN( *global_min, *min );
inf->valid = 1;
}
im_free( seq ); im_free( seq );
return( 0 ); return( 0 );
} }
/* Loop over region, adding to seq. #define LOOP( TYPE ) { \
*/ TYPE *p = (TYPE *) in; \
static int
scan_fn( REGION *reg, void *vseq, void *a, void *b )
{
Seq *seq = (Seq *) vseq;
Rect *r = &reg->valid;
IMAGE *im = reg->im;
int le = r->left;
int to = r->top;
int bo = IM_RECT_BOTTOM(r);
int nel = IM_REGION_N_ELEMENTS( reg );
int x, y;
double m;
#ifdef DEBUG
printf( "im_min: left = %d, top = %d, width = %d, height = %d\n",
r->left, r->top, r->width, r->height );
#endif /*DEBUG*/
#define loop(TYPE) { \
m = *((TYPE *) IM_REGION_ADDR( reg, le, to )); \
\ \
for( y = to; y < bo; y++ ) { \ for( x = 0; x < sz; x++ ) { \
TYPE *p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \
\
for( x = 0; x < nel; x++ ) { \
double v = p[x]; \ double v = p[x]; \
\ \
if( v < m ) \ if( v < m ) \
m = v; \ m = v; \
} \ } \
} \
} }
#define complex_loop(TYPE) { \ #define CLOOP( TYPE ) { \
TYPE *p = (TYPE *) IM_REGION_ADDR( reg, le, to ); \ TYPE *p = (TYPE *) in; \
double real = p[0]; \
double imag = p[1]; \
\ \
m = real * real + imag * imag; \ for( x = 0; x < sz; x++ ) { \
double mod, re, im; \
\ \
for( y = to; y < bo; y++ ) { \ re = p[0]; \
TYPE *p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \ im = p[1]; \
\ p += 2; \
for( x = 0; x < nel * 2; x += 2 ) { \ mod = re * re + im * im; \
double mod; \
\
real = p[x]; \
imag = p[x + 1]; \
mod = real * real + imag * imag; \
\ \
if( mod < m ) \ if( mod < m ) \
m = mod; \ 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 ) { switch( im->BandFmt ) {
case IM_BANDFMT_UCHAR: loop( unsigned char ); break; case IM_BANDFMT_UCHAR: LOOP( unsigned char ); break;
case IM_BANDFMT_CHAR: loop( signed char ); break; case IM_BANDFMT_CHAR: LOOP( signed char ); break;
case IM_BANDFMT_USHORT: loop( unsigned short ); break; case IM_BANDFMT_USHORT: LOOP( unsigned short ); break;
case IM_BANDFMT_SHORT: loop( signed short ); break; case IM_BANDFMT_SHORT: LOOP( signed short ); break;
case IM_BANDFMT_UINT: loop( unsigned int ); break; case IM_BANDFMT_UINT: LOOP( unsigned int ); break;
case IM_BANDFMT_INT: loop( signed int ); break; case IM_BANDFMT_INT: LOOP( signed int ); break;
case IM_BANDFMT_FLOAT: loop( float ); break; case IM_BANDFMT_FLOAT: LOOP( float ); break;
case IM_BANDFMT_DOUBLE: loop( double ); break; case IM_BANDFMT_DOUBLE: LOOP( double ); break;
case IM_BANDFMT_COMPLEX: complex_loop( float ); break; case IM_BANDFMT_COMPLEX: CLOOP( float ); break;
case IM_BANDFMT_DPCOMPLEX: complex_loop( double ); break; case IM_BANDFMT_DPCOMPLEX: CLOOP( double ); break;
default: default:
assert( 0 ); g_assert( 0 );
} }
if( seq->valid ) { *min = m;
seq->value = IM_MIN( seq->value, m );
}
else {
seq->value = m;
seq->valid = 1;
}
return( 0 ); return( 0 );
} }
@ -219,10 +166,10 @@ scan_fn( REGION *reg, void *vseq, void *a, void *b )
* @in: input #IMAGE * @in: input #IMAGE
* @out: output double * @out: output double
* *
* Finds the the minimum value of image @in and returns it at the * Finds the the minimum value of image #in and returns it at the
* location pointed by @out. If input is complex, the min square amplitude * location pointed by out. If input is complex, the min modulus
* (re*re+im*im) is returned. im_min() finds the maximum of all bands: if you * is returned. im_min() finds the minimum of all bands: if you
* want to find the maximum of each band separately, use im_stats(). * want to find the minimum of each band separately, use im_stats().
* *
* See also: im_minpos(), im_max(), im_stats(). * See also: im_minpos(), im_max(), im_stats().
* *
@ -231,20 +178,29 @@ scan_fn( REGION *reg, void *vseq, void *a, void *b )
int int
im_min( IMAGE *in, double *out ) im_min( IMAGE *in, double *out )
{ {
MinInfo inf; double global_min;
inf.in = in;
inf.out = out;
inf.valid = 0;
if( im_pincheck( in ) || if( im_pincheck( in ) ||
im_check_uncoded( "im_min", in ) ) im_check_uncoded( "im_min", in ) )
return( -1 ); return( -1 );
if( im_iterate( in, start_fn, scan_fn, stop_fn, &inf, NULL ) ) 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 ); return( -1 );
*out = inf.value; /* Back to modulus.
*/
if( im_iscomplex( in ) )
global_min = sqrt( global_min );
*out = global_min;
return( 0 ); return( 0 );
} }

281
libvips/arithmetic/im_stats.c
View File

@ -17,12 +17,12 @@
* 13/1/05 * 13/1/05
* - use 64 bit arithmetic * - use 64 bit arithmetic
* 1/9/09
1/9/09 * - argh nope min/max was broken again for >1CPU in short pipelines on
- argh nope min/max was broken again for >1CPU in short pipelines on * some architectures
some architectures * 7/9/09
* - rework based on new im__wrapscan() / im_max() ideas for a proper fix
*/ */
/* /*
@ -62,183 +62,119 @@
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <math.h> #include <math.h>
#include <assert.h>
#include <vips/vips.h> #include <vips/vips.h>
#include <vips/internal.h>
#ifdef WITH_DMALLOC #ifdef WITH_DMALLOC
#include <dmalloc.h> #include <dmalloc.h>
#endif /*WITH_DMALLOC*/ #endif /*WITH_DMALLOC*/
/* Make and initialise a DOUBLEMASK suitable for grabbing statistics. /* Track min/max/sum/sum-of-squares for each thread during a scan.
*/ */
static void * static void *
make_mask( IMAGE *im, void *a, void *b ) stats_start( IMAGE *im, void *a, void *b )
{ {
DOUBLEMASK *out; double *global_stats = (double *) b;
/* Make temp output. double *stats;
*/ int i;
if( !(out = im_create_dmask( "stats", 6, im->Bands + 1 )) )
if( !(stats = IM_ARRAY( NULL, 4 * im->Bands, double )) )
return( NULL ); return( NULL );
/* Set offset to magic value: 42 indicates we have not yet initialised for( i = 0; i < 4 * im->Bands; i++ )
* max and min for this mask. stats[i] = global_stats[i];
*/
out->offset = 42;
#ifdef DEBUG return( stats );
printf( "make_mask: created %p\n", out );
#endif /*DEBUG*/
return( out );
} }
/* Merge a temp DOUBLEMASK into the real DOUBLEMASK. Row 0 is unused, row 1 /* Merge a thread's output back into the global stats.
* has the stats for band 1. These are: (minimum, maximum, sum, sum^2). If the
* offset of out is 42, then it has not been inited yet and we just copy.
*/ */
static int static int
merge_mask( void *seq, void *a, void *b ) stats_stop( void *seq, void *a, void *b )
{ {
DOUBLEMASK *tmp = (DOUBLEMASK *) seq; const IMAGE *im = (IMAGE *) a;
DOUBLEMASK *out = (DOUBLEMASK *) a; double *global_stats = (double *) b;
double *rowi, *rowo; double *stats = (double *) seq;
int z;
#ifdef DEBUG int i;
printf( "merge_mask: tmp = %p, out = %p\n", tmp, out );
im_print_dmask( tmp );
im_print_dmask( out );
#endif /*DEBUG*/
/* Merge, or just copy? Also, tmp might be uninited, so allow for that for( i = 0; i < 4 * im->Bands; i += 4 ) {
* too. global_stats[0] = IM_MIN( global_stats[0], stats[0] );
*/ global_stats[1] = IM_MAX( global_stats[1], stats[1] );
if( out->offset == 42 && tmp->offset != 42 ) { global_stats[2] += stats[2];
#ifdef DEBUG global_stats[3] += stats[3];
printf( "merge_mask: copying\n", tmp );
#endif /*DEBUG*/
/* Copy info from tmp. global_stats += 4;
*/ stats += 4;
for( z = 1; z < tmp->ysize; z++ ) {
rowi = tmp->coeff + z * 6;
rowo = out->coeff + z * 6;
rowo[0] = rowi[0];
rowo[1] = rowi[1];
rowo[2] = rowi[2];
rowo[3] = rowi[3];
} }
out->offset = 0; im_free( seq );
}
else if( out->offset != 42 && tmp->offset != 42 ) {
#ifdef DEBUG
printf( "merge_mask: merging\n" );
#endif /*DEBUG*/
/* Add info from tmp.
*/
for( z = 1; z < tmp->ysize; z++ ) {
rowi = tmp->coeff + z * 6;
rowo = out->coeff + z * 6;
rowo[0] = IM_MIN( rowi[0], rowo[0] );
rowo[1] = IM_MAX( rowi[1], rowo[1] );
rowo[2] += rowi[2];
rowo[3] += rowi[3];
}
}
/* Can now free tmp.
*/
im_free_dmask( tmp );
return( 0 ); return( 0 );
} }
/* Loop over region, adding information to the appropriate fields of tmp. /* We scan lines bands times to avoid repeating band loops.
* We set max, min, sum, sum of squares. Our caller fills in the rest.
*/
static int
scan_fn( REGION *reg, void *seq, void *a, void *b )
{
DOUBLEMASK *tmp = (DOUBLEMASK *) seq;
Rect *r = &reg->valid;
IMAGE *im = reg->im;
int bands = im->Bands;
int le = r->left;
int ri = IM_RECT_RIGHT(r);
int to = r->top;
int bo = IM_RECT_BOTTOM(r);
int x, y, z;
/* What type? First define the loop we want to perform for all types.
* We scan lines bands times to avoid repeating band loops.
* Use temp variables of same type for min/max for faster comparisons. * Use temp variables of same type for min/max for faster comparisons.
* Use double to sum bands.
*/ */
#define non_complex_loop(TYPE) \ #define LOOP( TYPE ) { \
{ TYPE *p, *q; \ for( z = 0; z < im->Bands; z++ ) { \
TYPE value, small, big; \ TYPE *q = (TYPE *) in + z; \
double *row; \ double *row = stats + z * 4; \
TYPE small, big; \
double sum, sum2; \
\ \
/* Have min and max been initialised? \
*/ \
if( tmp->offset == 42 ) { \
/* Init min and max for each band. \
*/ \
p = (TYPE *) IM_REGION_ADDR( reg, le, to ); \
for( z = 0; z < bands; z++ ) { \
row = tmp->coeff + (z + 1) * 6; \
row[0] = p[z]; \
row[1] = p[z]; \
} \
tmp->offset = 0; \
} \
\
for( y = to; y < bo; y++ ) { \
p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \
\
for( z = 0; z < bands; z++ ) { \
q = p + z; \
row = tmp->coeff + (z + 1)*6; \
small = row[0]; \ small = row[0]; \
big = row[1]; \ big = row[1]; \
sum = row[2]; \
sum2 = row[3]; \
\ \
for( x = le; x < ri; x++ ) { \ for( x = 0; x < n; x++ ) { \
value = *q; \ TYPE value = *q; \
q += bands; \ \
row[2] += value;\ sum += value;\
row[3] += (double)value*(double)value;\ sum2 += (double) value * (double) value;\
if( value > big ) \ if( value > big ) \
big = value; \ big = value; \
else if( value < small ) \ else if( value < small ) \
small = value;\ small = value;\
\
q += im->Bands; \
}\ }\
\ \
row[0] = small; \ row[0] = small; \
row[1] = big; \ row[1] = big; \
row[2] = sum; \
row[3] = sum2; \
}\ }\
}\ }
}
/* Loop over region, adding to seq.
*/
static int
stats_scan( void *in, int n, void *seq, void *a, void *b )
{
const IMAGE *im = (IMAGE *) a;
double *stats = (double *) seq;
int x, z;
/* Now generate code for all types. /* Now generate code for all types.
*/ */
switch( im->BandFmt ) { switch( im->BandFmt ) {
case IM_BANDFMT_UCHAR: non_complex_loop(unsigned char); break; case IM_BANDFMT_UCHAR: LOOP( unsigned char ); break;
case IM_BANDFMT_CHAR: non_complex_loop(signed char); break; case IM_BANDFMT_CHAR: LOOP( signed char ); break;
case IM_BANDFMT_USHORT: non_complex_loop(unsigned short); break; case IM_BANDFMT_USHORT: LOOP( unsigned short ); break;
case IM_BANDFMT_SHORT: non_complex_loop(signed short); break; case IM_BANDFMT_SHORT: LOOP( signed short ); break;
case IM_BANDFMT_UINT: non_complex_loop(unsigned int); break; case IM_BANDFMT_UINT: LOOP( unsigned int ); break;
case IM_BANDFMT_INT: non_complex_loop(signed int); break; case IM_BANDFMT_INT: LOOP( signed int ); break;
case IM_BANDFMT_DOUBLE: non_complex_loop(double); break; case IM_BANDFMT_DOUBLE: LOOP( double ); break;
case IM_BANDFMT_FLOAT: non_complex_loop(float); break; case IM_BANDFMT_FLOAT: LOOP( float ); break;
default: default:
assert( 0 ); g_assert( 0 );
} }
return( 0 ); return( 0 );
@ -251,58 +187,69 @@ scan_fn( REGION *reg, void *seq, void *a, void *b )
* the figures for all bands together. * the figures for all bands together.
*/ */
DOUBLEMASK * DOUBLEMASK *
im_stats( IMAGE *in ) im_stats( IMAGE *im )
{ {
DOUBLEMASK *out; DOUBLEMASK *out;
double *row, *base; double *row;
gint64 pels, vals, z; gint64 pels, vals, z;
double *global_stats;
int i, j;
double value;
/* Check our args. if( im_pincheck( im ) ||
*/ im_check_noncomplex( "im_stats", im ) ||
if( im_pincheck( in ) ) im_check_uncoded( "im_stats", im ) )
return( NULL ); return( NULL );
if( im_iscomplex( in ) ) {
im_error( "im_stats", "%s", _( "bad input type" ) );
return( NULL );
}
if( in->Coding != IM_CODING_NONE ) {
im_error( "im_stats", "%s", _( "not uncoded" ) );
return( NULL );
}
/* Make output. if( !(global_stats = IM_ARRAY( im, 4 * im->Bands, double )) )
*/
pels = (gint64) in->Xsize * in->Ysize;
vals = pels * in->Bands;
if( !(out = make_mask( in, NULL, NULL )) )
return( NULL ); return( NULL );
if( im__value( im, &value ) )
return( NULL );
for( i = 0; i < 4 * im->Bands; i += 4 ) {
global_stats[i + 0] = value;
global_stats[i + 1] = value;
global_stats[i + 2] = 0.0;
global_stats[i + 3] = 0.0;
}
/* Loop over input, calculating min, max, sum, sum^2 for each band /* Loop over input, calculating min, max, sum, sum^2 for each band
* separately. * separately.
*/ */
if( im_iterate( in, make_mask, scan_fn, merge_mask, out, NULL ) ) { if( im__wrapscan( im, stats_start, stats_scan, stats_stop,
im_free_dmask( out ); im, &global_stats ) )
return( NULL ); return( NULL );
}
/* Calculate mean, deviation, plus overall stats. /* Calculate mean, deviation, plus overall stats.
*/ */
base = out->coeff; if( !(out = im_create_dmask( "stats", 6, im->Bands + 1 )) )
base[0] = base[6]; /* Init global max/min */ return( NULL );
base[1] = base[7];
for( z = 0; z < in->Bands; z++ ) { /* Init global max/min/sum/sum2.
row = base + (z + 1) * 6; */
base[0] = IM_MIN( base[0], row[0] ); out->coeff[0] = value;
base[1] = IM_MAX( base[1], row[1] ); out->coeff[1] = value;
base[2] += row[2]; out->coeff[2] = 0.0;
base[3] += row[3]; out->coeff[3] = 0.0;
pels = (gint64) im->Xsize * im->Ysize;
vals = pels * im->Bands;
for( i = 0; i < im->Bands; i++ ) {
row = out->coeff + (z + 1) * 6;
for( j = 0; j < 4; j++ )
row[j] = global_stats[i * 4 + j];
out->coeff[0] = IM_MIN( out->coeff[0], row[0] );
out->coeff[1] = IM_MAX( out->coeff[1], row[1] );
out->coeff[2] += row[2];
out->coeff[3] += row[3];
row[4] = row[2] / pels; row[4] = row[2] / pels;
row[5] = sqrt( fabs( row[3] - (row[2] * row[2] / pels) ) / row[5] = sqrt( fabs( row[3] - (row[2] * row[2] / pels) ) /
(pels - 1) ); (pels - 1) );
} }
base[4] = base[2] / vals; out->coeff[4] = out->coeff[2] / vals;
base[5] = sqrt( fabs( base[3] - (base[2] * base[2] / vals) ) / out->coeff[5] = sqrt( fabs( out->coeff[3] -
(vals - 1) ); (out->coeff[2] * out->coeff[2] / vals) ) / (vals - 1) );
#ifdef DEBUG #ifdef DEBUG
printf( "im_stats:\n" ); printf( "im_stats:\n" );

4
libvips/include/vips/internal.h
View File

@ -136,6 +136,10 @@ int im__cast_and_call( IMAGE *in1, IMAGE *in2, IMAGE *out,
VipsBandFmt im__format_common( IMAGE *in1, IMAGE *in2 ); VipsBandFmt im__format_common( IMAGE *in1, IMAGE *in2 );
int im__math( const char *name, IMAGE *in, IMAGE *out, im_wrapone_fn gen ); int im__math( const char *name, IMAGE *in, IMAGE *out, im_wrapone_fn gen );
int im__value( IMAGE *im, double *value ); int im__value( IMAGE *im, double *value );
typedef int (*im__wrapscan_fn)( void *p, int n, void *seq, void *a, void *b );
int im__wrapscan( IMAGE *in,
im_start_fn start, im__wrapscan_fn scan, im_stop_fn stop,
void *a, void *b );
int im__test_kill( IMAGE *im ); int im__test_kill( IMAGE *im );
void *im__mmap( int fd, int writeable, size_t length, gint64 offset ); void *im__mmap( int fd, int writeable, size_t length, gint64 offset );