From f0d47605602003d2cfb5ed998c6c9d17626828f8 Mon Sep 17 00:00:00 2001 From: John Cupitt Date: Fri, 8 Nov 2013 14:09:42 +0000 Subject: [PATCH] redo correlation funcs as classes im_fastcor() im_spcor() im_gradcor() though gradcor we just deprecate, it's complex and hardly used ... rework later is anyone complains --- ChangeLog | 3 +- TODO | 12 + libvips/arithmetic/add.c | 2 +- libvips/arithmetic/boolean.c | 2 +- libvips/arithmetic/divide.c | 2 +- libvips/arithmetic/linear.c | 2 +- libvips/arithmetic/math2.c | 2 +- libvips/arithmetic/multiply.c | 6 +- libvips/arithmetic/relational.c | 2 +- libvips/arithmetic/remainder.c | 2 +- libvips/arithmetic/subtract.c | 2 +- libvips/conversion/embed.c | 2 +- libvips/convolution/Makefile.am | 10 +- libvips/convolution/convolution.c | 4 + libvips/convolution/correlation.c | 165 ++++++++ libvips/convolution/correlation.h | 96 +++++ libvips/convolution/fastcor.c | 267 +++++++++++++ libvips/convolution/im_fastcor.c | 206 ---------- libvips/convolution/im_sharpen.c | 12 +- libvips/convolution/im_spcor.c | 348 ----------------- libvips/convolution/spcor.c | 367 ++++++++++++++++++ libvips/deprecated/Makefile.am | 2 + .../convol_dispatch.c | 0 .../{convolution => deprecated}/im_gradcor.c | 0 libvips/deprecated/vips7compat.c | 46 +++ libvips/foreign/radiance.c | 1 - libvips/include/vips/convolution.h | 7 - libvips/include/vips/vips7compat.h | 7 + 28 files changed, 991 insertions(+), 586 deletions(-) create mode 100644 libvips/convolution/correlation.c create mode 100644 libvips/convolution/correlation.h create mode 100644 libvips/convolution/fastcor.c delete mode 100644 libvips/convolution/im_fastcor.c delete mode 100644 libvips/convolution/im_spcor.c create mode 100644 libvips/convolution/spcor.c rename libvips/{convolution => deprecated}/convol_dispatch.c (100%) rename libvips/{convolution => deprecated}/im_gradcor.c (100%) diff --git a/ChangeLog b/ChangeLog index cba7026e..0a9bb333 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,7 +1,8 @@ 19/10/13 started 7.37.0 - redone im_rotate_*mask45(), im_gauss_*mask*(), im_log_*mask(), im_dilate(), im_erode(), im_rank_image(), im_compass(), im_linedet(), im_gradient(), - im_convsep(), im_convsep_f() as classes + im_convsep(), im_convsep_f(), im_fastcor(), im_spcor() as classes +- im_gradcor() deprecated - vips_init() now does some ABI compat checking, though this change requires an ABI break - add "interlace" option to vips_jpegsave() diff --git a/TODO b/TODO index cf305f7d..c3ae0d90 100644 --- a/TODO +++ b/TODO @@ -1,3 +1,15 @@ +- look again at gcc auto-vectorisation, what would we need to do to use this? + +- how about + + top -d 0.1 | grep vips > log + + to watch RES at 0.1s intervals? + +- add "crop" as a synonym for "extract_area" + +- add vips_gamma(), see rad etc. + - drag display control slider, broken - do conv and morph quickly as simple wrappers over the vips7 operations diff --git a/libvips/arithmetic/add.c b/libvips/arithmetic/add.c index e2e2cd56..30a13ea4 100644 --- a/libvips/arithmetic/add.c +++ b/libvips/arithmetic/add.c @@ -249,7 +249,7 @@ vips_add_init( VipsAdd *add ) * one-band image by joining n copies of the one-band image together, and then * the two n-band images are operated upon. * - * The two input images are cast up to the smallest common type (see table + * The two input images are cast up to the smallest common format (see table * Smallest common format in * arithmetic), then the * following table is used to determine the output type: diff --git a/libvips/arithmetic/boolean.c b/libvips/arithmetic/boolean.c index 1160ae68..9c2434cf 100644 --- a/libvips/arithmetic/boolean.c +++ b/libvips/arithmetic/boolean.c @@ -258,7 +258,7 @@ vips_booleanv( VipsImage *left, VipsImage *right, VipsImage **out, * one-band image by joining n copies of the one-band image together, and then * the two n-band images are operated upon. * - * The two input images are cast up to the smallest common type (see table + * The two input images are cast up to the smallest common format (see table * Smallest common format in * arithmetic). * diff --git a/libvips/arithmetic/divide.c b/libvips/arithmetic/divide.c index c92aa9b3..a82c5718 100644 --- a/libvips/arithmetic/divide.c +++ b/libvips/arithmetic/divide.c @@ -256,7 +256,7 @@ vips_divide_init( VipsDivide *divide ) * one-band image by joining n copies of the one-band image together, and then * the two n-band images are operated upon. * - * The two input images are cast up to the smallest common type (see table + * The two input images are cast up to the smallest common format (see table * Smallest common format in * arithmetic), then the * following table is used to determine the output type: diff --git a/libvips/arithmetic/linear.c b/libvips/arithmetic/linear.c index 155cd600..98408b51 100644 --- a/libvips/arithmetic/linear.c +++ b/libvips/arithmetic/linear.c @@ -180,7 +180,7 @@ vips_linear_build( VipsObject *object ) for( x = 0; x < width; x++ ) \ for( k = 0; k < nb; k++ ) { \ q[0] = a[k] * p[0] + b[k]; \ - q[1] = a[k] * p[1]; \ + q[1] = p[1]; \ q += 2; \ p += 2; \ } \ diff --git a/libvips/arithmetic/math2.c b/libvips/arithmetic/math2.c index 098d51b5..16d57af3 100644 --- a/libvips/arithmetic/math2.c +++ b/libvips/arithmetic/math2.c @@ -248,7 +248,7 @@ vips_math2v( VipsImage *left, VipsImage *right, VipsImage **out, * one-band image by joining n copies of the one-band image together, and then * the two n-band images are operated upon. * - * The two input images are cast up to the smallest common type (see table + * The two input images are cast up to the smallest common format (see table * Smallest common format in * arithmetic), and that format is the * result type. diff --git a/libvips/arithmetic/multiply.c b/libvips/arithmetic/multiply.c index 2a49cff4..5e2fd8f3 100644 --- a/libvips/arithmetic/multiply.c +++ b/libvips/arithmetic/multiply.c @@ -158,7 +158,7 @@ vips_multiply_buffer( VipsArithmetic *arithmetic, /* Type promotion for multiplication. Sign and value preserving. Make sure * these match the case statement in multiply_buffer() above. */ -static int vips_bandfmt_multiply[10] = { +static int vips_multiply_format_table[10] = { /* UC C US S UI I F X D DX */ US, S, UI, I, UI, I, F, X, D, DX }; @@ -172,7 +172,7 @@ vips_multiply_class_init( VipsMultiplyClass *class ) object_class->nickname = "multiply"; object_class->description = _( "multiply two images" ); - vips_arithmetic_set_format_table( aclass, vips_bandfmt_multiply ); + vips_arithmetic_set_format_table( aclass, vips_multiply_format_table ); aclass->process_line = vips_multiply_buffer; } @@ -198,7 +198,7 @@ vips_multiply_init( VipsMultiply *multiply ) * one-band image by joining n copies of the one-band image together, and then * the two n-band images are operated upon. * - * The two input images are cast up to the smallest common type (see table + * The two input images are cast up to the smallest common format (see table * Smallest common format in * arithmetic), then the * following table is used to determine the output type: diff --git a/libvips/arithmetic/relational.c b/libvips/arithmetic/relational.c index f9f3d395..a6064f2a 100644 --- a/libvips/arithmetic/relational.c +++ b/libvips/arithmetic/relational.c @@ -268,7 +268,7 @@ vips_relationalv( VipsImage *left, VipsImage *right, VipsImage **out, * one-band image by joining n copies of the one-band image together, and then * the two n-band images are operated upon. * - * The two input images are cast up to the smallest common type (see table + * The two input images are cast up to the smallest common format (see table * Smallest common format in * arithmetic). * diff --git a/libvips/arithmetic/remainder.c b/libvips/arithmetic/remainder.c index 82a46edc..f5ad0cfa 100644 --- a/libvips/arithmetic/remainder.c +++ b/libvips/arithmetic/remainder.c @@ -210,7 +210,7 @@ vips_remainder_init( VipsRemainder *remainder ) * one-band image by joining n copies of the one-band image together, and then * the two n-band images are operated upon. * - * The two input images are cast up to the smallest common type (see table + * The two input images are cast up to the smallest common format (see table * Smallest common format in * arithmetic), and that format is the * result type. diff --git a/libvips/arithmetic/subtract.c b/libvips/arithmetic/subtract.c index 46a92dc0..a5f6223d 100644 --- a/libvips/arithmetic/subtract.c +++ b/libvips/arithmetic/subtract.c @@ -188,7 +188,7 @@ vips_subtract_init( VipsSubtract *subtract ) * one-band image by joining n copies of the one-band image together, and then * the two n-band images are operated upon. * - * The two input images are cast up to the smallest common type (see table + * The two input images are cast up to the smallest common format (see table * Smallest common format in * arithmetic), then the * following table is used to determine the output type: diff --git a/libvips/conversion/embed.c b/libvips/conversion/embed.c index 46a2d557..408aec9e 100644 --- a/libvips/conversion/embed.c +++ b/libvips/conversion/embed.c @@ -615,7 +615,7 @@ vips_embed_init( VipsEmbed *embed ) * * Optional arguments: * - * @extend: how to generate the edge pixels + * @extend: how to generate the edge pixels (default: black) * @background: colour for edge pixels * * The opposite of vips_extract_area(): embed @in within an image of size diff --git a/libvips/convolution/Makefile.am b/libvips/convolution/Makefile.am index ddc1aca3..9c543d89 100644 --- a/libvips/convolution/Makefile.am +++ b/libvips/convolution/Makefile.am @@ -3,18 +3,18 @@ noinst_LTLIBRARIES = libconvolution.la libconvolution_la_SOURCES = \ convolution.c \ pconvolution.h \ + correlation.c \ + correlation.h \ conv.c \ convsep.c \ compass.c \ morph.c \ - convol_dispatch.c \ + fastcor.c \ + spcor.c \ im_aconv.c \ im_aconvsep.c \ im_conv.c \ im_conv_f.c \ - im_fastcor.c \ - im_gradcor.c \ - im_sharpen.c \ - im_spcor.c + im_sharpen.c AM_CPPFLAGS = -I${top_srcdir}/libvips/include @VIPS_CFLAGS@ @VIPS_INCLUDES@ diff --git a/libvips/convolution/convolution.c b/libvips/convolution/convolution.c index 12abd08f..dd488d63 100644 --- a/libvips/convolution/convolution.c +++ b/libvips/convolution/convolution.c @@ -150,9 +150,13 @@ vips_convolution_operation_init( void ) extern int vips_morph_get_type( void ); extern int vips_compass_get_type( void ); extern int vips_convsep_get_type( void ); + extern int vips_fastcor_get_type( void ); + extern int vips_spcor_get_type( void ); vips_conv_get_type(); vips_morph_get_type(); vips_compass_get_type(); vips_convsep_get_type(); + vips_fastcor_get_type(); + vips_spcor_get_type(); } diff --git a/libvips/convolution/correlation.c b/libvips/convolution/correlation.c new file mode 100644 index 00000000..e1a3facf --- /dev/null +++ b/libvips/convolution/correlation.c @@ -0,0 +1,165 @@ +/* base class for correlation + * + * 7/11/13 + * - from convolution.c + */ + +/* + + 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., 51 Franklin Street, Fifth Floor, Boston, MA + 02110-1301 USA + + */ + +/* + + These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk + + */ + +#ifdef HAVE_CONFIG_H +#include +#endif /*HAVE_CONFIG_H*/ +#include + +#include +#include + +#include +#include + +#include "pconvolution.h" +#include "correlation.h" + +G_DEFINE_ABSTRACT_TYPE( VipsCorrelation, vips_correlation, + VIPS_TYPE_OPERATION ); + +static int +vips_correlation_gen( VipsRegion *or, + void *seq, void *a, void *b, gboolean *stop ) +{ + VipsRegion *ir = (VipsRegion *) seq; + VipsCorrelation *correlation = (VipsCorrelation *) b; + VipsCorrelationClass *cclass = + VIPS_CORRELATION_GET_CLASS( correlation ); + VipsRect *r = &or->valid; + + VipsRect irect; + + /* What part of ir do we need? + */ + irect.left = r->left; + irect.top = r->top; + irect.width = r->width + correlation->ref_ready->Xsize - 1; + irect.height = r->height + correlation->ref_ready->Ysize - 1; + + if( vips_region_prepare( ir, &irect ) ) + return( -1 ); + + cclass->correlation( correlation, ir, or ); + + return( 0 ); +} + +static int +vips_correlation_build( VipsObject *object ) +{ + VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( object ); + VipsCorrelationClass *cclass = VIPS_CORRELATION_CLASS( class ); + VipsCorrelation *correlation = (VipsCorrelation *) object; + VipsImage **t = (VipsImage **) vips_object_local_array( object, 5 ); + + if( VIPS_OBJECT_CLASS( vips_correlation_parent_class )-> + build( object ) ) + return( -1 ); + + /* Stretch input out. + */ + if( vips_embed( correlation->in, &t[0], + correlation->ref->Xsize / 2, correlation->ref->Ysize / 2, + correlation->in->Xsize + correlation->ref->Xsize - 1, + correlation->in->Ysize + correlation->ref->Ysize - 1, + "extend", VIPS_EXTEND_COPY, + NULL ) || + vips__formatalike( t[0], correlation->ref, &t[1], &t[2] ) || + vips__bandalike( class->nickname, t[1], t[2], &t[3], &t[4] ) || + vips_image_wio_input( t[4] ) ) + return( -1 ); + correlation->in_ready = t[3]; + correlation->ref_ready = t[4]; + + g_object_set( object, "out", vips_image_new(), NULL ); + + /* FATSTRIP is good for us as THINSTRIP will cause + * too many recalculations on overlaps. + */ + if( vips_image_pipelinev( correlation->out, + VIPS_DEMAND_STYLE_FATSTRIP, + correlation->in_ready, correlation->ref_ready, NULL ) ) + return( -1 ); + correlation->out->Xsize = correlation->in->Xsize; + correlation->out->Ysize = correlation->in->Ysize; + correlation->out->BandFmt = + cclass->format_table[correlation->in_ready->BandFmt]; + if( cclass->pre_generate && + cclass->pre_generate( correlation ) ) + return( -1 ); + if( vips_image_generate( correlation->out, + vips_start_one, vips_correlation_gen, vips_stop_one, + correlation->in_ready, correlation ) ) + return( -1 ); + + return( 0 ); +} + +static void +vips_correlation_class_init( VipsCorrelationClass *class ) +{ + GObjectClass *gobject_class = G_OBJECT_CLASS( class ); + VipsObjectClass *object_class = (VipsObjectClass *) class; + + gobject_class->set_property = vips_object_set_property; + gobject_class->get_property = vips_object_get_property; + + object_class->nickname = "correlation"; + object_class->description = _( "correlation operation" ); + object_class->build = vips_correlation_build; + + VIPS_ARG_IMAGE( class, "in", 0, + _( "Input" ), + _( "Input image argument" ), + VIPS_ARGUMENT_REQUIRED_INPUT, + G_STRUCT_OFFSET( VipsCorrelation, in ) ); + + VIPS_ARG_IMAGE( class, "ref", 10, + _( "Mask" ), + _( "Input reference image" ), + VIPS_ARGUMENT_REQUIRED_INPUT, + G_STRUCT_OFFSET( VipsCorrelation, ref ) ); + + VIPS_ARG_IMAGE( class, "out", 20, + _( "Output" ), + _( "Output image" ), + VIPS_ARGUMENT_REQUIRED_OUTPUT, + G_STRUCT_OFFSET( VipsCorrelation, out ) ); + +} + +static void +vips_correlation_init( VipsCorrelation *correlation ) +{ +} diff --git a/libvips/convolution/correlation.h b/libvips/convolution/correlation.h new file mode 100644 index 00000000..65db60d9 --- /dev/null +++ b/libvips/convolution/correlation.h @@ -0,0 +1,96 @@ +/* base class for all correlation operations + */ + +/* + + Copyright (C) 1991-2005 The National Gallery + + This library 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.1 of the License, or (at your option) any later version. + + This library 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 library; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA + 02110-1301 USA + + */ + +/* + + These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk + + */ + +#ifndef VIPS_PCORRELATION_H +#define VIPS_PCORRELATION_H + +#ifdef __cplusplus +extern "C" { +#endif /*__cplusplus*/ + +#include + +#define VIPS_TYPE_CORRELATION (vips_correlation_get_type()) +#define VIPS_CORRELATION( obj ) \ + (G_TYPE_CHECK_INSTANCE_CAST( (obj), \ + VIPS_TYPE_CORRELATION, VipsCorrelation )) +#define VIPS_CORRELATION_CLASS( klass ) \ + (G_TYPE_CHECK_CLASS_CAST( (klass), \ + VIPS_TYPE_CORRELATION, VipsCorrelationClass)) +#define VIPS_IS_CORRELATION( obj ) \ + (G_TYPE_CHECK_INSTANCE_TYPE( (obj), VIPS_TYPE_CORRELATION )) +#define VIPS_IS_CORRELATION_CLASS( klass ) \ + (G_TYPE_CHECK_CLASS_TYPE( (klass), VIPS_TYPE_CORRELATION )) +#define VIPS_CORRELATION_GET_CLASS( obj ) \ + (G_TYPE_INSTANCE_GET_CLASS( (obj), \ + VIPS_TYPE_CORRELATION, VipsCorrelationClass )) + +typedef struct { + VipsOperation parent_instance; + + /* Params. + */ + VipsImage *in; + VipsImage *ref; + VipsImage *out; + + /* The two input images, upcast to the smallest common format. ref is + * a memory buffer. + */ + VipsImage *in_ready; + VipsImage *ref_ready; + +} VipsCorrelation; + +typedef struct { + VipsOperationClass parent_class; + + /* For each upcast input format, what output format. + */ + const VipsBandFormat *format_table; + + /* Run just before generate. The subclass can fill in some stuff. + */ + int (*pre_generate)( VipsCorrelation * ); + + void (*correlation)( VipsCorrelation *, + VipsRegion *in, VipsRegion *out ); + +} VipsCorrelationClass; + +GType vips_correlation_get_type( void ); + +#ifdef __cplusplus +} +#endif /*__cplusplus*/ + +#endif /*VIPS_PCORRELATION_H*/ + + diff --git a/libvips/convolution/fastcor.c b/libvips/convolution/fastcor.c new file mode 100644 index 00000000..8740e23f --- /dev/null +++ b/libvips/convolution/fastcor.c @@ -0,0 +1,267 @@ +/* fastcor + * + * Copyright: 1990, N. Dessipris. + * + * Author: Nicos Dessipris + * Written on: 02/05/1990 + * Modified on : 15/03/1991 + * 20/2/95 JC + * - ANSIfied + * - in1 and in2 swapped, to match order for im_spcor + * - memory leaks fixed + * 21/2/95 JC + * - partialed + * - speed-ups + * 7/4/04 + * - now uses im_embed() with edge stretching on the output + * - sets Xoffset / Yoffset + * 8/3/06 JC + * - use im_embed() with edge stretching on the input, not the output + * - calculate sum of squares of differences, rather than abs of + * difference + * 3/2/10 + * - gtkdoc + * - cleanups + * 7/11/13 + * - redone 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., 51 Franklin Street, Fifth Floor, Boston, MA + 02110-1301 USA + + */ + +/* + + These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk + + */ + +#ifdef HAVE_CONFIG_H +#include +#endif /*HAVE_CONFIG_H*/ +#include + +#include +#include + +#include + +#include "pconvolution.h" +#include "correlation.h" + +typedef VipsCorrelationClass VipsFastcor; +typedef VipsCorrelationClass VipsFastcorClass; + +G_DEFINE_TYPE( VipsFastcor, vips_fastcor, VIPS_TYPE_CORRELATION ); + +#define CORR_INT( TYPE ) { \ + for( y = 0; y < r->height; y++ ) { \ + unsigned int *q = (unsigned int *) \ + VIPS_REGION_ADDR( out, r->left, r->top + y ); \ + \ + for( x = 0; x < r->width; x++ ) \ + for( b = 0; b < bands; b++ ) { \ + TYPE *p1 = (TYPE *) ref->data; \ + TYPE *p2 = (TYPE *) VIPS_REGION_ADDR( in, \ + r->left + x, r->top + y ); \ + \ + unsigned int sum; \ + \ + sum = 0; \ + for( j = 0; j < ref->Ysize; j++ ) { \ + for( i = b; i < sz; i += bands ) { \ + int t = p1[i] - p2[i]; \ + \ + sum += t * t; \ + } \ + \ + p1 += sz; \ + p2 += lsk; \ + } \ + \ + *q++ = sum; \ + } \ + } \ +} + +#define CORR_FLOAT( TYPE ) { \ + for( y = 0; y < r->height; y++ ) { \ + TYPE *q = (TYPE *) \ + VIPS_REGION_ADDR( out, r->left, r->top + y ); \ + \ + for( x = 0; x < r->width; x++ ) \ + for( b = 0; b < bands; b++ ) { \ + TYPE *p1 = (TYPE *) ref->data; \ + TYPE *p2 = (TYPE *) VIPS_REGION_ADDR( in, \ + r->left + x, r->top + y ); \ + \ + TYPE sum; \ + \ + sum = 0; \ + for( j = 0; j < ref->Ysize; j++ ) { \ + for( i = b; i < sz; i += bands ) { \ + TYPE t = p1[i] - p2[i]; \ + \ + sum += t * t; \ + } \ + \ + p1 += sz; \ + p2 += lsk; \ + } \ + \ + *q++ = sum; \ + } \ + } \ +} + +static void +vips_fastcor_correlation( VipsCorrelation *correlation, + VipsRegion *in, VipsRegion *out ) +{ + VipsRect *r = &out->valid; + VipsImage *ref = correlation->ref_ready; + int bands = vips_band_format_iscomplex( ref->BandFmt ) ? + ref->Bands * 2 : ref->Bands; + int sz = ref->Xsize * bands; + int lsk = VIPS_REGION_LSKIP( in ); + + int x, y, i, j, b; + + switch( vips_image_get_format( ref ) ) { + case VIPS_FORMAT_CHAR: + CORR_INT( signed char ); + break; + + case VIPS_FORMAT_UCHAR: + CORR_INT( unsigned char ); + break; + + case VIPS_FORMAT_SHORT: + CORR_INT( signed short ); + break; + + case VIPS_FORMAT_USHORT: + CORR_INT( unsigned short ); + break; + + case VIPS_FORMAT_INT: + CORR_INT( signed int ); + break; + + case VIPS_FORMAT_UINT: + CORR_INT( unsigned int ); + break; + + case VIPS_FORMAT_FLOAT: + case VIPS_FORMAT_COMPLEX: + CORR_FLOAT( float ); + break; + + case VIPS_FORMAT_DOUBLE: + case VIPS_FORMAT_DPCOMPLEX: + CORR_FLOAT( double ); + break; + + default: + g_assert( 0 ); + } +} + +/* Save a bit of typing. + */ +#define UC VIPS_FORMAT_UCHAR +#define C VIPS_FORMAT_CHAR +#define US VIPS_FORMAT_USHORT +#define S VIPS_FORMAT_SHORT +#define UI VIPS_FORMAT_UINT +#define I VIPS_FORMAT_INT +#define F VIPS_FORMAT_FLOAT +#define X VIPS_FORMAT_COMPLEX +#define D VIPS_FORMAT_DOUBLE +#define DX VIPS_FORMAT_DPCOMPLEX + +/* Type promotion for multiplication. Sign and value preserving. Make sure + * these match the case statement in multiply_buffer() above. + */ +static int vips_fastcor_format_table[10] = { +/* UC C US S UI I F X D DX */ + UI, UI, UI, UI, UI, UI,F, X, D, DX +}; + +static void +vips_fastcor_class_init( VipsFastcorClass *class ) +{ + VipsObjectClass *object_class = (VipsObjectClass *) class; + VipsCorrelationClass *cclass = VIPS_CORRELATION_CLASS( class ); + + object_class->nickname = "fastcor"; + object_class->description = _( "fast correlation" ); + + cclass->format_table = vips_fastcor_format_table; + cclass->correlation = vips_fastcor_correlation; +} + +static void +vips_fastcor_init( VipsFastcor *fastcor ) +{ +} + +/** + * vips_fastcor: + * @in: input image + * @ref: reference image + * @out: output image + * @...: %NULL-terminated list of optional named arguments + * + * Calculate a fast correlation surface. + * + * @ref is placed at every position in @in and the sum of squares of + * differences calculated. + * + * The output + * image is the same size as the input. Extra input edge pixels are made by + * copying the existing edges outwards. + * + * If the number of bands differs, one of the images + * must have one band. In this case, an n-band image is formed from the + * one-band image by joining n copies of the one-band image together, and then + * the two n-band images are operated upon. + * + * The output type is uint if both inputs are integer, float if both are float + * or complex, and double if either is double or double complex. + * In other words, the output type is just large enough to hold the whole + * range of possible values. + * + * See also: vips_spcor(). + * + * Returns: 0 on success, -1 on error + */ +int +vips_fastcor( VipsImage *in, VipsImage *ref, VipsImage **out, ... ) +{ + va_list ap; + int result; + + va_start( ap, out ); + result = vips_call_split( "fastcor", ap, in, ref, out ); + va_end( ap ); + + return( result ); +} diff --git a/libvips/convolution/im_fastcor.c b/libvips/convolution/im_fastcor.c deleted file mode 100644 index 91962687..00000000 --- a/libvips/convolution/im_fastcor.c +++ /dev/null @@ -1,206 +0,0 @@ -/* im_fastcor - * - * Copyright: 1990, N. Dessipris. - * - * Author: Nicos Dessipris - * Written on: 02/05/1990 - * Modified on : 15/03/1991 - * 20/2/95 JC - * - ANSIfied - * - in1 and in2 swapped, to match order for im_spcor - * - memory leaks fixed - * 21/2/95 JC - * - partialed - * - speed-ups - * 7/4/04 - * - now uses im_embed() with edge stretching on the output - * - sets Xoffset / Yoffset - * 8/3/06 JC - * - use im_embed() with edge stretching on the input, not the output - * - calculate sum of squares of differences, rather than abs of - * difference - * 3/2/10 - * - gtkdoc - * - cleanups - */ - -/* - - 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., 51 Franklin Street, Fifth Floor, Boston, MA - 02110-1301 USA - - */ - -/* - - These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk - - */ - -#ifdef HAVE_CONFIG_H -#include -#endif /*HAVE_CONFIG_H*/ -#include - -#include -#include - -#include - -/* Fastcor generate function. - */ -static int -fastcor_gen( REGION *or, void *seq, void *a, void *b ) -{ - REGION *ir = (REGION *) seq; - IMAGE *ref = (IMAGE *) b; - Rect irect; - Rect *r = &or->valid; - - int x, y, i, j; - int lsk; - - /* What part of ir do we need? - */ - irect.left = or->valid.left; - irect.top = or->valid.top; - irect.width = or->valid.width + ref->Xsize - 1; - irect.height = or->valid.height + ref->Ysize - 1; - - if( im_prepare( ir, &irect ) ) - return( -1 ); - lsk = IM_REGION_LSKIP( ir ); - - /* Loop over or. - */ - for( y = 0; y < r->height; y++ ) { - unsigned int *q = (unsigned int *) - IM_REGION_ADDR( or, r->left, r->top + y ); - - for( x = 0; x < r->width; x++ ) { - VipsPel *b = ref->data; - VipsPel *a = - IM_REGION_ADDR( ir, r->left + x, r->top + y ); - - int sum; - - sum = 0; - for( j = 0; j < ref->Ysize; j++ ) { - for( i = 0; i < ref->Xsize; i++ ) { - int t = b[i] - a[i]; - - sum += t * t; - } - - a += lsk; - b += ref->Xsize; - } - - q[x] = sum; - } - } - - return( 0 ); -} - -/* Raw fastcor, with no borders. - */ -int -im_fastcor_raw( IMAGE *in, IMAGE *ref, IMAGE *out ) -{ - /* PIO between in and out; WIO from ref. - */ - if( im_piocheck( in, out ) || - im_incheck( ref ) ) - return( -1 ); - - /* Check sizes. - */ - if( in->Xsize < ref->Xsize || in->Ysize < ref->Ysize ) { - im_error( "im_fastcor", "%s", - _( "ref not smaller than or equal to in" ) ); - return( -1 ); - } - - /* Check types. - */ - if( im_check_uncoded( "im_fastcor", in ) || - im_check_mono( "im_fastcor", in ) || - im_check_format( "im_fastcor", in, IM_BANDFMT_UCHAR ) || - im_check_coding_same( "im_fastcor", in, ref ) || - im_check_bands_same( "im_fastcor", in, ref ) || - im_check_format_same( "im_fastcor", in, ref ) ) - return( -1 ); - - /* Prepare the output image. - */ - if( im_cp_descv( out, in, ref, NULL ) ) - return( -1 ); - out->BandFmt = IM_BANDFMT_UINT; - out->Xsize = in->Xsize - ref->Xsize + 1; - out->Ysize = in->Ysize - ref->Ysize + 1; - - /* FATSTRIP is good for us, as THINSTRIP will cause - * too many recalculations on overlaps. - */ - if( im_demand_hint( out, IM_FATSTRIP, in, NULL ) || - im_generate( out, - im_start_one, fastcor_gen, im_stop_one, in, ref ) ) - return( -1 ); - - out->Xoffset = -ref->Xsize / 2; - out->Yoffset = -ref->Ysize / 2; - - return( 0 ); -} - -/** - * im_fastcor: - * @in: input image - * @ref: reference image - * @out: output image - * - * Calculate a fast correlation surface. - * - * @ref is placed at every position in @in and the sum of squares of - * differences calculated. One-band, 8-bit unsigned images only. The output - * image is always %IM_BANDFMT_UINT. @ref must be smaller than or equal to - * @in. The output - * image is the same size as the input. - * - * See also: im_spcor(). - * - * Returns: 0 on success, -1 on error - */ -int -im_fastcor( IMAGE *in, IMAGE *ref, IMAGE *out ) -{ - IMAGE *t1 = im_open_local( out, "im_fastcor intermediate", "p" ); - - if( !t1 || - im_embed( in, t1, 1, - ref->Xsize / 2, ref->Ysize / 2, - in->Xsize + ref->Xsize - 1, - in->Ysize + ref->Ysize - 1 ) || - im_fastcor_raw( t1, ref, out ) ) - return( -1 ); - - out->Xoffset = 0; - out->Yoffset = 0; - - return( 0 ); -} diff --git a/libvips/convolution/im_sharpen.c b/libvips/convolution/im_sharpen.c index aab29405..97ccad14 100644 --- a/libvips/convolution/im_sharpen.c +++ b/libvips/convolution/im_sharpen.c @@ -96,10 +96,10 @@ typedef struct { static SharpenLut * build_lut( IMAGE *out, int x1, int x2, int x3, double m1, double m2 ) { + SharpenLut *slut; int i; - SharpenLut *slut = IM_NEW( out, SharpenLut ); - if( !slut ) + if( !(slut = IM_NEW( out, SharpenLut )) ) return( NULL ); if( !(slut->lut = IM_ARRAY( out, x2 + x3 + 1, int )) ) @@ -109,13 +109,13 @@ build_lut( IMAGE *out, int x1, int x2, int x3, double m1, double m2 ) slut->x3 = x3; for( i = 0; i < x1; i++ ) { - slut->lut[x3 + i] = i*m1; - slut->lut[x3 - i] = -i*m1; + slut->lut[x3 + i] = i * m1; + slut->lut[x3 - i] = -i * m1; } for( i = x1; i <= x2; i++ ) - slut->lut[x3 + i] = x1*m1 + (i-x1)*m2; + slut->lut[x3 + i] = x1 * m1 + (i - x1) * m2; for( i = x1; i <= x3; i++ ) - slut->lut[x3 - i] = -(x1*m1 + (i-x1)*m2); + slut->lut[x3 - i] = -(x1 * m1 + (i - x1) * m2); return( slut ); } diff --git a/libvips/convolution/im_spcor.c b/libvips/convolution/im_spcor.c deleted file mode 100644 index c2448bd9..00000000 --- a/libvips/convolution/im_spcor.c +++ /dev/null @@ -1,348 +0,0 @@ -/* im_spcor - * - * Copyright: 1990, N. Dessipris; 2006, 2007 Nottingham Trent University. - * - * - * Author: Nicos Dessipris - * Written on: 02/05/1990 - * Modified on : - * 20/2/95 JC - * - updated - * - ANSIfied, a little - * 21/2/95 JC - * - rewritten - * - partialed - * - speed-ups - * - new correlation coefficient (see above), from Niblack "An - * Introduction to Digital Image Processing", Prentice/Hall, pp 138. - * 4/9/97 JC - * - now does short/ushort as well - * 13/2/03 JC - * - oops, could segv for short images - * 14/4/04 JC - * - sets Xoffset / Yoffset - * 8/3/06 JC - * - use im_embed() with edge stretching on the input, not the output - * - * 2006-10-24 tcv - * - add im_spcor2 - * - * 2007-11-12 tcv - * - make im_spcor a wrapper selecting either im__spcor or im__spcor2 - * 2008-09-09 JC - * - roll back the windowed version for now, it has some tile edge effects - * 3/2/10 - * - gtkdoc - * - cleanups - */ - -/* - - 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., 51 Franklin Street, Fifth Floor, Boston, MA - 02110-1301 USA - - */ - -/* - - These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk - - */ - -#ifdef HAVE_CONFIG_H -#include -#endif /*HAVE_CONFIG_H*/ -#include - -#include -#include -#include -#include - -#include - -/* Hold per-call state here. - */ -typedef struct { - IMAGE *ref; /* Image we are searching for */ - double rmean; /* Mean of search window */ - double c1; /* sqrt(sumij (ref(i,j)-mean(ref))^2) */ -} Spcor; - -#define LOOP( IN ) { \ - IN *a = (IN *) p; \ - IN *b = (IN *) ref->data; \ - int in_lsk = lsk / sizeof( IN ); \ - IN *a1; \ - IN *b1; \ - \ - /* For each pel in or, loop over ref. First, \ - * calculate mean of area in ir corresponding to ref. \ - */ \ - a1 = a; \ - sum1 = 0; \ - for( j = 0; j < ref->Ysize; j++ ) { \ - for( i = 0; i < ref->Xsize; i++ ) \ - sum1 += a1[i]; \ - a1 += in_lsk; \ - } \ - imean = sum1 / VIPS_IMAGE_N_PELS( ref ); \ - \ - /* Loop over ir again, this time calculating \ - * sum-of-squares-of-differences for this window on \ - * ir, and also sum-of-products-of-differences from mean. \ - */ \ - a1 = a; \ - b1 = b; \ - sum2 = 0.0; \ - sum3 = 0.0; \ - for( j = 0; j < ref->Ysize; j++ ) { \ - for( i = 0; i < ref->Xsize; i++ ) { \ - /* Reference pel, and input pel. \ - */ \ - IN rp = b1[i]; \ - IN ip = a1[i]; \ - \ - /* Accumulate sum-of-squares-of- \ - * differences for input image. \ - */ \ - double t = ip - imean; \ - sum2 += t * t; \ - \ - /* Accumulate product-of-difference from mean. \ - */ \ - sum3 += (rp - spcor->rmean) * (ip - imean); \ - } \ - a1 += in_lsk; \ - b1 += ref->Xsize; \ - } \ -} - -/* spcor generate function. - */ -static int -spcor_gen( REGION *or, void *vseq, void *a, void *b ) -{ - REGION *ir = (REGION *) vseq; - Spcor *spcor = (Spcor *) b; - IMAGE *ref = spcor->ref; - Rect irect; - Rect *r = &or->valid; - int le = r->left; - int to = r->top; - int bo = IM_RECT_BOTTOM( r ); - int ri = IM_RECT_RIGHT( r ); - - int x, y, i, j; - int lsk; - - double imean; - double sum1; - double sum2, sum3; - double c2, cc; - - /* What part of ir do we need? - */ - irect.left = or->valid.left; - irect.top = or->valid.top; - irect.width = or->valid.width + ref->Xsize - 1; - irect.height = or->valid.height + ref->Ysize - 1; - - if( im_prepare( ir, &irect ) ) - return( -1 ); - lsk = IM_REGION_LSKIP( ir ); - - /* Loop over or. - */ - for( y = to; y < bo; y++ ) { - float *q = (float *) IM_REGION_ADDR( or, le, y ); - - for( x = le; x < ri; x++ ) { - VipsPel *p = IM_REGION_ADDR( ir, x, y ); - - /* Find sums for this position. - */ - switch( ref->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; - default: - g_assert( 0 ); - - /* Keep -Wall happy. - */ - return( 0 ); - } - - /* Now: calculate correlation coefficient! - */ - c2 = sqrt( sum2 ); - cc = sum3 / (spcor->c1 * c2); - - *q++ = cc; - } - } - - return( 0 ); -} - -/* Pre-calculate stuff for our reference image. - */ -static Spcor * -spcor_new( IMAGE *out, IMAGE *ref ) -{ - Spcor *spcor; - guint64 sz = VIPS_IMAGE_N_PELS( ref ); - VipsPel *p = ref->data; - double s; - size_t i; - - if( !(spcor = IM_NEW( out, Spcor )) ) - return( NULL ); - - /* Pre-calculate stuff on our reference image. - */ - spcor->ref = ref; - if( im_avg( spcor->ref, &spcor->rmean ) ) - return( NULL ); - - /* Find sqrt-of-sum-of-squares-of-differences. - */ - s = 0.0; - for( i = 0; i < sz; i++ ) { - double t = (int) p[i] - spcor->rmean; - s += t * t; - } - spcor->c1 = sqrt( s ); - - return( spcor ); -} - -int -im_spcor_raw( IMAGE *in, IMAGE *ref, IMAGE *out ) -{ - Spcor *spcor; - - /* PIO between in and out; WIO from ref, since it's probably tiny. - */ - if( im_piocheck( in, out ) || - im_incheck( ref ) ) - return( -1 ); - - /* Check sizes. - */ - if( in->Xsize < ref->Xsize || - in->Ysize < ref->Ysize ) { - im_error( "im_spcor_raw", - "%s", _( "ref not smaller than or equal to in" ) ); - return( -1 ); - } - - /* Check types. - */ - if( im_check_uncoded( "im_spcor", in ) || - im_check_mono( "im_spcor", in ) || - im_check_8or16( "im_spcor", in ) || - im_check_coding_same( "im_spcor", in, ref ) || - im_check_bands_same( "im_spcor", in, ref ) || - im_check_format_same( "im_spcor", in, ref ) ) - return( -1 ); - - /* Prepare the output image. - */ - if( im_cp_descv( out, in, ref, NULL ) ) - return( -1 ); - out->BandFmt = IM_BANDFMT_FLOAT; - out->Xsize = in->Xsize - ref->Xsize + 1; - out->Ysize = in->Ysize - ref->Ysize + 1; - - /* Pre-calculate some stuff. - */ - if( !(spcor = spcor_new( out, ref )) ) - return( -1 ); - - /* Set demand hints. FATSTRIP is good for us, as THINSTRIP will cause - * too many recalculations on overlaps. - */ - if( im_demand_hint( out, IM_FATSTRIP, in, NULL ) ) - return( -1 ); - - /* Write the correlation. - */ - if( im_generate( out, - im_start_one, spcor_gen, im_stop_one, in, spcor ) ) - return( -1 ); - - out->Xoffset = -ref->Xsize / 2; - out->Yoffset = -ref->Ysize / 2; - - return( 0 ); -} - -/** - * im_spcor: - * @in: input image - * @ref: reference image - * @out: output image - * - * Calculate a correlation surface. - * - * @ref is placed at every position in @in and the correlation coefficient - * calculated. One-band, 8 or 16-bit images only. @in and @ref must have the - * same #VipsBandFmt. The output - * image is always %IM_BANDFMT_FLOAT. @ref must be smaller than or equal to - * @in. The output - * image is the same size as the input. - * - * The correlation coefficient is calculated as: - * - * |[ - * sumij (ref(i,j)-mean(ref))(inkl(i,j)-mean(inkl)) - * c(k,l) = ------------------------------------------------ - * sqrt(sumij (ref(i,j)-mean(ref))^2) * - * sqrt(sumij (inkl(i,j)-mean(inkl))^2) - * ]| - * - * where inkl is the area of @in centred at position (k,l). - * - * from Niblack "An Introduction to Digital Image Processing", - * Prentice/Hall, pp 138. - * - * See also: im_gradcor(), im_fastcor(). - * - * Returns: 0 on success, -1 on error - */ -int -im_spcor( IMAGE *in, IMAGE *ref, IMAGE *out ) -{ - IMAGE *t1 = im_open_local( out, "im_spcor intermediate", "p" ); - - if( !t1 || - im_embed( in, t1, 1, - ref->Xsize / 2, ref->Ysize / 2, - in->Xsize + ref->Xsize - 1, - in->Ysize + ref->Ysize - 1 ) || - im_spcor_raw( t1, ref, out ) ) - return( -1 ); - - out->Xoffset = 0; - out->Yoffset = 0; - - return( 0 ); -} - diff --git a/libvips/convolution/spcor.c b/libvips/convolution/spcor.c new file mode 100644 index 00000000..edf854b7 --- /dev/null +++ b/libvips/convolution/spcor.c @@ -0,0 +1,367 @@ +/* spcor + * + * Copyright: 1990, N. Dessipris; 2006, 2007 Nottingham Trent University. + * + * + * Author: Nicos Dessipris + * Written on: 02/05/1990 + * Modified on : + * 20/2/95 JC + * - updated + * - ANSIfied, a little + * 21/2/95 JC + * - rewritten + * - partialed + * - speed-ups + * - new correlation coefficient (see above), from Niblack "An + * Introduction to Digital Image Processing", Prentice/Hall, pp 138. + * 4/9/97 JC + * - now does short/ushort as well + * 13/2/03 JC + * - oops, could segv for short images + * 14/4/04 JC + * - sets Xoffset / Yoffset + * 8/3/06 JC + * - use im_embed() with edge stretching on the input, not the output + * + * 2006-10-24 tcv + * - add im_spcor2 + * + * 2007-11-12 tcv + * - make im_spcor a wrapper selecting either im__spcor or im__spcor2 + * 2008-09-09 JC + * - roll back the windowed version for now, it has some tile edge effects + * 3/2/10 + * - gtkdoc + * - cleanups + * 7/11/13 + * - redone 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., 51 Franklin Street, Fifth Floor, Boston, MA + 02110-1301 USA + + */ + +/* + + These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk + + */ + +#ifdef HAVE_CONFIG_H +#include +#endif /*HAVE_CONFIG_H*/ +#include + +#include +#include + +#include + +#include "pconvolution.h" +#include "correlation.h" + +typedef struct _VipsSpcor { + VipsCorrelation parent_instance; + + /* Per-band mean of ref images. + */ + double *rmean; + + /* Per band sqrt(sumij (ref(i,j)-mean(ref))^2) + */ + double *c1; +} VipsSpcor; + +typedef VipsCorrelationClass VipsSpcorClass; + +G_DEFINE_TYPE( VipsSpcor, vips_spcor, VIPS_TYPE_CORRELATION ); + +static int +vips_spcor_pre_generate( VipsCorrelation *correlation ) +{ + VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( correlation ); + VipsSpcor *spcor = (VipsSpcor *) correlation; + VipsImage *ref = correlation->ref_ready; + int bands = ref->Bands; + VipsImage **b = (VipsImage **) + vips_object_local_array( VIPS_OBJECT( spcor ), bands ); + VipsImage **t = (VipsImage **) + vips_object_local_array( VIPS_OBJECT( spcor ), 2 ); + VipsImage **b2 = (VipsImage **) + vips_object_local_array( VIPS_OBJECT( spcor ), bands ); + + int i; + double *offset; + double *scale; + + if( vips_check_noncomplex( class->nickname, ref ) ) + return( -1 ); + + /* Per-band mean. + */ + if( !(spcor->rmean = VIPS_ARRAY( spcor, bands, double )) || + !(spcor->c1 = VIPS_ARRAY( spcor, bands, double )) ) + return( -1 ); + for( i = 0; i < bands; i++ ) + if( vips_extract_band( ref, &b[i], i, NULL ) || + vips_avg( b[i], &spcor->rmean[i], NULL ) ) + return( -1 ); + + /* Per band sqrt(sumij (ref(i,j)-mean(ref))^2) + */ + if( !(offset = VIPS_ARRAY( spcor, bands, double )) || + !(scale = VIPS_ARRAY( spcor, bands, double )) ) + return( -1 ); + for( i = 0; i < bands; i++ ) { + offset[i] = -spcor->rmean[i]; + scale[i] = 1.0; + } + if( vips_linear( ref, &t[0], scale, offset, bands, NULL ) || + vips_multiply( t[0], t[0], &t[1], NULL ) ) + return( -1 ); + for( i = 0; i < bands; i++ ) + if( vips_extract_band( t[1], &b2[i], i, NULL ) || + vips_avg( b2[i], &spcor->c1[i], NULL ) ) + return( -1 ); + for( i = 0; i < bands; i++ ) { + spcor->c1[i] *= ref->Xsize * ref->Ysize; + spcor->c1[i] = sqrt( spcor->c1[i] ); + } + + return( 0 ); +} + +#define LOOP( IN ) { \ + IN *r1 = ((IN *) ref->data) + b; \ + IN *p1 = ((IN *) p) + b; \ + int in_lsk = lsk / sizeof( IN ); \ + IN *r1a; \ + IN *p1a; \ + \ + /* Mean of area of in corresponding to ref. \ + */ \ + p1a = p1; \ + sum1 = 0.0; \ + for( j = 0; j < ref->Ysize; j++ ) { \ + for( i = 0; i < sz; i += bands ) \ + sum1 += p1a[i]; \ + p1a += in_lsk; \ + } \ + imean = sum1 / VIPS_IMAGE_N_PELS( ref ); \ + \ + /* Calculate sum-of-squares-of-differences for this window on \ + * in, and also sum-of-products-of-differences from mean. \ + */ \ + p1a = p1; \ + r1a = r1; \ + sum2 = 0.0; \ + sum3 = 0.0; \ + for( j = 0; j < ref->Ysize; j++ ) { \ + for( i = 0; i < sz; i += bands ) { \ + /* Reference pel, and input pel. \ + */ \ + IN ip = p1a[i]; \ + IN rp = r1a[i]; \ + \ + /* Accumulate sum-of-squares-of- \ + * differences for input image. \ + */ \ + double t = ip - imean; \ + sum2 += t * t; \ + \ + /* Accumulate product-of-difference from mean. \ + */ \ + sum3 += (rp - spcor->rmean[b]) * (ip - imean); \ + } \ + \ + p1a += in_lsk; \ + r1a += sz; \ + } \ +} + +static void +vips_spcor_correlation( VipsCorrelation *correlation, + VipsRegion *in, VipsRegion *out ) +{ + VipsSpcor *spcor = (VipsSpcor *) correlation; + VipsRect *r = &out->valid; + VipsImage *ref = correlation->ref_ready; + int bands = vips_band_format_iscomplex( ref->BandFmt ) ? + ref->Bands * 2 : ref->Bands; + int sz = ref->Xsize * bands; + int lsk = VIPS_REGION_LSKIP( in ); + + int x, y, b, j, i; + + double imean; + double sum1; + double sum2, sum3; + double c2, cc; + + for( y = 0; y < r->height; y++ ) { + float *q = (float *) + VIPS_REGION_ADDR( out, r->left, r->top + y ); + + for( x = 0; x < r->width; x++ ) { + VipsPel *p = + VIPS_REGION_ADDR( in, r->left + x, r->top + y ); + + for( b = 0; b < bands; b++ ) { + switch( vips_image_get_format( ref ) ) { + 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: + case VIPS_FORMAT_COMPLEX: + LOOP( float ); + break; + + case VIPS_FORMAT_DOUBLE: + case VIPS_FORMAT_DPCOMPLEX: + LOOP( double ); + break; + + default: + g_assert( 0 ); + } + + c2 = sqrt( sum2 ); + cc = sum3 / (spcor->c1[b] * c2); + + *q++ = cc; + } + } + } +} + +/* Save a bit of typing. + */ +#define UC VIPS_FORMAT_UCHAR +#define C VIPS_FORMAT_CHAR +#define US VIPS_FORMAT_USHORT +#define S VIPS_FORMAT_SHORT +#define UI VIPS_FORMAT_UINT +#define I VIPS_FORMAT_INT +#define F VIPS_FORMAT_FLOAT +#define X VIPS_FORMAT_COMPLEX +#define D VIPS_FORMAT_DOUBLE +#define DX VIPS_FORMAT_DPCOMPLEX + +static int vips_spcor_format_table[10] = { +/* UC C US S UI I F X D DX */ + F, F, F, F, F, F, F, F, F, F +}; + +static void +vips_spcor_class_init( VipsSpcorClass *class ) +{ + VipsObjectClass *object_class = (VipsObjectClass *) class; + VipsCorrelationClass *cclass = VIPS_CORRELATION_CLASS( class ); + + object_class->nickname = "spcor"; + object_class->description = _( "spatial correlation" ); + + cclass->format_table = vips_spcor_format_table; + cclass->pre_generate = vips_spcor_pre_generate; + cclass->correlation = vips_spcor_correlation; +} + +static void +vips_spcor_init( VipsSpcor *spcor ) +{ +} + +/** + * vips_spcor: + * @in: input image + * @ref: reference image + * @out: output image + * @...: %NULL-terminated list of optional named arguments + * + * Calculate a correlation surface. + * + * @ref is placed at every position in @in and the correlation coefficient + * calculated. The output + * image is always float. + * + * The output + * image is the same size as the input. Extra input edge pixels are made by + * copying the existing edges outwards. + * + * The correlation coefficient is calculated as: + * + * |[ + * sumij (ref(i,j)-mean(ref))(inkl(i,j)-mean(inkl)) + * c(k,l) = ------------------------------------------------ + * sqrt(sumij (ref(i,j)-mean(ref))^2) * + * sqrt(sumij (inkl(i,j)-mean(inkl))^2) + * ]| + * + * where inkl is the area of @in centred at position (k,l). + * + * from Niblack "An Introduction to Digital Image Processing", + * Prentice/Hall, pp 138. + * + * If the number of bands differs, one of the images + * must have one band. In this case, an n-band image is formed from the + * one-band image by joining n copies of the one-band image together, and then + * the two n-band images are operated upon. + * + * The output image is always float, unless either of the two inputs is + * double, in which case the output is also double. + * + * See also: vips_gradcor(), vips_fastcor(). + * + * Returns: 0 on success, -1 on error + */ +int +vips_spcor( VipsImage *in, VipsImage *ref, VipsImage **out, ... ) +{ + va_list ap; + int result; + + va_start( ap, out ); + result = vips_call_split( "spcor", ap, in, ref, out ); + va_end( ap ); + + return( result ); +} diff --git a/libvips/deprecated/Makefile.am b/libvips/deprecated/Makefile.am index d0110ab6..4e67065e 100644 --- a/libvips/deprecated/Makefile.am +++ b/libvips/deprecated/Makefile.am @@ -9,6 +9,7 @@ libdeprecated_la_SOURCES = \ im_lab_morph.c \ deprecated_dispatch.c \ colour_dispatch.c \ + convol_dispatch.c \ arith_dispatch.c \ hist_dispatch.c \ im_maxpos_avg.c \ @@ -24,6 +25,7 @@ libdeprecated_la_SOURCES = \ im_fav4.c \ im_gadd.c \ im_gaddim.c \ + im_gradcor.c \ im_cmulnorm.c \ im_printlines.c \ im_convsub.c \ diff --git a/libvips/convolution/convol_dispatch.c b/libvips/deprecated/convol_dispatch.c similarity index 100% rename from libvips/convolution/convol_dispatch.c rename to libvips/deprecated/convol_dispatch.c diff --git a/libvips/convolution/im_gradcor.c b/libvips/deprecated/im_gradcor.c similarity index 100% rename from libvips/convolution/im_gradcor.c rename to libvips/deprecated/im_gradcor.c diff --git a/libvips/deprecated/vips7compat.c b/libvips/deprecated/vips7compat.c index 8e89736b..266ad0c1 100644 --- a/libvips/deprecated/vips7compat.c +++ b/libvips/deprecated/vips7compat.c @@ -2278,6 +2278,52 @@ im_contrast_surface( IMAGE *in, IMAGE *out, int half_win_size, int spacing ) return( 0 ); } +int +im_spcor_raw( IMAGE *in, IMAGE *ref, IMAGE *out ) +{ + im_error( "im_spcor_raw", "no compat function" ); + return( -1 ); +} + +int +im_spcor( IMAGE *in, IMAGE *ref, IMAGE *out ) +{ + VipsImage *x; + + if( vips_call( "spcor", in, ref, &x, NULL ) ) + return( -1 ); + if( im_copy( x, out ) ) { + g_object_unref( x ); + return( -1 ); + } + g_object_unref( x ); + + return( 0 ); +} + +int +im_fastcor_raw( IMAGE *in, IMAGE *ref, IMAGE *out ) +{ + im_error( "im_fastcor_raw", "no compat function" ); + return( -1 ); +} + +int +im_fastcor( IMAGE *in, IMAGE *ref, IMAGE *out ) +{ + VipsImage *x; + + if( vips_call( "fastcor", in, ref, &x, NULL ) ) + return( -1 ); + if( im_copy( x, out ) ) { + g_object_unref( x ); + return( -1 ); + } + g_object_unref( x ); + + return( 0 ); +} + static int vips__round( VipsImage *in, VipsImage *out, VipsOperationRound round ) { diff --git a/libvips/foreign/radiance.c b/libvips/foreign/radiance.c index 0d90063f..a540f2ef 100644 --- a/libvips/foreign/radiance.c +++ b/libvips/foreign/radiance.c @@ -788,7 +788,6 @@ scanline_write( COLR *scanline, int width, FILE *fp ) } int i, j, beg, cnt; - int c2; if( width < MINELEN || width > MAXELEN ) diff --git a/libvips/include/vips/convolution.h b/libvips/include/vips/convolution.h index 798766b4..0f0fcf60 100644 --- a/libvips/include/vips/convolution.h +++ b/libvips/include/vips/convolution.h @@ -88,13 +88,6 @@ int im_sharpen( VipsImage *in, VipsImage *out, double x1, double y2, double y3, double m1, double m2 ); -int im_grad_x( VipsImage *in, VipsImage *out ); -int im_grad_y( VipsImage *in, VipsImage *out ); - -int im_fastcor( VipsImage *in, VipsImage *ref, VipsImage *out ); -int im_spcor( VipsImage *in, VipsImage *ref, VipsImage *out ); -int im_gradcor( VipsImage *in, VipsImage *ref, VipsImage *out ); - #ifdef __cplusplus } #endif /*__cplusplus*/ diff --git a/libvips/include/vips/vips7compat.h b/libvips/include/vips/vips7compat.h index cfeae70d..616f5039 100644 --- a/libvips/include/vips/vips7compat.h +++ b/libvips/include/vips/vips7compat.h @@ -923,6 +923,13 @@ int im_contrast_surface_raw( IMAGE *in, IMAGE *out, int im_contrast_surface( VipsImage *in, VipsImage *out, int half_win_size, int spacing ); +int im_grad_x( VipsImage *in, VipsImage *out ); +int im_grad_y( VipsImage *in, VipsImage *out ); + +int im_fastcor( VipsImage *in, VipsImage *ref, VipsImage *out ); +int im_spcor( VipsImage *in, VipsImage *ref, VipsImage *out ); +int im_gradcor( VipsImage *in, VipsImage *ref, VipsImage *out ); + #ifdef __cplusplus } #endif /*__cplusplus*/