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
This commit is contained in:
parent
9a06b2cea7
commit
f0d4760560
@ -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()
|
||||
|
12
TODO
12
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
|
||||
|
@ -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
|
||||
* <link linkend="VIPS-arithmetic">arithmetic</link>), then the
|
||||
* following table is used to determine the output type:
|
||||
|
@ -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
|
||||
* <link linkend="VIPS-arithmetic">arithmetic</link>).
|
||||
*
|
||||
|
@ -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
|
||||
* <link linkend="VIPS-arithmetic">arithmetic</link>), then the
|
||||
* following table is used to determine the output type:
|
||||
|
@ -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; \
|
||||
} \
|
||||
|
@ -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
|
||||
* <link linkend="VIPS-arithmetic">arithmetic</link>), and that format is the
|
||||
* result type.
|
||||
|
@ -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
|
||||
* <link linkend="VIPS-arithmetic">arithmetic</link>), then the
|
||||
* following table is used to determine the output type:
|
||||
|
@ -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
|
||||
* <link linkend="VIPS-arithmetic">arithmetic</link>).
|
||||
*
|
||||
|
@ -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
|
||||
* <link linkend="VIPS-arithmetic">arithmetic</link>), and that format is the
|
||||
* result type.
|
||||
|
@ -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
|
||||
* <link linkend="VIPS-arithmetic">arithmetic</link>), then the
|
||||
* following table is used to determine the output type:
|
||||
|
@ -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
|
||||
|
@ -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@
|
||||
|
@ -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();
|
||||
}
|
||||
|
165
libvips/convolution/correlation.c
Normal file
165
libvips/convolution/correlation.c
Normal file
@ -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 <config.h>
|
||||
#endif /*HAVE_CONFIG_H*/
|
||||
#include <vips/intl.h>
|
||||
|
||||
#include <stdio.h>
|
||||
#include <math.h>
|
||||
|
||||
#include <vips/vips.h>
|
||||
#include <vips/internal.h>
|
||||
|
||||
#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 )
|
||||
{
|
||||
}
|
96
libvips/convolution/correlation.h
Normal file
96
libvips/convolution/correlation.h
Normal file
@ -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 <vips/vector.h>
|
||||
|
||||
#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*/
|
||||
|
||||
|
267
libvips/convolution/fastcor.c
Normal file
267
libvips/convolution/fastcor.c
Normal file
@ -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 <config.h>
|
||||
#endif /*HAVE_CONFIG_H*/
|
||||
#include <vips/intl.h>
|
||||
|
||||
#include <stdio.h>
|
||||
#include <math.h>
|
||||
|
||||
#include <vips/vips.h>
|
||||
|
||||
#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 );
|
||||
}
|
@ -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 <config.h>
|
||||
#endif /*HAVE_CONFIG_H*/
|
||||
#include <vips/intl.h>
|
||||
|
||||
#include <stdio.h>
|
||||
#include <math.h>
|
||||
|
||||
#include <vips/vips.h>
|
||||
|
||||
/* 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 );
|
||||
}
|
@ -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 );
|
||||
}
|
||||
|
@ -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 <config.h>
|
||||
#endif /*HAVE_CONFIG_H*/
|
||||
#include <vips/intl.h>
|
||||
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <stdio.h>
|
||||
#include <math.h>
|
||||
|
||||
#include <vips/vips.h>
|
||||
|
||||
/* 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 );
|
||||
}
|
||||
|
367
libvips/convolution/spcor.c
Normal file
367
libvips/convolution/spcor.c
Normal file
@ -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 <config.h>
|
||||
#endif /*HAVE_CONFIG_H*/
|
||||
#include <vips/intl.h>
|
||||
|
||||
#include <stdio.h>
|
||||
#include <math.h>
|
||||
|
||||
#include <vips/vips.h>
|
||||
|
||||
#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 );
|
||||
}
|
@ -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 \
|
||||
|
@ -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 )
|
||||
{
|
||||
|
@ -788,7 +788,6 @@ scanline_write( COLR *scanline, int width, FILE *fp )
|
||||
}
|
||||
|
||||
int i, j, beg, cnt;
|
||||
int c2;
|
||||
|
||||
if( width < MINELEN ||
|
||||
width > MAXELEN )
|
||||
|
@ -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*/
|
||||
|
@ -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*/
|
||||
|
Loading…
Reference in New Issue
Block a user