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

Merge branch 'master' into add-magick7

This commit is contained in:
John Cupitt 2016-08-14 10:25:03 +01:00
parent 8c05a2bbe8 933080a265
commit 9143bda915
29 changed files with 4178 additions and 3802 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.travis.yml
View File

@ -1,7 +1,7 @@
language: cpp
before_script:
- ./bootstrap.sh
- ./autogen.sh
- ./configure
--disable-dependency-tracking
--with-jpeg-includes=$JPEG/include

5
ChangeLog
View File

@ -29,9 +29,14 @@
- better upsizing with vips_resize()
- add imagemagick v7 support, thanks sachinwalia2k8
- added vips_worley(), vips_perlin() noise generators
- added vips_convf(), vips_convi(), vips_convasep(), vips_conva() ...
im_conv*() functions rewritten as classes
- vips_convsep() calls vips_convasep() for the approximate case
- new fixed-point vector path for convi is up to about 2x faster
- gif loader can write 1, 2, 3, or 4 bands depending on file contents
- support --strip for pngsave
- add svgz support [Felix Bünemann]
- rename boostrap.sh -> autogen.sh to help snapcraft
30/7/16 started 8.3.3
- fix performance regression in 8.3.2, thanks Lovell

2
README.md
View File

@ -80,7 +80,7 @@ and `gobject-introspection`, see the dependencies section below. For example:
Then build the build system with:
$ ./bootstrap.sh
$ ./autogen.sh
Debug build:

13
TODO
View File

@ -16,18 +16,7 @@
all with file and buffer source
- try:
$ vips avg broken.jpg[fail]
about 50% of the time it'll trigger a range of out-of-order reads and lock
for 10s or so while seq times out
- add more webp tests to py suite
- try moving some more of the CLI tests to py
- add APPROX convsep test?
- add more webp tests to py suite

0
bootstrap.sh → autogen.sh
View File

14
libvips/conversion/sequential.c
View File

@ -184,7 +184,7 @@ vips_sequential_generate( VipsRegion *or,
/* Exit the loop on timeout or condition passes. We have to
* be wary of spurious wakeups.
*/
while( r->top > sequential->y_pos )
while( r->top > sequential->y_pos ) {
#ifdef HAVE_COND_INIT
if( !g_cond_wait_until( sequential->ready,
sequential->lock, time ) )
@ -195,6 +195,14 @@ vips_sequential_generate( VipsRegion *or,
break;
#endif
/* We may have woken up because of an eval error.
*/
if( sequential->error ) {
g_mutex_unlock( sequential->lock );
return( -1 );
}
}
VIPS_GATE_STOP( "vips_sequential_generate: wait" );
VIPS_DEBUG_MSG_GREEN( "thread %p awake again ...\n",
@ -220,7 +228,7 @@ vips_sequential_generate( VipsRegion *or,
area.width = 1;
area.height = r->top - sequential->y_pos;
if( vips_region_prepare( ir, &area ) ) {
VIPS_DEBUG_MSG( "thread %p error, unlocking ...\n",
VIPS_DEBUG_MSG( "thread %p error, unlocking #1 ...\n",
g_thread_self() );
sequential->error = -1;
g_cond_broadcast( sequential->ready );
@ -237,7 +245,7 @@ vips_sequential_generate( VipsRegion *or,
VIPS_DEBUG_MSG_GREEN( "thread %p reading ...\n", g_thread_self() );
if( vips_region_prepare( ir, r ) ||
vips_region_region( or, ir, r, r->left, r->top ) ) {
VIPS_DEBUG_MSG( "thread %p error, unlocking ...\n",
VIPS_DEBUG_MSG( "thread %p error, unlocking #2 ...\n",
g_thread_self() );
sequential->error = -1;
g_cond_broadcast( sequential->ready );

10
libvips/convolution/Makefile.am
View File

@ -6,15 +6,15 @@ libconvolution_la_SOURCES = \
correlation.c \
correlation.h \
conv.c \
conva.c \
convf.c \
convi.c \
convasep.c \
convsep.c \
compass.c \
fastcor.c \
spcor.c \
sharpen.c \
gaussblur.c \
im_aconv.c \
im_aconvsep.c \
im_conv.c \
im_conv_f.c
gaussblur.c
AM_CPPFLAGS = -I${top_srcdir}/libvips/include @VIPS_CFLAGS@ @VIPS_INCLUDES@

73
libvips/convolution/conv.c
View File

@ -66,18 +66,15 @@ vips_conv_build( VipsObject *object )
VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( object );
VipsConvolution *convolution = (VipsConvolution *) object;
VipsConv *conv = (VipsConv *) object;
VipsImage **t = (VipsImage **)
vips_object_local_array( object, 4 );
VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
VipsImage *in;
INTMASK *imsk;
DOUBLEMASK *dmsk;
g_object_set( conv, "out", vips_image_new(), NULL );
if( VIPS_OBJECT_CLASS( vips_conv_parent_class )->build( object ) )
return( -1 );
g_object_set( conv, "out", vips_image_new(), NULL );
in = convolution->in;
/*
@ -85,13 +82,6 @@ vips_conv_build( VipsObject *object )
vips_matrixprint( convolution->M, NULL );
*/
if( !(imsk = im_vips2imask( convolution->M, class->nickname )) ||
!im_local_imask( convolution->out, imsk ) )
return( -1 );
if( !(dmsk = im_vips2mask( convolution->M, class->nickname )) ||
!im_local_dmask( convolution->out, dmsk ) )
return( -1 );
/* Unpack for processing.
*/
if( vips_image_decode( in, &t[0] ) )
@ -99,20 +89,30 @@ vips_conv_build( VipsObject *object )
in = t[0];
switch( conv->precision ) {
case VIPS_PRECISION_INTEGER:
if( im_conv( in, convolution->out, imsk ) )
case VIPS_PRECISION_FLOAT:
if( vips_convf( in, &t[1], convolution->M, NULL ) ||
vips_image_write( t[1], convolution->out ) )
return( -1 );
break;
case VIPS_PRECISION_FLOAT:
if( im_conv_f( in, convolution->out, dmsk ) )
case VIPS_PRECISION_INTEGER:
if( vips_convi( in, &t[1], convolution->M, NULL ) ||
vips_image_write( t[1], convolution->out ) )
return( -1 );
break;
case VIPS_PRECISION_APPROXIMATE:
{
DOUBLEMASK *dmsk;
if( !(dmsk = im_vips2mask( convolution->M, class->nickname )) ||
!im_local_dmask( convolution->out, dmsk ) )
return( -1 );
if( im_aconv( in, convolution->out, dmsk,
conv->layers, conv->cluster ) )
return( -1 );
}
break;
default:
@ -175,31 +175,40 @@ vips_conv_init( VipsConv *conv )
*
* Optional arguments:
*
* * @precision: calculation accuracy
* * @layers: number of layers for approximation
* * @cluster: cluster lines closer than this distance
* * @precision: #VipsPrecision, calculation accuracy
* * @layers: %gint, number of layers for approximation
* * @cluster: %gint, cluster lines closer than this distance
*
* Convolution.
*
* Perform a convolution of @in with @mask.
* Each output pixel is
* calculated as sigma[i]{pixel[i] * mask[i]} / scale + offset, where scale
* and offset are part of @mask.
* Each output pixel is calculated as:
*
* If @precision is #VIPS_PRECISION_INTEGER then the convolution is performed
* with integer arithmetic and the output image
* |[
* sigma[i]{pixel[i] * mask[i]} / scale + offset
* ]|
*
* where scale and offset are part of @mask.
*
* If @precision is #VIPS_PRECISION_INTEGER, then
* elements of @mask are converted to
* integers before convolution, using rint(),
* and the output image
* always has the same #VipsBandFormat as the input image.
*
* Convolutions on unsigned 8-bit images are calculated with the
* processor's vector unit, if possible. Disable this with --vips-novector or
* IM_NOVECTOR.
* For #VIPS_FORMAT_UCHAR images, vips_conv() uses a fast vector path based on
* fixed-point arithmetic. This can produce slightly different results.
* Disable the vector path with `--vips-novector` or `VIPS_NOVECTOR` or
* vips_vector_set_enabled().
*
* If @precision is #VIPS_PRECISION_FLOAT then the convolution is performed
* with floating-point arithmetic. The output image
* is always %VIPS_FORMAT_FLOAT unless @in is %VIPS_FORMAT_DOUBLE, in which case
* @out is also %VIPS_FORMAT_DOUBLE.
* is always #VIPS_FORMAT_FLOAT unless @in is #VIPS_FORMAT_DOUBLE, in which case
* @out is also #VIPS_FORMAT_DOUBLE.
*
* If @precision is #VIPS_PRECISION_APPROXIMATE then the output image
* If @precision is #VIPS_PRECISION_APPROXIMATE then, like
* #VIPS_PRECISION_INTEGER, @mask is converted to int before convolution, and
* the output image
* always has the same #VipsBandFormat as the input image.
*
* Larger values for @layers give more accurate
@ -211,6 +220,8 @@ vips_conv_init( VipsConv *conv )
* Smaller values of @cluster will give more accurate results, but be slower
* and use more memory. 10% of the mask radius is a good rule of thumb.
*
* See also: vips_convsep().
*
* Returns: 0 on success, -1 on error
*/
int

1360
libvips/convolution/conva.c Normal file
View File

File diff suppressed because it is too large Load Diff

958
libvips/convolution/convasep.c Normal file
View File

@ -0,0 +1,958 @@
/* convasep ... separable approximate convolution
*
* This operation does an approximate, seperable convolution.
*
* Author: John Cupitt & Nicolas Robidoux
* Written on: 31/5/11
* Modified on:
* 31/5/11
* - from im_conv()
* 5/7/16
* - 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
*/
/*
See:
http://incubator.quasimondo.com/processing/stackblur.pde
This thing is a little like stackblur, but generalised to any separable
mask.
*/
/*
TODO
- how about making a cumulative image and then subtracting points in
that, rather than keeping a set of running totals
faster?
we could then use orc to write a bit of code to implement this set
of lines
stackoverflow has an algorithm for cumulativization using SIMD and
threads, see that font rasterization with rust piece on medium by
ralph levien
*/
/*
#define DEBUG
#define VIPS_DEBUG
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif /*HAVE_CONFIG_H*/
#include <vips/intl.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <math.h>
#include <vips/vips.h>
#include <vips/vector.h>
#include <vips/debug.h>
#include <vips/internal.h>
#include "pconvolution.h"
/* Maximum number of lines we can break the mask into.
*/
#define MAX_LINES (1000)
/* Euclid's algorithm. Use this to common up mults.
*/
static int
gcd( int a, int b )
{
if( b == 0 )
return( abs( a ) );
else
return( gcd( b, a % b ) );
}
typedef struct {
VipsConvolution parent_instance;
int layers;
int divisor;
int rounding;
int offset;
/* The "width" of the mask, ie. n for our 1xn or nx1 argument, plus
* an int version of our mask.
*/
int width;
VipsImage *iM;
/* The mask broken into a set of lines.
*
* Start is the left-most pixel in the line, end is one beyond the
* right-most pixel.
*/
int n_lines;
int start[MAX_LINES];
int end[MAX_LINES];
int factor[MAX_LINES];
} VipsConvasep;
typedef VipsConvolutionClass VipsConvasepClass;
G_DEFINE_TYPE( VipsConvasep, vips_convasep, VIPS_TYPE_CONVOLUTION );
static void
vips_convasep_line_start( VipsConvasep *convasep, int x, int factor )
{
convasep->start[convasep->n_lines] = x;
convasep->factor[convasep->n_lines] = factor;
}
static int
vips_convasep_line_end( VipsConvasep *convasep, int x )
{
VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( convasep );
convasep->end[convasep->n_lines] = x;
if( convasep->n_lines >= MAX_LINES - 1 ) {
vips_error( class->nickname, "%s", _( "mask too complex" ) );
return( -1 );
}
convasep->n_lines += 1;
return( 0 );
}
/* Break a mask into lines.
*/
static int
vips_convasep_decompose( VipsConvasep *convasep )
{
VipsImage *iM = convasep->iM;
double *coeff = (double *) VIPS_IMAGE_ADDR( iM, 0, 0 );
double scale = vips_image_get_scale( iM );
double offset = vips_image_get_offset( iM );
double max;
double min;
double depth;
double sum;
double area;
int layers;
int layers_above;
int layers_below;
int z, n, x;
VIPS_DEBUG_MSG( "vips_convasep_decompose: "
"breaking into %d layers ...\n", convasep->layers );
/* Find mask range. We must always include the zero axis in the mask.
*/
max = 0;
min = 0;
for( x = 0; x < convasep->width; x++ ) {
if( coeff[x] > max )
max = coeff[x];
if( coeff[x] < min )
min = coeff[x];
}
/* The zero axis must fall on a layer boundary. Estimate the
* depth, find n-lines-above-zero, get exact depth, then calculate a
* fixed n-lines which includes any negative parts.
*/
depth = (max - min) / convasep->layers;
layers_above = ceil( max / depth );
depth = max / layers_above;
layers_below = floor( min / depth );
layers = layers_above - layers_below;
VIPS_DEBUG_MSG( "depth = %g, layers = %d\n", depth, layers );
/* For each layer, generate a set of lines which are inside the
* perimeter. Work down from the top.
*/
for( z = 0; z < layers; z++ ) {
double y = max - (1 + z) * depth;
/* y plus half depth ... ie. the layer midpoint.
*/
double y_ph = y + depth / 2;
/* Odd, but we must avoid rounding errors that make us miss 0
* in the line above.
*/
int y_positive = z < layers_above;
int inside;
/* Start outside the perimeter.
*/
inside = 0;
for( x = 0; x < convasep->width; x++ ) {
/* The vertical line from mask[z] to 0 is inside. Is
* our current square (x, y) part of that line?
*/
if( (y_positive && coeff[x] >= y_ph) ||
(!y_positive && coeff[x] <= y_ph) ) {
if( !inside ) {
vips_convasep_line_start( convasep, x,
y_positive ? 1 : -1 );
inside = 1;
}
}
else if( inside ) {
if( vips_convasep_line_end( convasep, x ) )
return( -1 );
inside = 0;
}
}
if( inside &&
vips_convasep_line_end( convasep, convasep->width ) )
return( -1 );
}
/* Can we common up any lines? Search for lines with identical
* start/end.
*/
for( z = 0; z < convasep->n_lines; z++ ) {
for( n = z + 1; n < convasep->n_lines; n++ ) {
if( convasep->start[z] == convasep->start[n] &&
convasep->end[z] == convasep->end[n] ) {
convasep->factor[z] += convasep->factor[n];
/* n can be deleted. Do this in a separate
* pass below.
*/
convasep->factor[n] = 0;
}
}
}
/* Now we can remove all factor 0 lines.
*/
for( z = 0; z < convasep->n_lines; z++ ) {
if( convasep->factor[z] == 0 ) {
for( x = z; x < convasep->n_lines; x++ ) {
convasep->start[x] = convasep->start[x + 1];
convasep->end[x] = convasep->end[x + 1];
convasep->factor[x] = convasep->factor[x + 1];
}
convasep->n_lines -= 1;
}
}
/* Find the area of the lines.
*/
area = 0;
for( z = 0; z < convasep->n_lines; z++ )
area += convasep->factor[z] *
(convasep->end[z] - convasep->start[z]);
/* Strength reduction: if all lines are divisible by n, we can move
* that n out into the ->area factor. The aim is to produce as many
* factor 1 lines as we can and to reduce the chance of overflow.
*/
x = convasep->factor[0];
for( z = 1; z < convasep->n_lines; z++ )
x = gcd( x, convasep->factor[z] );
for( z = 0; z < convasep->n_lines; z++ )
convasep->factor[z] /= x;
area *= x;
/* Find the area of the original mask.
*/
sum = 0;
for( z = 0; z < convasep->width; z++ )
sum += coeff[z];
convasep->divisor = VIPS_RINT( sum * area / scale );
if( convasep->divisor == 0 )
convasep->divisor = 1;
convasep->rounding = (convasep->divisor + 1) / 2;
convasep->offset = offset;
#ifdef DEBUG
/* ASCII-art layer drawing.
*/
printf( "lines:\n" );
for( z = 0; z < convasep->n_lines; z++ ) {
printf( "%3d - %2d x ", z, convasep->factor[z] );
for( x = 0; x < 55; x++ ) {
int rx = x * (convasep->width + 1) / 55;
if( rx >= convasep->start[z] && rx < convasep->end[z] )
printf( "#" );
else
printf( " " );
}
printf( " %3d .. %3d\n", convasep->start[z], convasep->end[z] );
}
printf( "divisor = %d\n", convasep->divisor );
printf( "rounding = %d\n", convasep->rounding );
printf( "offset = %d\n", convasep->offset );
#endif /*DEBUG*/
return( 0 );
}
/* Our sequence value.
*/
typedef struct {
VipsConvasep *convasep;
VipsRegion *ir; /* Input region */
int *start; /* Offsets for start and stop */
int *end;
/* The sums for each line. int for integer types, double for floating
* point types.
*/
int *isum;
double *dsum;
int last_stride; /* Avoid recalcing offsets, if we can */
} VipsConvasepSeq;
/* Free a sequence value.
*/
static int
vips_convasep_stop( void *vseq, void *a, void *b )
{
VipsConvasepSeq *seq = (VipsConvasepSeq *) vseq;
VIPS_UNREF( seq->ir );
VIPS_FREE( seq->start );
VIPS_FREE( seq->end );
VIPS_FREE( seq->isum );
VIPS_FREE( seq->dsum );
return( 0 );
}
/* Convolution start function.
*/
static void *
vips_convasep_start( VipsImage *out, void *a, void *b )
{
VipsImage *in = (IMAGE *) a;
VipsConvasep *convasep = (VipsConvasep *) b;
VipsConvasepSeq *seq;
if( !(seq = VIPS_NEW( out, VipsConvasepSeq )) )
return( NULL );
/* Init!
*/
seq->convasep = convasep;
seq->ir = vips_region_new( in );
seq->start = VIPS_ARRAY( NULL, convasep->n_lines, int );
seq->end = VIPS_ARRAY( NULL, convasep->n_lines, int );
seq->isum = NULL;
seq->dsum = NULL;
if( vips_band_format_isint( out->BandFmt ) )
seq->isum = VIPS_ARRAY( NULL, convasep->n_lines, int );
else
seq->dsum = VIPS_ARRAY( NULL, convasep->n_lines, double );
seq->last_stride = -1;
if( !seq->ir ||
!seq->start ||
!seq->end ||
(!seq->isum && !seq->dsum) ) {
vips_convasep_stop( seq, in, convasep );
return( NULL );
}
return( seq );
}
#define CLIP_UCHAR( V ) \
G_STMT_START { \
if( (V) < 0 ) \
(V) = 0; \
else if( (V) > UCHAR_MAX ) \
(V) = UCHAR_MAX; \
} G_STMT_END
#define CLIP_CHAR( V ) \
G_STMT_START { \
if( (V) < SCHAR_MIN ) \
(V) = SCHAR_MIN; \
else if( (V) > SCHAR_MAX ) \
(V) = SCHAR_MAX; \
} G_STMT_END
#define CLIP_USHORT( V ) \
G_STMT_START { \
if( (V) < 0 ) \
(V) = 0; \
else if( (V) > USHRT_MAX ) \
(V) = USHRT_MAX; \
} G_STMT_END
#define CLIP_SHORT( V ) \
G_STMT_START { \
if( (V) < SHRT_MIN ) \
(V) = SHRT_MIN; \
else if( (V) > SHRT_MAX ) \
(V) = SHRT_MAX; \
} G_STMT_END
#define CLIP_NONE( V ) {}
/* The h and v loops are very similar, but also annoyingly different. Keep
* them separate for easy debugging.
*/
#define HCONV_INT( TYPE, CLIP ) { \
for( i = 0; i < bands; i++ ) { \
int *isum = seq->isum; \
\
TYPE *q; \
TYPE *p; \
int sum; \
\
p = i + (TYPE *) VIPS_REGION_ADDR( ir, r->left, r->top + y ); \
q = i + (TYPE *) VIPS_REGION_ADDR( or, r->left, r->top + y ); \
\
sum = 0; \
for( z = 0; z < n_lines; z++ ) { \
isum[z] = 0; \
for( x = seq->start[z]; x < seq->end[z]; x += istride ) \
isum[z] += p[x]; \
sum += convasep->factor[z] * isum[z]; \
} \
\
/* Don't add offset ... we only want to do that once, do it on \
* the vertical pass. \
*/ \
sum = (sum + convasep->rounding) / convasep->divisor; \
CLIP( sum ); \
*q = sum; \
q += ostride; \
\
for( x = 1; x < r->width; x++ ) { \
sum = 0; \
for( z = 0; z < n_lines; z++ ) { \
isum[z] += p[seq->end[z]]; \
isum[z] -= p[seq->start[z]]; \
sum += convasep->factor[z] * isum[z]; \
} \
p += istride; \
sum = (sum + convasep->rounding) / convasep->divisor; \
CLIP( sum ); \
*q = sum; \
q += ostride; \
} \
} \
}
#define HCONV_FLOAT( TYPE ) { \
for( i = 0; i < bands; i++ ) { \
double *dsum = seq->dsum; \
\
TYPE *q; \
TYPE *p; \
double sum; \
\
p = i + (TYPE *) VIPS_REGION_ADDR( ir, r->left, r->top + y ); \
q = i + (TYPE *) VIPS_REGION_ADDR( or, r->left, r->top + y ); \
\
sum = 0; \
for( z = 0; z < n_lines; z++ ) { \
dsum[z] = 0; \
for( x = seq->start[z]; x < seq->end[z]; x += istride ) \
dsum[z] += p[x]; \
sum += convasep->factor[z] * dsum[z]; \
} \
\
/* Don't add offset ... we only want to do that once, do it on \
* the vertical pass. \
*/ \
sum = sum / convasep->divisor; \
*q = sum; \
q += ostride; \
\
for( x = 1; x < r->width; x++ ) { \
sum = 0; \
for( z = 0; z < n_lines; z++ ) { \
dsum[z] += p[seq->end[z]]; \
dsum[z] -= p[seq->start[z]]; \
sum += convasep->factor[z] * dsum[z]; \
} \
p += istride; \
sum = sum / convasep->divisor; \
*q = sum; \
q += ostride; \
} \
} \
}
/* Do horizontal masks ... we scan the mask along scanlines.
*/
static int
vips_convasep_generate_horizontal( VipsRegion *or,
void *vseq, void *a, void *b, gboolean *stop )
{
VipsConvasepSeq *seq = (VipsConvasepSeq *) vseq;
VipsImage *in = (VipsImage *) a;
VipsConvasep *convasep = (VipsConvasep *) b;
VipsConvolution *convolution = (VipsConvolution *) convasep;
VipsRegion *ir = seq->ir;
const int n_lines = convasep->n_lines;
VipsRect *r = &or->valid;
/* Double the bands (notionally) for complex.
*/
int bands = vips_band_format_iscomplex( in->BandFmt ) ?
2 * in->Bands : in->Bands;
VipsRect s;
int x, y, z, i;
int istride;
int ostride;
/* Prepare the section of the input image we need. A little larger
* than the section of the output image we are producing.
*/
s = *r;
s.width += convasep->width - 1;
if( vips_region_prepare( ir, &s ) )
return( -1 );
/* Stride can be different for the vertical case, keep this here for
* ease of direction change.
*/
istride = VIPS_IMAGE_SIZEOF_PEL( in ) /
VIPS_IMAGE_SIZEOF_ELEMENT( in );
ostride = VIPS_IMAGE_SIZEOF_PEL( convolution->out ) /
VIPS_IMAGE_SIZEOF_ELEMENT( convolution->out );
/* Init offset array.
*/
if( seq->last_stride != istride ) {
seq->last_stride = istride;
for( z = 0; z < n_lines; z++ ) {
seq->start[z] = convasep->start[z] * istride;
seq->end[z] = convasep->end[z] * istride;
}
}
for( y = 0; y < r->height; y++ ) {
switch( in->BandFmt ) {
case VIPS_FORMAT_UCHAR:
HCONV_INT( unsigned char, CLIP_UCHAR );
break;
case VIPS_FORMAT_CHAR:
HCONV_INT( signed char, CLIP_CHAR );
break;
case VIPS_FORMAT_USHORT:
HCONV_INT( unsigned short, CLIP_USHORT );
break;
case VIPS_FORMAT_SHORT:
HCONV_INT( signed short, CLIP_SHORT );
break;
case VIPS_FORMAT_UINT:
HCONV_INT( unsigned int, CLIP_NONE );
break;
case VIPS_FORMAT_INT:
HCONV_INT( signed int, CLIP_NONE );
break;
case VIPS_FORMAT_FLOAT:
case VIPS_FORMAT_COMPLEX:
HCONV_FLOAT( float );
break;
case VIPS_FORMAT_DOUBLE:
case VIPS_FORMAT_DPCOMPLEX:
HCONV_FLOAT( double );
break;
default:
g_assert_not_reached();
}
}
return( 0 );
}
#define VCONV_INT( TYPE, CLIP ) { \
for( x = 0; x < sz; x++ ) { \
int *isum = seq->isum; \
\
TYPE *q; \
TYPE *p; \
int sum; \
\
p = x + (TYPE *) VIPS_REGION_ADDR( ir, r->left, r->top ); \
q = x + (TYPE *) VIPS_REGION_ADDR( or, r->left, r->top ); \
\
sum = 0; \
for( z = 0; z < n_lines; z++ ) { \
isum[z] = 0; \
for( y = seq->start[z]; y < seq->end[z]; y += istride ) \
isum[z] += p[y]; \
sum += convasep->factor[z] * isum[z]; \
} \
sum = (sum + convasep->rounding) / convasep->divisor + \
convasep->offset; \
CLIP( sum ); \
*q = sum; \
q += ostride; \
\
for( y = 1; y < r->height; y++ ) { \
sum = 0; \
for( z = 0; z < n_lines; z++ ) { \
isum[z] += p[seq->end[z]]; \
isum[z] -= p[seq->start[z]]; \
sum += convasep->factor[z] * isum[z]; \
} \
p += istride; \
sum = (sum + convasep->rounding) / convasep->divisor + \
convasep->offset; \
CLIP( sum ); \
*q = sum; \
q += ostride; \
} \
} \
}
#define VCONV_FLOAT( TYPE ) { \
for( x = 0; x < sz; x++ ) { \
double *dsum = seq->dsum; \
\
TYPE *q; \
TYPE *p; \
double sum; \
\
p = x + (TYPE *) VIPS_REGION_ADDR( ir, r->left, r->top ); \
q = x + (TYPE *) VIPS_REGION_ADDR( or, r->left, r->top ); \
\
sum = 0; \
for( z = 0; z < n_lines; z++ ) { \
dsum[z] = 0; \
for( y = seq->start[z]; y < seq->end[z]; y += istride ) \
dsum[z] += p[y]; \
sum += convasep->factor[z] * dsum[z]; \
} \
sum = sum / convasep->divisor + convasep->offset; \
*q = sum; \
q += ostride; \
\
for( y = 1; y < r->height; y++ ) { \
sum = 0; \
for( z = 0; z < n_lines; z++ ) { \
dsum[z] += p[seq->end[z]]; \
dsum[z] -= p[seq->start[z]]; \
sum += convasep->factor[z] * dsum[z]; \
} \
p += istride; \
sum = sum / convasep->divisor + convasep->offset; \
*q = sum; \
q += ostride; \
} \
} \
}
/* Do vertical masks ... we scan the mask down columns of pixels. Copy-paste
* from above with small changes.
*/
static int
vips_convasep_generate_vertical( VipsRegion *or,
void *vseq, void *a, void *b, gboolean *stop )
{
VipsConvasepSeq *seq = (VipsConvasepSeq *) vseq;
VipsImage *in = (VipsImage *) a;
VipsConvasep *convasep = (VipsConvasep *) b;
VipsConvolution *convolution = (VipsConvolution *) convasep;
VipsRegion *ir = seq->ir;
const int n_lines = convasep->n_lines;
VipsRect *r = &or->valid;
/* Double the width (notionally) for complex.
*/
int sz = vips_band_format_iscomplex( in->BandFmt ) ?
2 * VIPS_REGION_N_ELEMENTS( or ) : VIPS_REGION_N_ELEMENTS( or );
VipsRect s;
int x, y, z;
int istride;
int ostride;
/* Prepare the section of the input image we need. A little larger
* than the section of the output image we are producing.
*/
s = *r;
s.height += convasep->width - 1;
if( vips_region_prepare( ir, &s ) )
return( -1 );
/* Stride can be different for the vertical case, keep this here for
* ease of direction change.
*/
istride = VIPS_REGION_LSKIP( ir ) / VIPS_IMAGE_SIZEOF_ELEMENT( in );
ostride = VIPS_REGION_LSKIP( or ) /
VIPS_IMAGE_SIZEOF_ELEMENT( convolution->out );
/* Init offset array.
*/
if( seq->last_stride != istride ) {
seq->last_stride = istride;
for( z = 0; z < n_lines; z++ ) {
seq->start[z] = convasep->start[z] * istride;
seq->end[z] = convasep->end[z] * istride;
}
}
switch( in->BandFmt ) {
case VIPS_FORMAT_UCHAR:
VCONV_INT( unsigned char, CLIP_UCHAR );
break;
case VIPS_FORMAT_CHAR:
VCONV_INT( signed char, CLIP_CHAR );
break;
case VIPS_FORMAT_USHORT:
VCONV_INT( unsigned short, CLIP_USHORT );
break;
case VIPS_FORMAT_SHORT:
VCONV_INT( signed short, CLIP_SHORT );
break;
case VIPS_FORMAT_UINT:
VCONV_INT( unsigned int, CLIP_NONE );
break;
case VIPS_FORMAT_INT:
VCONV_INT( signed int, CLIP_NONE );
break;
case VIPS_FORMAT_FLOAT:
case VIPS_FORMAT_COMPLEX:
VCONV_FLOAT( float );
break;
case VIPS_FORMAT_DOUBLE:
case VIPS_FORMAT_DPCOMPLEX:
VCONV_FLOAT( double );
break;
default:
g_assert_not_reached();
}
return( 0 );
}
static int
vips_convasep_pass( VipsConvasep *convasep,
VipsImage *in, VipsImage **out, VipsDirection direction )
{
VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( convasep );
VipsGenerateFn gen;
*out = vips_image_new();
if( vips_image_pipelinev( *out,
VIPS_DEMAND_STYLE_SMALLTILE, in, NULL ) )
return( -1 );
if( direction == VIPS_DIRECTION_HORIZONTAL ) {
(*out)->Xsize -= convasep->width - 1;
gen = vips_convasep_generate_horizontal;
}
else {
(*out)->Ysize -= convasep->width - 1;
gen = vips_convasep_generate_vertical;
}
if( (*out)->Xsize <= 0 ||
(*out)->Ysize <= 0 ) {
vips_error( class->nickname,
"%s", _( "image too small for mask" ) );
return( -1 );
}
if( vips_image_generate( *out,
vips_convasep_start, gen, vips_convasep_stop, in, convasep ) )
return( -1 );
return( 0 );
}
static int
vips_convasep_build( VipsObject *object )
{
VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( object );
VipsConvolution *convolution = (VipsConvolution *) object;
VipsConvasep *convasep = (VipsConvasep *) object;
VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
VipsImage *in;
if( VIPS_OBJECT_CLASS( vips_convasep_parent_class )->build( object ) )
return( -1 );
if( vips_check_separable( class->nickname, convolution->M ) )
return( -1 );
/* An int version of our mask.
*/
if( vips__image_intize( convolution->M, &t[3] ) )
return( -1 );
convasep->iM = t[3];
convasep->width = convasep->iM->Xsize * convasep->iM->Ysize;
in = convolution->in;
if( vips_convasep_decompose( convasep ) )
return( -1 );
g_object_set( convasep, "out", vips_image_new(), NULL );
if(
vips_embed( in, &t[0],
convasep->width / 2,
convasep->width / 2,
in->Xsize + convasep->width - 1,
in->Ysize + convasep->width - 1,
"extend", VIPS_EXTEND_COPY,
NULL ) ||
vips_convasep_pass( convasep,
t[0], &t[1], VIPS_DIRECTION_HORIZONTAL ) ||
vips_convasep_pass( convasep,
t[1], &t[2], VIPS_DIRECTION_VERTICAL ) ||
vips_image_write( t[2], convolution->out ) )
return( -1 );
convolution->out->Xoffset = 0;
convolution->out->Yoffset = 0;
return( 0 );
}
static void
vips_convasep_class_init( VipsConvasepClass *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 = "convasep";
object_class->description =
_( "approximate separable integer convolution" );
object_class->build = vips_convasep_build;
VIPS_ARG_INT( class, "layers", 104,
_( "Layers" ),
_( "Use this many layers in approximation" ),
VIPS_ARGUMENT_OPTIONAL_INPUT,
G_STRUCT_OFFSET( VipsConvasep, layers ),
1, 1000, 5 );
}
static void
vips_convasep_init( VipsConvasep *convasep )
{
convasep->layers = 5;
convasep->n_lines = 0;
}
/**
* vips_convasep:
* @in: input image
* @out: output image
* @mask: convolve with this mask
* @...: %NULL-terminated list of optional named arguments
*
* Optional arguments:
*
* * @layers: %gint, number of layers for approximation
*
* Approximate separable integer convolution. This is a low-level operation, see
* vips_convsep() for something more convenient.
*
* The image is convolved twice: once with @mask and then again with @mask
* rotated by 90 degrees.
* @mask must be 1xn or nx1 elements.
* Elements of @mask are converted to
* integers before convolution.
*
* Larger values for @layers give more accurate
* results, but are slower. As @layers approaches the mask radius, the
* accuracy will become close to exact convolution and the speed will drop to
* match. For many large masks, such as Gaussian, @layers need be only 10% of
* this value and accuracy will still be good.
*
* The output image
* always has the same #VipsBandFormat as the input image.
*
* See also: vips_convsep().
*
* Returns: 0 on success, -1 on error
*/
int
vips_convasep( VipsImage *in, VipsImage **out, VipsImage *mask, ... )
{
va_list ap;
int result;
va_start( ap, mask );
result = vips_call_split( "convasep", ap, in, out, mask );
va_end( ap );
return( result );
}

406
libvips/convolution/convf.c Normal file
View File

@ -0,0 +1,406 @@
/* convf
*
* Copyright: 1990, N. Dessipris.
*
* Author: Nicos Dessipris & Kirk Martinez
* Written on: 29/04/1991
* Modified on: 19/05/1991
* 8/7/93 JC
* - adapted for partial v2
* - memory leaks fixed
* - ANSIfied
* 12/7/93 JC
* - adapted im_convbi() to im_convbf()
* 7/10/94 JC
* - new IM_ARRAY() macro
* - evalend callbacks
* - more typedef
* 9/3/01 JC
* - redone from im_conv()
* 27/7/01 JC
* - rejects masks with scale == 0
* 7/4/04
* - now uses im_embed() with edge stretching on the input, not
* the output
* - sets Xoffset / Yoffset
* 11/11/05
* - simpler inner loop avoids gcc4 bug
* 12/11/09
* - only rebuild the buffer offsets if bpl changes
* - tiny speedups and cleanups
* - add restrict, though it doesn't seem to help gcc
* - add mask-all-zero check
* 13/11/09
* - rename as im_conv_f() to make it easier for vips.c to make the
* overloaded version
* 3/2/10
* - gtkdoc
* - more cleanups
* 1/10/10
* - support complex (just double the bands)
* 29/10/10
* - get rid of im_convsep_f(), just call this twice, no longer worth
* keeping two versions
* 15/10/11 Nicolas
* - handle offset correctly in seperable convolutions
* 26/1/16 Lovell Fuller
* - remove Duff for a 25% speedup
* 23/6/16
* - 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 <stdlib.h>
#include <limits.h>
#include <vips/vips.h>
#include "pconvolution.h"
typedef struct {
VipsConvolution parent_instance;
/* We make a smaller version of the mask with the zeros squeezed out.
*/
int nnz; /* Number of non-zero mask elements */
double *coeff; /* Array of non-zero mask coefficients */
int *coeff_pos; /* Index of each nnz element in mask->coeff */
} VipsConvf;
typedef VipsConvolutionClass VipsConvfClass;
G_DEFINE_TYPE( VipsConvf, vips_convf, VIPS_TYPE_CONVOLUTION );
/* Our sequence value.
*/
typedef struct {
VipsConvf *convf;
VipsRegion *ir; /* Input region */
int *offsets; /* Offsets for each non-zero matrix element */
VipsPel **pts; /* Per-non-zero mask element image pointers */
int last_bpl; /* Avoid recalcing offsets, if we can */
} VipsConvfSequence;
/* Free a sequence value.
*/
static int
vips_convf_stop( void *vseq, void *a, void *b )
{
VipsConvfSequence *seq = (VipsConvfSequence *) vseq;
VIPS_UNREF( seq->ir );
return( 0 );
}
/* Convolution start function.
*/
static void *
vips_convf_start( VipsImage *out, void *a, void *b )
{
VipsImage *in = (VipsImage *) a;
VipsConvf *convf = (VipsConvf *) b;
VipsConvfSequence *seq;
if( !(seq = VIPS_NEW( out, VipsConvfSequence )) )
return( NULL );
seq->convf = convf;
seq->ir = NULL;
seq->pts = NULL;
seq->last_bpl = -1;
seq->ir = vips_region_new( in );
if( !(seq->offsets = VIPS_ARRAY( out, convf->nnz, int )) ||
!(seq->pts = VIPS_ARRAY( out, convf->nnz, VipsPel * )) ) {
vips_convf_stop( seq, in, convf );
return( NULL );
}
return( (void *) seq );
}
#define CONV_FLOAT( ITYPE, OTYPE ) { \
ITYPE ** restrict p = (ITYPE **) seq->pts; \
OTYPE * restrict q = (OTYPE *) VIPS_REGION_ADDR( or, le, y ); \
\
for( x = 0; x < sz; x++ ) { \
double sum; \
int i; \
\
sum = 0; \
for ( i = 0; i < nnz; i++ ) \
sum += t[i] * p[i][x]; \
\
sum = (sum / scale) + offset; \
\
q[x] = sum; \
} \
}
/* Convolve!
*/
static int
vips_convf_gen( REGION *or, void *vseq, void *a, void *b, gboolean *stop )
{
VipsConvfSequence *seq = (VipsConvfSequence *) vseq;
VipsConvf *convf = (VipsConvf *) b;
VipsConvolution *convolution = (VipsConvolution *) convf;
VipsImage *M = convolution->M;
double scale = vips_image_get_scale( M );
double offset = vips_image_get_offset( M );
VipsImage *in = (VipsImage *) a;
VipsRegion *ir = seq->ir;
double * restrict t = convf->coeff;
const int nnz = convf->nnz;
VipsRect *r = &or->valid;
int le = r->left;
int to = r->top;
int bo = VIPS_RECT_BOTTOM( r );
int sz = VIPS_REGION_N_ELEMENTS( or ) *
(vips_band_format_iscomplex( in->BandFmt ) ? 2 : 1);
VipsRect s;
int x, y, z, i;
/* Prepare the section of the input image we need. A little larger
* than the section of the output image we are producing.
*/
s = *r;
s.width += M->Xsize - 1;
s.height += M->Ysize - 1;
if( vips_region_prepare( ir, &s ) )
return( -1 );
/* Fill offset array. Only do this if the bpl has changed since the
* previous vips_region_prepare().
*/
if( seq->last_bpl != VIPS_REGION_LSKIP( ir ) ) {
seq->last_bpl = VIPS_REGION_LSKIP( ir );
for( i = 0; i < nnz; i++ ) {
z = convf->coeff_pos[i];
x = z % M->Xsize;
y = z / M->Xsize;
seq->offsets[i] =
VIPS_REGION_ADDR( ir, x + le, y + to ) -
VIPS_REGION_ADDR( ir, le, to );
}
}
for( y = to; y < bo; y++ ) {
/* Init pts for this line of PELs.
*/
for( z = 0; z < nnz; z++ )
seq->pts[z] = seq->offsets[z] +
VIPS_REGION_ADDR( ir, le, y );
switch( in->BandFmt ) {
case VIPS_FORMAT_UCHAR:
CONV_FLOAT( unsigned char, float );
break;
case VIPS_FORMAT_CHAR:
CONV_FLOAT( signed char, float );
break;
case VIPS_FORMAT_USHORT:
CONV_FLOAT( unsigned short, float );
break;
case VIPS_FORMAT_SHORT:
CONV_FLOAT( signed short, float );
break;
case VIPS_FORMAT_UINT:
CONV_FLOAT( unsigned int, float );
break;
case VIPS_FORMAT_INT:
CONV_FLOAT( signed int, float );
break;
case VIPS_FORMAT_FLOAT:
case VIPS_FORMAT_COMPLEX:
CONV_FLOAT( float, float );
break;
case VIPS_FORMAT_DOUBLE:
case VIPS_FORMAT_DPCOMPLEX:
CONV_FLOAT( double, double );
break;
default:
g_assert_not_reached();
}
}
return( 0 );
}
static int
vips_convf_build( VipsObject *object )
{
VipsConvolution *convolution = (VipsConvolution *) object;
VipsConvf *convf = (VipsConvf *) object;
VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
VipsImage *in;
VipsImage *M;
double *coeff;
int ne;
int i;
if( VIPS_OBJECT_CLASS( vips_convf_parent_class )->build( object ) )
return( -1 );
M = convolution->M;
coeff = (double *) VIPS_IMAGE_ADDR( M, 0, 0 );
ne = M->Xsize * M->Ysize;
if( !(convf->coeff = VIPS_ARRAY( object, ne, double )) ||
!(convf->coeff_pos = VIPS_ARRAY( object, ne, int )) )
return( -1 );
/* Find non-zero mask elements.
*/
for( i = 0; i < ne; i++ )
if( coeff[i] ) {
convf->coeff[convf->nnz] = coeff[i];
convf->coeff_pos[convf->nnz] = i;
convf->nnz += 1;
}
/* Was the whole mask zero? We must have at least 1 element in there:
* set it to zero.
*/
if( convf->nnz == 0 ) {
convf->coeff[0] = 0;
convf->coeff_pos[0] = 0;
convf->nnz = 1;
}
in = convolution->in;
if( vips_embed( in, &t[0],
M->Xsize / 2, M->Ysize / 2,
in->Xsize + M->Xsize - 1, in->Ysize + M->Ysize - 1,
"extend", VIPS_EXTEND_COPY,
NULL ) )
return( -1 );
in = t[0];
g_object_set( convf, "out", vips_image_new(), NULL );
if( vips_image_pipelinev( convolution->out,
VIPS_DEMAND_STYLE_SMALLTILE, in, NULL ) )
return( -1 );
convolution->out->Xoffset = 0;
convolution->out->Yoffset = 0;
/* Prepare output. Consider a 7x7 mask and a 7x7 image --- the output
* would be 1x1.
*/
if( vips_bandfmt_isint( in->BandFmt ) )
convolution->out->BandFmt = IM_BANDFMT_FLOAT;
convolution->out->Xsize -= M->Xsize - 1;
convolution->out->Ysize -= M->Ysize - 1;
if( vips_image_generate( convolution->out,
vips_convf_start, vips_convf_gen, vips_convf_stop, in, convf ) )
return( -1 );
convolution->out->Xoffset = -M->Xsize / 2;
convolution->out->Yoffset = -M->Ysize / 2;
return( 0 );
}
static void
vips_convf_class_init( VipsConvfClass *class )
{
VipsObjectClass *object_class = (VipsObjectClass *) class;
object_class->nickname = "convf";
object_class->description = _( "float convolution operation" );
object_class->build = vips_convf_build;
}
static void
vips_convf_init( VipsConvf *convf )
{
convf->nnz = 0;
convf->coeff = NULL;
convf->coeff_pos = NULL;
}
/**
* vips_convf:
* @in: input image
* @out: output image
* @mask: convolve with this mask
* @...: %NULL-terminated list of optional named arguments
*
* Convolution. This is a low-level operation, see vips_conv() for something
* more convenient.
*
* Perform a convolution of @in with @mask.
* Each output pixel is
* calculated as sigma[i]{pixel[i] * mask[i]} / scale + offset, where scale
* and offset are part of @mask.
*
* The convolution is performed with floating-point arithmetic. The output image
* is always #VIPS_FORMAT_FLOAT unless @in is #VIPS_FORMAT_DOUBLE, in which case
* @out is also #VIPS_FORMAT_DOUBLE.
*
* See also: vips_conv().
*
* Returns: 0 on success, -1 on error
*/
int
vips_convf( VipsImage *in, VipsImage **out, VipsImage *mask, ... )
{
va_list ap;
int result;
va_start( ap, mask );
result = vips_call_split( "convf", ap, in, out, mask );
va_end( ap );
return( result );
}

1066
libvips/convolution/convi.c Normal file
View File

File diff suppressed because it is too large Load Diff

10
libvips/convolution/convolution.c
View File

@ -157,16 +157,24 @@ void
vips_convolution_operation_init( void )
{
extern int vips_conv_get_type( void );
extern int vips_compass_get_type( void );
extern int vips_conva_get_type( void );
extern int vips_convf_get_type( void );
extern int vips_convi_get_type( void );
extern int vips_convsep_get_type( void );
extern int vips_convasep_get_type( void );
extern int vips_compass_get_type( void );
extern int vips_fastcor_get_type( void );
extern int vips_spcor_get_type( void );
extern int vips_sharpen_get_type( void );
extern int vips_gaussblur_get_type( void );
vips_conv_get_type();
vips_conva_get_type();
vips_convf_get_type();
vips_convi_get_type();
vips_compass_get_type();
vips_convsep_get_type();
vips_convasep_get_type();
vips_fastcor_get_type();
vips_spcor_get_type();
vips_sharpen_get_type();

46
libvips/convolution/convsep.c
View File

@ -63,6 +63,8 @@ vips_convsep_build( VipsObject *object )
VipsImage **t = (VipsImage **)
vips_object_local_array( object, 3 );
VipsImage *in;
g_object_set( convsep, "out", vips_image_new(), NULL );
if( VIPS_OBJECT_CLASS( vips_convsep_parent_class )->build( object ) )
@ -71,20 +73,38 @@ vips_convsep_build( VipsObject *object )
if( vips_check_separable( class->nickname, convolution->M ) )
return( -1 );
if( vips_rot( convolution->M, &t[0], VIPS_ANGLE_D90, NULL ) ||
vips_conv( convolution->in, &t[1], convolution->M,
"precision", convsep->precision,
"layers", convsep->layers,
"cluster", convsep->cluster,
NULL ) ||
vips_conv( t[1], &t[2], t[0],
"precision", convsep->precision,
"layers", convsep->layers,
"cluster", convsep->cluster,
NULL ) )
return( -1 );
in = convolution->in;
if( vips_image_write( t[2], convolution->out ) )
if( convsep->precision == VIPS_PRECISION_APPROXIMATE ) {
if( vips_convasep( convolution->in, &t[0], convolution->M,
"layers", convsep->layers,
NULL ) )
return( -1 );
in = t[0];
}
else {
if( vips_rot( convolution->M, &t[0], VIPS_ANGLE_D90, NULL ) )
return( -1 );
/* We must only add the offset once.
*/
vips_image_set_double( t[0], "offset", 0 );
if( vips_conv( convolution->in, &t[1], convolution->M,
"precision", convsep->precision,
"layers", convsep->layers,
"cluster", convsep->cluster,
NULL ) ||
vips_conv( t[1], &t[2], t[0],
"precision", convsep->precision,
"layers", convsep->layers,
"cluster", convsep->cluster,
NULL ) )
return( -1 );
in = t[2];
}
if( vips_image_write( in, convolution->out ) )
return( -1 );
return( 0 );

1273
libvips/convolution/im_aconv.c
View File

File diff suppressed because it is too large Load Diff

879
libvips/convolution/im_aconvsep.c
View File

@ -1,879 +0,0 @@
/* im_aconvsep ... separable approximate convolution
*
* This operation does an approximate, seperable convolution.
*
* Author: John Cupitt & Nicolas Robidoux
* Written on: 31/5/11
* Modified on:
* 31/5/11
* - from im_conv()
*/
/*
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
*/
/*
See:
http://incubator.quasimondo.com/processing/stackblur.pde
This thing is a little like stackblur, but generalised to any separable
mask.
*/
/*
TODO
- are we handling mask offset correctly?
*/
/* Show sample pixels as they are transformed.
#define DEBUG_PIXELS
*/
/*
#define DEBUG
#define VIPS_DEBUG
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif /*HAVE_CONFIG_H*/
#include <vips/intl.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <math.h>
#include <vips/vips.h>
#include <vips/vector.h>
#include <vips/debug.h>
/* Maximum number of lines we can break the mask into.
*/
#define MAX_LINES (1000)
/* Euclid's algorithm. Use this to common up mults.
*/
static int
gcd( int a, int b )
{
if( b == 0 )
return( abs( a ) );
else
return( gcd( b, a % b ) );
}
/* A set of lines.
*/
typedef struct _Lines {
/* Copy of our arguments.
*/
IMAGE *in;
IMAGE *out;
DOUBLEMASK *mask;
int n_layers;
int area;
int rounding;
/* Start is the left-most pixel in the line, end is one beyond the
* right-most pixel.
*/
int n_lines;
int start[MAX_LINES];
int end[MAX_LINES];
int factor[MAX_LINES];
} Lines;
static void
lines_start( Lines *lines, int x, int factor )
{
lines->start[lines->n_lines] = x;
lines->factor[lines->n_lines] = factor;
}
static int
lines_end( Lines *lines, int x )
{
lines->end[lines->n_lines] = x;
if( lines->n_lines >= MAX_LINES - 1 ) {
vips_error( "im_aconvsep", "%s", _( "mask too complex" ) );
return( -1 );
}
lines->n_lines += 1;
return( 0 );
}
/* Break a mask into lines.
*/
static Lines *
lines_new( IMAGE *in, IMAGE *out, DOUBLEMASK *mask, int n_layers )
{
const int width = mask->xsize * mask->ysize;
Lines *lines;
double max;
double min;
double depth;
double sum;
int layers_above;
int layers_below;
int z, n, x;
/* Check parameters.
*/
if( im_piocheck( in, out ) ||
im_check_uncoded( "im_aconvsep", in ) ||
vips_check_dmask_1d( "im_aconvsep", mask ) )
return( NULL );
lines = VIPS_NEW( out, Lines );
lines->in = in;
lines->out = out;
if( !(lines->mask = (DOUBLEMASK *) im_local( out,
(im_construct_fn) im_dup_dmask,
(im_callback_fn) im_free_dmask, mask, mask->filename, NULL )) )
return( NULL );
lines->n_layers = n_layers;
lines->n_lines = 0;
VIPS_DEBUG_MSG( "lines_new: breaking into %d layers ...\n", n_layers );
/* Find mask range. We must always include the zero axis in the mask.
*/
max = 0;
min = 0;
for( x = 0; x < width; x++ ) {
if( mask->coeff[x] > max )
max = mask->coeff[x];
if( mask->coeff[x] < min )
min = mask->coeff[x];
}
/* The zero axis must fall on a layer boundary. Estimate the
* depth, find n-lines-above-zero, get exact depth, then calculate a
* fixed n-lines which includes any negative parts.
*/
depth = (max - min) / n_layers;
layers_above = ceil( max / depth );
depth = max / layers_above;
layers_below = floor( min / depth );
n_layers = layers_above - layers_below;
VIPS_DEBUG_MSG( "depth = %g, n_layers = %d\n", depth, n_layers );
/* For each layer, generate a set of lines which are inside the
* perimeter. Work down from the top.
*/
for( z = 0; z < n_layers; z++ ) {
double y = max - (1 + z) * depth;
/* y plus half depth ... ie. the layer midpoint.
*/
double y_ph = y + depth / 2;
/* Odd, but we must avoid rounding errors that make us miss 0
* in the line above.
*/
int y_positive = z < layers_above;
int inside;
/* Start outside the perimeter.
*/
inside = 0;
for( x = 0; x < width; x++ ) {
/* The vertical line from mask[z] to 0 is inside. Is
* our current square (x, y) part of that line?
*/
if( (y_positive && mask->coeff[x] >= y_ph) ||
(!y_positive && mask->coeff[x] <= y_ph) ) {
if( !inside ) {
lines_start( lines, x,
y_positive ? 1 : -1 );
inside = 1;
}
}
else {
if( inside ) {
if( lines_end( lines, x ) )
return( NULL );
inside = 0;
}
}
}
if( inside &&
lines_end( lines, width ) )
return( NULL );
}
/* Can we common up any lines? Search for lines with identical
* start/end.
*/
for( z = 0; z < lines->n_lines; z++ ) {
for( n = z + 1; n < lines->n_lines; n++ ) {
if( lines->start[z] == lines->start[n] &&
lines->end[z] == lines->end[n] ) {
lines->factor[z] += lines->factor[n];
/* n can be deleted. Do this in a separate
* pass below.
*/
lines->factor[n] = 0;
}
}
}
/* Now we can remove all factor 0 lines.
*/
for( z = 0; z < lines->n_lines; z++ ) {
if( lines->factor[z] == 0 ) {
for( x = z; x < lines->n_lines; x++ ) {
lines->start[x] = lines->start[x + 1];
lines->end[x] = lines->end[x + 1];
lines->factor[x] = lines->factor[x + 1];
}
lines->n_lines -= 1;
}
}
/* Find the area of the lines.
*/
lines->area = 0;
for( z = 0; z < lines->n_lines; z++ )
lines->area += lines->factor[z] *
(lines->end[z] - lines->start[z]);
/* Strength reduction: if all lines are divisible by n, we can move
* that n out into the ->area factor. The aim is to produce as many
* factor 1 lines as we can and to reduce the chance of overflow.
*/
x = lines->factor[0];
for( z = 1; z < lines->n_lines; z++ )
x = gcd( x, lines->factor[z] );
for( z = 0; z < lines->n_lines; z++ )
lines->factor[z] /= x;
lines->area *= x;
/* Find the area of the original mask.
*/
sum = 0;
for( z = 0; z < width; z++ )
sum += mask->coeff[z];
lines->area = rint( sum * lines->area / mask->scale );
lines->rounding = (lines->area + 1) / 2 + mask->offset * lines->area;
/* ASCII-art layer drawing.
printf( "lines:\n" );
for( z = 0; z < lines->n_lines; z++ ) {
printf( "%3d - %2d x ", z, lines->factor[z] );
for( x = 0; x < 55; x++ ) {
int rx = x * (width + 1) / 55;
if( rx >= lines->start[z] && rx < lines->end[z] )
printf( "#" );
else
printf( " " );
}
printf( " %3d .. %3d\n", lines->start[z], lines->end[z] );
}
printf( "area = %d\n", lines->area );
printf( "rounding = %d\n", lines->rounding );
*/
return( lines );
}
/* Our sequence value.
*/
typedef struct {
Lines *lines;
REGION *ir; /* Input region */
int *start; /* Offsets for start and stop */
int *end;
/* The sums for each line. int for integer types, double for floating
* point types.
*/
void *sum;
int last_stride; /* Avoid recalcing offsets, if we can */
} AConvSep;
/* Free a sequence value.
*/
static int
aconvsep_stop( void *vseq, void *a, void *b )
{
AConvSep *seq = (AConvSep *) vseq;
IM_FREEF( im_region_free, seq->ir );
return( 0 );
}
/* Convolution start function.
*/
static void *
aconvsep_start( IMAGE *out, void *a, void *b )
{
IMAGE *in = (IMAGE *) a;
Lines *lines = (Lines *) b;
AConvSep *seq;
if( !(seq = IM_NEW( out, AConvSep )) )
return( NULL );
/* Init!
*/
seq->lines = lines;
seq->ir = im_region_create( in );
seq->start = IM_ARRAY( out, lines->n_lines, int );
seq->end = IM_ARRAY( out, lines->n_lines, int );
if( vips_band_format_isint( out->BandFmt ) )
seq->sum = IM_ARRAY( out, lines->n_lines, int );
else
seq->sum = IM_ARRAY( out, lines->n_lines, double );
seq->last_stride = -1;
if( !seq->ir || !seq->start || !seq->end || !seq->sum ) {
aconvsep_stop( seq, in, lines );
return( NULL );
}
return( seq );
}
#define CLIP_UCHAR( V ) \
G_STMT_START { \
if( (V) < 0 ) \
(V) = 0; \
else if( (V) > UCHAR_MAX ) \
(V) = UCHAR_MAX; \
} G_STMT_END
#define CLIP_CHAR( V ) \
G_STMT_START { \
if( (V) < SCHAR_MIN ) \
(V) = SCHAR_MIN; \
else if( (V) > SCHAR_MAX ) \
(V) = SCHAR_MAX; \
} G_STMT_END
#define CLIP_USHORT( V ) \
G_STMT_START { \
if( (V) < 0 ) \
(V) = 0; \
else if( (V) > USHRT_MAX ) \
(V) = USHRT_MAX; \
} G_STMT_END
#define CLIP_SHORT( V ) \
G_STMT_START { \
if( (V) < SHRT_MIN ) \
(V) = SHRT_MIN; \
else if( (V) > SHRT_MAX ) \
(V) = SHRT_MAX; \
} G_STMT_END
#define CLIP_NONE( V ) {}
/* The h and v loops are very similar, but also annoyingly different. Keep
* them separate for easy debugging.
*/
#define HCONV_INT( TYPE, CLIP ) { \
for( i = 0; i < bands; i++ ) { \
int *seq_sum = (int *) seq->sum; \
\
TYPE *q; \
TYPE *p; \
int sum; \
\
p = i + (TYPE *) IM_REGION_ADDR( ir, r->left, r->top + y ); \
q = i + (TYPE *) IM_REGION_ADDR( or, r->left, r->top + y ); \
\
sum = 0; \
for( z = 0; z < n_lines; z++ ) { \
seq_sum[z] = 0; \
for( x = lines->start[z]; x < lines->end[z]; x++ ) \
seq_sum[z] += p[x * istride]; \
sum += lines->factor[z] * seq_sum[z]; \
} \
sum = (sum + lines->rounding) / lines->area; \
CLIP( sum ); \
*q = sum; \
q += ostride; \
\
for( x = 1; x < r->width; x++ ) { \
sum = 0; \
for( z = 0; z < n_lines; z++ ) { \
seq_sum[z] += p[seq->end[z]]; \
seq_sum[z] -= p[seq->start[z]]; \
sum += lines->factor[z] * seq_sum[z]; \
} \
p += istride; \
sum = (sum + lines->rounding) / lines->area; \
CLIP( sum ); \
*q = sum; \
q += ostride; \
} \
} \
}
#define HCONV_FLOAT( TYPE ) { \
for( i = 0; i < bands; i++ ) { \
double *seq_sum = (double *) seq->sum; \
\
TYPE *q; \
TYPE *p; \
double sum; \
\
p = i + (TYPE *) IM_REGION_ADDR( ir, r->left, r->top + y ); \
q = i + (TYPE *) IM_REGION_ADDR( or, r->left, r->top + y ); \
\
sum = 0; \
for( z = 0; z < lines->n_lines; z++ ) { \
seq_sum[z] = 0; \
for( x = lines->start[z]; x < lines->end[z]; x++ ) \
seq_sum[z] += p[x * istride]; \
sum += lines->factor[z] * seq_sum[z]; \
} \
sum = sum / lines->area + mask->offset; \
*q = sum; \
q += ostride; \
\
for( x = 1; x < r->width; x++ ) { \
sum = 0; \
for( z = 0; z < lines->n_lines; z++ ) { \
seq_sum[z] += p[seq->end[z]]; \
seq_sum[z] -= p[seq->start[z]]; \
sum += lines->factor[z] * seq_sum[z]; \
} \
p += istride; \
sum = sum / lines->area + mask->offset; \
*q = sum; \
q += ostride; \
} \
} \
}
/* Do horizontal masks ... we scan the mask along scanlines.
*/
static int
aconvsep_generate_horizontal( REGION *or, void *vseq, void *a, void *b )
{
AConvSep *seq = (AConvSep *) vseq;
IMAGE *in = (IMAGE *) a;
Lines *lines = (Lines *) b;
REGION *ir = seq->ir;
const int n_lines = lines->n_lines;
DOUBLEMASK *mask = lines->mask;
Rect *r = &or->valid;
/* Double the bands (notionally) for complex.
*/
int bands = vips_band_format_iscomplex( in->BandFmt ) ?
2 * in->Bands : in->Bands;
Rect s;
int x, y, z, i;
int istride;
int ostride;
/* Prepare the section of the input image we need. A little larger
* than the section of the output image we are producing.
*/
s = *r;
s.width += mask->xsize - 1;
s.height += mask->ysize - 1;
if( im_prepare( ir, &s ) )
return( -1 );
/* Stride can be different for the vertical case, keep this here for
* ease of direction change.
*/
istride = IM_IMAGE_SIZEOF_PEL( in ) /
IM_IMAGE_SIZEOF_ELEMENT( in );
ostride = IM_IMAGE_SIZEOF_PEL( lines->out ) /
IM_IMAGE_SIZEOF_ELEMENT( lines->out );
/* Init offset array.
*/
if( seq->last_stride != istride ) {
seq->last_stride = istride;
for( z = 0; z < n_lines; z++ ) {
seq->start[z] = lines->start[z] * istride;
seq->end[z] = lines->end[z] * istride;
}
}
for( y = 0; y < r->height; y++ ) {
switch( in->BandFmt ) {
case IM_BANDFMT_UCHAR:
HCONV_INT( unsigned char, CLIP_UCHAR );
break;
case IM_BANDFMT_CHAR:
HCONV_INT( signed char, CLIP_UCHAR );
break;
case IM_BANDFMT_USHORT:
HCONV_INT( unsigned short, CLIP_USHORT );
break;
case IM_BANDFMT_SHORT:
HCONV_INT( signed short, CLIP_SHORT );
break;
case IM_BANDFMT_UINT:
HCONV_INT( unsigned int, CLIP_NONE );
break;
case IM_BANDFMT_INT:
HCONV_INT( signed int, CLIP_NONE );
break;
case IM_BANDFMT_FLOAT:
HCONV_FLOAT( float );
break;
case IM_BANDFMT_DOUBLE:
HCONV_FLOAT( double );
break;
case IM_BANDFMT_COMPLEX:
HCONV_FLOAT( float );
break;
case IM_BANDFMT_DPCOMPLEX:
HCONV_FLOAT( double );
break;
default:
g_assert_not_reached();
}
}
return( 0 );
}
#define VCONV_INT( TYPE, CLIP ) { \
for( x = 0; x < sz; x++ ) { \
int *seq_sum = (int *) seq->sum; \
\
TYPE *q; \
TYPE *p; \
int sum; \
\
p = x + (TYPE *) IM_REGION_ADDR( ir, r->left, r->top ); \
q = x + (TYPE *) IM_REGION_ADDR( or, r->left, r->top ); \
\
sum = 0; \
for( z = 0; z < lines->n_lines; z++ ) { \
seq_sum[z] = 0; \
for( y = lines->start[z]; y < lines->end[z]; y++ ) \
seq_sum[z] += p[y * istride]; \
sum += lines->factor[z] * seq_sum[z]; \
} \
sum = (sum + lines->rounding) / lines->area; \
CLIP( sum ); \
*q = sum; \
q += ostride; \
\
for( y = 1; y < r->height; y++ ) { \
sum = 0; \
for( z = 0; z < lines->n_lines; z++ ) { \
seq_sum[z] += p[seq->end[z]]; \
seq_sum[z] -= p[seq->start[z]]; \
sum += lines->factor[z] * seq_sum[z]; \
} \
p += istride; \
sum = (sum + lines->rounding) / lines->area; \
CLIP( sum ); \
*q = sum; \
q += ostride; \
} \
} \
}
#define VCONV_FLOAT( TYPE ) { \
for( x = 0; x < sz; x++ ) { \
double *seq_sum = (double *) seq->sum; \
\
TYPE *q; \
TYPE *p; \
double sum; \
\
p = x + (TYPE *) IM_REGION_ADDR( ir, r->left, r->top ); \
q = x + (TYPE *) IM_REGION_ADDR( or, r->left, r->top ); \
\
sum = 0; \
for( z = 0; z < lines->n_lines; z++ ) { \
seq_sum[z] = 0; \
for( y = lines->start[z]; y < lines->end[z]; y++ ) \
seq_sum[z] += p[y * istride]; \
sum += lines->factor[z] * seq_sum[z]; \
} \
sum = sum / lines->area + mask->offset; \
*q = sum; \
q += ostride; \
\
for( y = 1; y < r->height; y++ ) { \
sum = 0; \
for( z = 0; z < lines->n_lines; z++ ) { \
seq_sum[z] += p[seq->end[z]]; \
seq_sum[z] -= p[seq->start[z]]; \
sum += lines->factor[z] * seq_sum[z]; \
} \
p += istride; \
sum = sum / lines->area + mask->offset; \
*q = sum; \
q += ostride; \
} \
} \
}
/* Do vertical masks ... we scan the mask down columns of pixels. Copy-paste
* from above with small changes.
*/
static int
aconvsep_generate_vertical( REGION *or, void *vseq, void *a, void *b )
{
AConvSep *seq = (AConvSep *) vseq;
IMAGE *in = (IMAGE *) a;
Lines *lines = (Lines *) b;
REGION *ir = seq->ir;
const int n_lines = lines->n_lines;
DOUBLEMASK *mask = lines->mask;
Rect *r = &or->valid;
/* Double the width (notionally) for complex.
*/
int sz = vips_band_format_iscomplex( in->BandFmt ) ?
2 * IM_REGION_N_ELEMENTS( or ) : IM_REGION_N_ELEMENTS( or );
Rect s;
int x, y, z;
int istride;
int ostride;
/* Prepare the section of the input image we need. A little larger
* than the section of the output image we are producing.
*/
s = *r;
s.width += mask->xsize - 1;
s.height += mask->ysize - 1;
if( im_prepare( ir, &s ) )
return( -1 );
/* Stride can be different for the vertical case, keep this here for
* ease of direction change.
*/
istride = IM_REGION_LSKIP( ir ) /
IM_IMAGE_SIZEOF_ELEMENT( lines->in );
ostride = IM_REGION_LSKIP( or ) /
IM_IMAGE_SIZEOF_ELEMENT( lines->out );
/* Init offset array.
*/
if( seq->last_stride != istride ) {
seq->last_stride = istride;
for( z = 0; z < n_lines; z++ ) {
seq->start[z] = lines->start[z] * istride;
seq->end[z] = lines->end[z] * istride;
}
}
switch( in->BandFmt ) {
case IM_BANDFMT_UCHAR:
VCONV_INT( unsigned char, CLIP_UCHAR );
break;
case IM_BANDFMT_CHAR:
VCONV_INT( signed char, CLIP_UCHAR );
break;
case IM_BANDFMT_USHORT:
VCONV_INT( unsigned short, CLIP_USHORT );
break;
case IM_BANDFMT_SHORT:
VCONV_INT( signed short, CLIP_SHORT );
break;
case IM_BANDFMT_UINT:
VCONV_INT( unsigned int, CLIP_NONE );
break;
case IM_BANDFMT_INT:
VCONV_INT( signed int, CLIP_NONE );
break;
case IM_BANDFMT_FLOAT:
VCONV_FLOAT( float );
break;
case IM_BANDFMT_DOUBLE:
VCONV_FLOAT( double );
break;
case IM_BANDFMT_COMPLEX:
VCONV_FLOAT( float );
break;
case IM_BANDFMT_DPCOMPLEX:
VCONV_FLOAT( double );
break;
default:
g_assert_not_reached();
}
return( 0 );
}
static int
aconvsep_raw( IMAGE *in, IMAGE *out, DOUBLEMASK *mask, int n_layers )
{
Lines *lines;
im_generate_fn generate;
#ifdef DEBUG
printf( "aconvsep_raw: starting with matrix:\n" );
im_print_dmask( mask );
#endif /*DEBUG*/
if( !(lines = lines_new( in, out, mask, n_layers )) )
return( -1 );
/* Prepare output. Consider a 7x7 mask and a 7x7 image --- the output
* would be 1x1.
*/
if( im_cp_desc( out, in ) )
return( -1 );
out->Xsize -= mask->xsize - 1;
out->Ysize -= mask->ysize - 1;
if( out->Xsize <= 0 || out->Ysize <= 0 ) {
im_error( "im_aconvsep", "%s", _( "image too small for mask" ) );
return( -1 );
}
if( mask->xsize == 1 )
generate = aconvsep_generate_vertical;
else
generate = aconvsep_generate_horizontal;
if( im_demand_hint( out, IM_SMALLTILE, in, NULL ) ||
im_generate( out,
aconvsep_start, generate, aconvsep_stop, in, lines ) )
return( -1 );
out->Xoffset = -mask->xsize / 2;
out->Yoffset = -mask->ysize / 2;
return( 0 );
}
/**
* im_aconvsep:
* @in: input image
* @out: output image
* @mask: convolution mask
* @n_layers: number of layers for approximation
*
* Perform an approximate separable convolution of @in with @mask.
*
* The mask must be 1xn or nx1 elements.
* The output image
* always has the same #VipsBandFormat as the input image.
*
* The image is convolved twice: once with @mask and then again with @mask
* rotated by 90 degrees.
*
* Larger values for @n_layers give more accurate
* results, but are slower. As @n_layers approaches the mask radius, the
* accuracy will become close to exact convolution and the speed will drop to
* match. For many large masks, such as Gaussian, @n_layers need be only 10% of
* this value and accuracy will still be good.
*
* See also: im_convsep_f(), im_create_dmaskv().
*
* Returns: 0 on success, -1 on error
*/
int
im_aconvsep( IMAGE *in, IMAGE *out, DOUBLEMASK *mask, int n_layers )
{
IMAGE *t[2];
const int n_mask = mask->xsize * mask->ysize;
DOUBLEMASK *rmask;
if( im_open_local_array( out, t, 2, "im_aconvsep", "p" ) ||
!(rmask = (DOUBLEMASK *) im_local( out,
(im_construct_fn) im_dup_dmask,
(im_callback_fn) im_free_dmask, mask, mask->filename, NULL )) )
return( -1 );
rmask->xsize = mask->ysize;
rmask->ysize = mask->xsize;
/*
*/
if( im_embed( in, t[0], 1, n_mask / 2, n_mask / 2,
in->Xsize + n_mask - 1, in->Ysize + n_mask - 1 ) ||
aconvsep_raw( t[0], t[1], mask, n_layers ) ||
aconvsep_raw( t[1], out, rmask, n_layers ) )
return( -1 );
/* For testing .. just try one direction.
if( aconvsep_raw( in, out, mask, n_layers ) )
return( -1 );
*/
out->Xoffset = 0;
out->Yoffset = 0;
return( 0 );
}

1106
libvips/convolution/im_conv.c
View File

File diff suppressed because it is too large Load Diff

389
libvips/convolution/im_conv_f.c
View File

@ -1,389 +0,0 @@
/* im_conv_f
*
* Copyright: 1990, N. Dessipris.
*
* Author: Nicos Dessipris & Kirk Martinez
* Written on: 29/04/1991
* Modified on: 19/05/1991
* 8/7/93 JC
* - adapted for partial v2
* - memory leaks fixed
* - ANSIfied
* 12/7/93 JC
* - adapted im_convbi() to im_convbf()
* 7/10/94 JC
* - new IM_ARRAY() macro
* - evalend callbacks
* - more typedef
* 9/3/01 JC
* - redone from im_conv()
* 27/7/01 JC
* - rejects masks with scale == 0
* 7/4/04
* - now uses im_embed() with edge stretching on the input, not
* the output
* - sets Xoffset / Yoffset
* 11/11/05
* - simpler inner loop avoids gcc4 bug
* 12/11/09
* - only rebuild the buffer offsets if bpl changes
* - tiny speedups and cleanups
* - add restrict, though it doesn't seem to help gcc
* - add mask-all-zero check
* 13/11/09
* - rename as im_conv_f() to make it easier for vips.c to make the
* overloaded version
* 3/2/10
* - gtkdoc
* - more cleanups
* 1/10/10
* - support complex (just double the bands)
* 29/10/10
* - get rid of im_convsep_f(), just call this twice, no longer worth
* keeping two versions
* 15/10/11 Nicolas
* - handle offset correctly in seperable convolutions
* 26/1/16 Lovell Fuller
* - remove Duff for a 25% speedup
*/
/*
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 <stdlib.h>
#include <limits.h>
#include <vips/vips.h>
/* Our parameters ... we take a copy of the mask argument, plus we make a
* smaller version with the zeros squeezed out.
*/
typedef struct {
IMAGE *in;
IMAGE *out;
DOUBLEMASK *mask; /* Copy of mask arg */
int nnz; /* Number of non-zero mask elements */
double *coeff; /* Array of non-zero mask coefficients */
int *coeff_pos; /* Index of each nnz element in mask->coeff */
} Conv;
static int
conv_close( Conv *conv )
{
IM_FREEF( im_free_dmask, conv->mask );
return( 0 );
}
static Conv *
conv_new( IMAGE *in, IMAGE *out, DOUBLEMASK *mask )
{
Conv *conv = IM_NEW( out, Conv );
const int ne = mask->xsize * mask->ysize;
int i;
if( !conv )
return( NULL );
conv->in = in;
conv->out = out;
conv->mask = NULL;
conv->nnz = 0;
conv->coeff = NULL;
if( im_add_close_callback( out,
(im_callback_fn) conv_close, conv, NULL ) ||
!(conv->coeff = IM_ARRAY( out, ne, double )) ||
!(conv->coeff_pos = IM_ARRAY( out, ne, int )) ||
!(conv->mask = im_dup_dmask( mask, "conv_mask" )) )
return( NULL );
/* Find non-zero mask elements.
*/
for( i = 0; i < ne; i++ )
if( mask->coeff[i] ) {
conv->coeff[conv->nnz] = mask->coeff[i];
conv->coeff_pos[conv->nnz] = i;
conv->nnz += 1;
}
/* Was the whole mask zero? We must have at least 1 element in there:
* set it to zero.
*/
if( conv->nnz == 0 ) {
conv->coeff[0] = mask->coeff[0];
conv->coeff_pos[0] = 0;
conv->nnz = 1;
}
return( conv );
}
/* Our sequence value.
*/
typedef struct {
Conv *conv;
REGION *ir; /* Input region */
int *offsets; /* Offsets for each non-zero matrix element */
VipsPel **pts; /* Per-non-zero mask element image pointers */
int last_bpl; /* Avoid recalcing offsets, if we can */
} ConvSequence;
/* Free a sequence value.
*/
static int
conv_stop( void *vseq, void *a, void *b )
{
ConvSequence *seq = (ConvSequence *) vseq;
IM_FREEF( im_region_free, seq->ir );
return( 0 );
}
/* Convolution start function.
*/
static void *
conv_start( IMAGE *out, void *a, void *b )
{
IMAGE *in = (IMAGE *) a;
Conv *conv = (Conv *) b;
ConvSequence *seq;
if( !(seq = IM_NEW( out, ConvSequence )) )
return( NULL );
/* Init!
*/
seq->conv = conv;
seq->ir = NULL;
seq->pts = NULL;
seq->last_bpl = -1;
/* Attach region and arrays.
*/
seq->ir = im_region_create( in );
seq->offsets = IM_ARRAY( out, conv->nnz, int );
seq->pts = IM_ARRAY( out, conv->nnz, VipsPel * );
if( !seq->ir || !seq->offsets || !seq->pts ) {
conv_stop( seq, in, conv );
return( NULL );
}
return( (void *) seq );
}
#define CONV_FLOAT( ITYPE, OTYPE ) { \
ITYPE ** restrict p = (ITYPE **) seq->pts; \
OTYPE * restrict q = (OTYPE *) IM_REGION_ADDR( or, le, y ); \
\
for( x = 0; x < sz; x++ ) { \
double sum; \
int i; \
\
sum = 0; \
for ( i = 0; i < nnz; i++ ) \
sum += t[i] * p[i][x]; \
\
sum = (sum / mask->scale) + mask->offset; \
\
q[x] = sum; \
} \
}
/* Convolve!
*/
static int
conv_gen( REGION *or, void *vseq, void *a, void *b )
{
ConvSequence *seq = (ConvSequence *) vseq;
IMAGE *in = (IMAGE *) a;
Conv *conv = (Conv *) b;
REGION *ir = seq->ir;
DOUBLEMASK *mask = conv->mask;
double * restrict t = conv->coeff;
const int nnz = conv->nnz;
Rect *r = &or->valid;
Rect s;
int le = r->left;
int to = r->top;
int bo = IM_RECT_BOTTOM(r);
int sz = IM_REGION_N_ELEMENTS( or ) *
(vips_band_format_iscomplex( in->BandFmt ) ? 2 : 1);
int x, y, z, i;
/* Prepare the section of the input image we need. A little larger
* than the section of the output image we are producing.
*/
s = *r;
s.width += mask->xsize - 1;
s.height += mask->ysize - 1;
if( im_prepare( ir, &s ) )
return( -1 );
/* Fill offset array. Only do this if the bpl has changed since the
* previous im_prepare().
*/
if( seq->last_bpl != IM_REGION_LSKIP( ir ) ) {
seq->last_bpl = IM_REGION_LSKIP( ir );
for( i = 0; i < nnz; i++ ) {
z = conv->coeff_pos[i];
x = z % conv->mask->xsize;
y = z / conv->mask->xsize;
seq->offsets[i] =
IM_REGION_ADDR( ir, x + le, y + to ) -
IM_REGION_ADDR( ir, le, to );
}
}
for( y = to; y < bo; y++ ) {
/* Init pts for this line of PELs.
*/
for( z = 0; z < nnz; z++ )
seq->pts[z] = seq->offsets[z] +
IM_REGION_ADDR( ir, le, y );
switch( in->BandFmt ) {
case IM_BANDFMT_UCHAR:
CONV_FLOAT( unsigned char, float ); break;
case IM_BANDFMT_CHAR:
CONV_FLOAT( signed char, float ); break;
case IM_BANDFMT_USHORT:
CONV_FLOAT( unsigned short, float ); break;
case IM_BANDFMT_SHORT:
CONV_FLOAT( signed short, float ); break;
case IM_BANDFMT_UINT:
CONV_FLOAT( unsigned int, float ); break;
case IM_BANDFMT_INT:
CONV_FLOAT( signed int, float ); break;
case IM_BANDFMT_FLOAT:
case IM_BANDFMT_COMPLEX:
CONV_FLOAT( float, float ); break;
case IM_BANDFMT_DOUBLE:
case IM_BANDFMT_DPCOMPLEX:
CONV_FLOAT( double, double ); break;
default:
g_assert_not_reached();
}
}
return( 0 );
}
int
im_conv_f_raw( IMAGE *in, IMAGE *out, DOUBLEMASK *mask )
{
Conv *conv;
/* Check parameters.
*/
if( im_piocheck( in, out ) ||
im_check_uncoded( "im_conv", in ) ||
im_check_dmask( "im_conv", mask ) )
return( -1 );
if( mask->scale == 0 ) {
im_error( "im_conv_f", "%s", "mask scale must be non-zero" );
return( -1 );
}
if( !(conv = conv_new( in, out, mask )) )
return( -1 );
/* Prepare output. Consider a 7x7 mask and a 7x7 image --- the output
* would be 1x1.
*/
if( im_cp_desc( out, in ) )
return( -1 );
if( vips_bandfmt_isint( in->BandFmt ) )
out->BandFmt = IM_BANDFMT_FLOAT;
out->Xsize -= mask->xsize - 1;
out->Ysize -= mask->ysize - 1;
if( out->Xsize <= 0 || out->Ysize <= 0 ) {
im_error( "im_conv_f", "%s", _( "image too small for mask" ) );
return( -1 );
}
if( im_demand_hint( out, IM_SMALLTILE, in, NULL ) )
return( -1 );
if( im_generate( out, conv_start, conv_gen, conv_stop, in, conv ) )
return( -1 );
out->Xoffset = -mask->xsize / 2;
out->Yoffset = -mask->ysize / 2;
return( 0 );
}
/**
* im_conv_f:
* @in: input image
* @out: output image
* @mask: convolution mask
*
* Convolve @in with @mask using floating-point arithmetic. The output image
* is always %IM_BANDFMT_FLOAT unless @in is %IM_BANDFMT_DOUBLE, in which case
* @out is also %IM_BANDFMT_DOUBLE.
*
* Each output pixel is
* calculated as sigma[i]{pixel[i] * mask[i]} / scale + offset, where scale
* and offset are part of @mask.
*
* See also: im_conv(), im_convsep_f(), im_create_dmaskv().
*
* Returns: 0 on success, -1 on error
*/
int
im_conv_f( IMAGE *in, IMAGE *out, DOUBLEMASK *mask )
{
IMAGE *t1 = im_open_local( out, "im_conv_f intermediate", "p" );
if( !t1 ||
im_embed( in, t1, 1, mask->xsize / 2, mask->ysize / 2,
in->Xsize + mask->xsize - 1,
in->Ysize + mask->ysize - 1 ) ||
im_conv_f_raw( t1, out, mask ) )
return( -1 );
out->Xoffset = 0;
out->Yoffset = 0;
return( 0 );
}

17
libvips/create/worley.c
View File

@ -1,6 +1,9 @@
/* Worley noise generator.
*
* 19/7/16
*
* 11/8/16
* - float output
*/
/*
@ -215,7 +218,7 @@ vips_worley_start( VipsImage *out, void *a, void *b )
return( seq );
}
static int
static float
vips_hypot( int x, int y )
{
/* Faster than hypot() for int args.
@ -223,10 +226,10 @@ vips_hypot( int x, int y )
return( sqrt( x * x + y * y ) );
}
static int
static float
vips_worley_distance( VipsWorley *worley, Cell cells[9], int x, int y )
{
int distance;
float distance;
int i, j;
@ -236,7 +239,7 @@ vips_worley_distance( VipsWorley *worley, Cell cells[9], int x, int y )
Cell *cell = &cells[i];
for( j = 0; j < cell->n_features; j++ ) {
int d = vips_hypot(
float d = vips_hypot(
x - cell->feature_x[j],
y - cell->feature_y[j] );
@ -258,7 +261,7 @@ vips_worley_gen( VipsRegion *or, void *vseq, void *a, void *b,
int x, y;
for( y = 0; y < r->height; y++ ) {
int *q = (int *) VIPS_REGION_ADDR( or, r->left, r->top + y );
float *q = (float *) VIPS_REGION_ADDR( or, r->left, r->top + y );
for( x = 0; x < r->width; x++ ) {
int cell_x = (r->left + x) / worley->cell_size;
@ -300,7 +303,7 @@ vips_worley_build( VipsObject *object )
vips_image_init_fields( create->out,
worley->width, worley->height, 1,
VIPS_FORMAT_INT, VIPS_CODING_NONE, VIPS_INTERPRETATION_B_W,
VIPS_FORMAT_FLOAT, VIPS_CODING_NONE, VIPS_INTERPRETATION_B_W,
1.0, 1.0 );
vips_image_pipelinev( create->out,
VIPS_DEMAND_STYLE_ANY, NULL );
@ -365,7 +368,7 @@ vips_worley_init( VipsWorley *worley )
*
* * @cell_size: %gint, size of Worley cells
*
* Create a one-band int image of Worley noise. See:
* Create a one-band float image of Worley noise. See:
*
* https://en.wikipedia.org/wiki/Worley_noise
*

110
libvips/deprecated/vips7compat.c
View File

@ -2404,6 +2404,116 @@ im_convsep_f( IMAGE *in, IMAGE *out, DOUBLEMASK *mask )
return( 0 );
}
int
im_conv( VipsImage *in, VipsImage *out, INTMASK *mask )
{
VipsImage *t1, *t2;
if( !(t1 = vips_image_new()) ||
im_imask2vips( mask, t1 ) )
return( -1 );
if( vips_convi( in, &t2, t1,
NULL ) ) {
g_object_unref( t1 );
return( -1 );
}
g_object_unref( t1 );
if( vips_image_write( t2, out ) ) {
g_object_unref( t2 );
return( -1 );
}
g_object_unref( t2 );
return( 0 );
}
int
im_conv_raw( VipsImage *in, VipsImage *out, INTMASK *mask )
{
im_error( "im_conv_raw", "no compat function" );
return( -1 );
}
int
im_conv_f( VipsImage *in, VipsImage *out, DOUBLEMASK *mask )
{
VipsImage *t1, *t2;
if( !(t1 = vips_image_new()) ||
im_mask2vips( mask, t1 ) )
return( -1 );
if( vips_convf( in, &t2, t1,
NULL ) ) {
g_object_unref( t1 );
return( -1 );
}
g_object_unref( t1 );
if( vips_image_write( t2, out ) ) {
g_object_unref( t2 );
return( -1 );
}
g_object_unref( t2 );
return( 0 );
}
int
im_aconvsep( VipsImage *in, VipsImage *out, DOUBLEMASK *mask, int n_layers )
{
VipsImage *t1, *t2;
if( !(t1 = vips_image_new()) ||
im_mask2vips( mask, t1 ) )
return( -1 );
if( vips_convasep( in, &t2, t1,
"layers", n_layers,
NULL ) ) {
g_object_unref( t1 );
return( -1 );
}
g_object_unref( t1 );
if( vips_image_write( t2, out ) ) {
g_object_unref( t2 );
return( -1 );
}
g_object_unref( t2 );
return( 0 );
}
int
im_aconv( VipsImage *in, VipsImage *out,
DOUBLEMASK *mask, int n_layers, int cluster )
{
VipsImage *t1, *t2;
if( !(t1 = vips_image_new()) ||
im_mask2vips( mask, t1 ) )
return( -1 );
if( vips_conva( in, &t2, t1,
"layers", n_layers,
"cluster", cluster,
NULL ) ) {
g_object_unref( t1 );
return( -1 );
}
g_object_unref( t1 );
if( vips_image_write( t2, out ) ) {
g_object_unref( t2 );
return( -1 );
}
g_object_unref( t2 );
return( 0 );
}
int
im_conv_f_raw( VipsImage *in, VipsImage *out, DOUBLEMASK *mask )
{
im_error( "im_conv_f_raw", "no compat function" );
return( -1 );
}
int
im_addgnoise( IMAGE *in, IMAGE *out, double sigma )
{

20
libvips/foreign/csv.c
View File

@ -30,6 +30,8 @@
* 4/6/15
* - try to support DOS files under linux ... we have to look for \r\n
* linebreaks
* 12/8/16
* - allow missing offset in matrix header
*/
/*
@ -522,6 +524,8 @@ read_ascii_double( FILE *fp, const char whitemap[256], double *out )
/* Read the header. Two numbers for width and height, and two optional
* numbers for scale and offset.
*
* We can have scale and no offset, in which case we assume offset = 0.
*/
static int
vips__matrix_header( char *whitemap, FILE *fp,
@ -532,6 +536,11 @@ vips__matrix_header( char *whitemap, FILE *fp,
int i;
int ch;
/* Offset defaults to zero.
*/
header[2] = 1.0;
header[3] = 0.0;
for( i = 0; i < 4 &&
(ch = read_ascii_double( fp, whitemap, &header[i] )) == 0;
i++ )
@ -556,22 +565,17 @@ vips__matrix_header( char *whitemap, FILE *fp,
"%s", _( "width / height out of range" ) );
return( -1 );
}
if( i == 3 ) {
vips_error( "mask2vips", "%s", _( "bad scale / offset" ) );
return( -1 );
}
if( (ch = read_ascii_double( fp, whitemap, &d )) != '\n' ) {
vips_error( "mask2vips", "%s", _( "extra chars in header" ) );
return( -1 );
}
if( i > 2 &&
header[2] == 0.0 ) {
if( header[2] == 0.0 ) {
vips_error( "mask2vips", "%s", _( "zero scale" ) );
return( -1 );
}
*scale = i > 2 ? header[2] : 1.0;
*offset = i > 2 ? header[3] : 0.0;
*scale = header[2];
*offset = header[3];
skip_line( fp );

18
libvips/include/vips/convolution.h
View File

@ -31,8 +31,8 @@
*/
#ifndef IM_CONVOLUTION_H
#define IM_CONVOLUTION_H
#ifndef VIPS_CONVOLUTION_H
#define VIPS_CONVOLUTION_H
#ifdef __cplusplus
extern "C" {
@ -46,15 +46,23 @@ typedef enum {
int vips_conv( VipsImage *in, VipsImage **out, VipsImage *mask, ... )
__attribute__((sentinel));
int vips_compass( VipsImage *in, VipsImage **out, VipsImage *mask, ... )
int vips_convf( VipsImage *in, VipsImage **out, VipsImage *mask, ... )
__attribute__((sentinel));
int vips_convi( VipsImage *in, VipsImage **out, VipsImage *mask, ... )
__attribute__((sentinel));
int vips_conva( VipsImage *in, VipsImage **out, VipsImage *mask, ... )
__attribute__((sentinel));
int vips_convsep( VipsImage *in, VipsImage **out, VipsImage *mask, ... )
__attribute__((sentinel));
int vips_convasep( VipsImage *in, VipsImage **out, VipsImage *mask, ... )
__attribute__((sentinel));
int vips_sharpen( VipsImage *in, VipsImage **out, ... )
int vips_compass( VipsImage *in, VipsImage **out, VipsImage *mask, ... )
__attribute__((sentinel));
int vips_gaussblur( VipsImage *in, VipsImage **out, double sigma, ... )
__attribute__((sentinel));
int vips_sharpen( VipsImage *in, VipsImage **out, ... )
__attribute__((sentinel));
int vips_spcor( VipsImage *in, VipsImage *ref, VipsImage **out, ... )
__attribute__((sentinel));
@ -65,4 +73,4 @@ int vips_fastcor( VipsImage *in, VipsImage *ref, VipsImage **out, ... )
}
#endif /*__cplusplus*/
#endif /*IM_CONVOLUTION_H*/
#endif /*VIPS_CONVOLUTION_H*/

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

@ -246,11 +246,12 @@ typedef struct _VipsImagePixels {
gint64 npels; /* Number of pels calculated so far */
} VipsImagePixels;
int
vips__foreign_convert_saveable( VipsImage *in, VipsImage **ready,
int vips__foreign_convert_saveable( VipsImage *in, VipsImage **ready,
VipsSaveable saveable, VipsBandFormat *format, VipsCoding *coding,
VipsArrayDouble *background );
int vips__image_intize( VipsImage *in, VipsImage **out );
#ifdef __cplusplus
}
#endif /*__cplusplus*/

2
libvips/include/vips/vector.h
View File

@ -136,6 +136,8 @@ void vips_executor_set_array( VipsExecutor *executor, int var, void *value );
void vips_executor_run( VipsExecutor *executor );
void vips_vector_to_fixed_point( double *in, int *out, int n, int scale );
#ifdef __cplusplus
}
#endif /*__cplusplus*/

24
libvips/iofuncs/memory.c
View File

@ -109,7 +109,15 @@ static GMutex *vips_tracked_mutex = NULL;
* @OBJ: allocate memory local to @OBJ, or %NULL for no auto-free
* @T: type of thing to allocate
*
* Returns: A pointer of type @T *, or %NULL on error.
* Allocate memory for a thing of type @T. The memory is not
* cleared.
*
* This macro cannot fail. See vips_tracked_malloc() if you are
* allocating large amounts of memory.
*
* See also: vips_malloc().
*
* Returns: A pointer of type @T *.
*/
/**
@ -118,7 +126,15 @@ static GMutex *vips_tracked_mutex = NULL;
* @N: number of @T 's to allocate
* @T: type of thing to allocate
*
* Returns: A pointer of type @T *, or %NULL on error.
* Allocate memory for an array of objects of type @T. The memory is not
* cleared.
*
* This macro cannot fail. See vips_tracked_malloc() if you are
* allocating large amounts of memory.
*
* See also: vips_malloc().
*
* Returns: A pointer of type @T *.
*/
static void
@ -141,7 +157,7 @@ vips_malloc_cb( VipsObject *object, char *buf )
*
* See also: vips_tracked_malloc().
*
* Returns: (transfer full): a pointer to the allocated memory
* Returns: (transfer full): a pointer to the allocated memory.
*/
void *
vips_malloc( VipsObject *object, size_t size )
@ -166,7 +182,7 @@ vips_malloc( VipsObject *object, size_t size )
*
* g_strdup() a string. When @object is freed, the string will be freed for
* you. If @object is %NULL, you need to
* free the memory explicitly with g_free().
* free the memory yourself with g_free().
*
* This function cannot fail.
*

81
libvips/iofuncs/vector.c
View File

@ -522,3 +522,84 @@ vips_executor_run( VipsExecutor *executor )
orc_executor_run( &executor->executor );
#endif /*HAVE_ORC*/
}
/* Make a fixed-point version of a matrix. Each
* out[i] = rint(in[i] * adj_scale), where adj_scale is selected so that
* sum(out) = sum(in) * scale.
*
* Because of the vagaries of rint(), we can't just calc this, we have to
* iterate and converge on the best value for adj_scale.
*/
void
vips_vector_to_fixed_point( double *in, int *out, int n, int scale )
{
double fsum;
int i;
int target;
int sum;
double high;
double low;
double guess;
fsum = 0.0;
for( i = 0; i < n; i++ )
fsum += in[i];
target = VIPS_RINT( fsum * scale );
/* As we rint() each scale element, we can get up to 0.5 error.
* Therefore, by the end of the mask, we can be off by up to n/2. Our
* high and low guesses are therefore n/2 either side of the obvious
* answer.
*/
high = scale + (n + 1) / 2;
low = scale - (n + 1) / 2;
do {
guess = (high + low) / 2.0;
for( i = 0; i < n; i++ )
out[i] = VIPS_RINT( in[i] * guess );
sum = 0;
for( i = 0; i < n; i++ )
sum += out[i];
if( sum == target )
break;
if( sum < target )
low = guess;
if( sum > target )
high = guess;
/* This will typically produce about 5 iterations.
*/
} while( high - low > 0.01 );
if( sum != target ) {
/* Spread the error out thinly over the whole array. For
* example, consider the matrix:
*
* 3 3 9 0
* 1 1 1
* 1 1 1
* 1 1 1
*
* being converted with scale = 64 (convi does this). We want
* to generate a mix of 7s and 8s.
*/
int each_error = (target - sum) / n;
int extra_error = (target - sum) % n;
/* To share the residual error, we add or subtract 1 from the
* first abs(extra_error) elements.
*/
int direction = extra_error > 0 ? 1 : -1;
int n_elements = VIPS_ABS( extra_error );
for( i = 0; i < n; i++ )
out[i] += each_error;
for( i = 0; i < n_elements; i++ )
out[i] += direction;
}
}

63
libvips/resample/reducev.cpp
View File

@ -702,65 +702,6 @@ vips_reducev_vector_gen( VipsRegion *out_region, void *vseq,
return( 0 );
}
/* Make a fixed-point version of a mask. Each out[i] = rint(in[i] * adj_scale),
* where adj_scale is selected so that sum(out) = sum(in) * scale.
*
* Because of the vagaries of rint(), we can't just calc this, we have to
* iterate and converge on the best value for adj_scale.
*/
static void
vips_reducev_intize( double *in, int *out, int n, int scale )
{
double fsum;
int target;
int sum;
double high;
double low;
double guess;
fsum = 0.0;
for( int i = 0; i < n; i++ )
fsum += in[i];
target = VIPS_RINT( fsum * scale );
/* As we rint() each scale element, we can get up to 0.5 error.
* Therefore, by the end of the mask, we can be off by up to n/2. Our
* high and low guesses are therefore n/2 either side of the obvious
* answer.
*/
high = scale + (n + 1) / 2;
low = scale - (n + 1) / 2;
do {
guess = (high + low) / 2.0;
for( int i = 0; i < n; i++ )
out[i] = VIPS_RINT( in[i] * guess );
sum = 0;
for( int i = 0; i < n; i++ )
sum += out[i];
if( sum == target )
break;
if( sum < target )
low = guess;
if( sum > target )
high = guess;
/* This will typically produce about 5 iterations.
*/
} while( high - low > 0.01 );
if( sum != target )
/* We're as close as we can get ... add the remaining error to
* the centre element. Hopefully we'll get slight sharpness
* changes rather than slight brightness changes and it'll
* be less visible.
*/
out[n / 2] += target - sum;
}
static int
vips_reducev_raw( VipsReducev *reducev, VipsImage *in )
{
@ -798,7 +739,7 @@ vips_reducev_raw( VipsReducev *reducev, VipsImage *in )
if( !reducev->matrixi[y] )
return( -1 );
vips_reducev_intize(
vips_vector_to_fixed_point(
reducev->matrixf[y], reducev->matrixi[y],
reducev->n_point, VIPS_INTERPOLATE_SCALE );
}
@ -813,7 +754,7 @@ vips_reducev_raw( VipsReducev *reducev, VipsImage *in )
if( !reducev->matrixo[y] )
return( -1 );
vips_reducev_intize(
vips_vector_to_fixed_point(
reducev->matrixf[y], reducev->matrixo[y],
reducev->n_point, 64 );
}

26
test/test_convolution.py
View File

@ -96,6 +96,12 @@ class TestConvolution(unittest.TestCase):
for x, y in zip_expand(a, b):
self.assertAlmostEqual(x, y, places = places, msg = msg)
# test a pair of things which can be lists for difference less than a
# threshold
def assertLessThreshold(self, a, b, diff):
for x, y in zip_expand(a, b):
self.assertLess(abs(x - y), diff)
def setUp(self):
im = Vips.Image.mask_ideal(100, 100, 0.5, reject = True, optical = True)
self.colour = im * [1, 2, 3] + [2, 3, 4]
@ -131,6 +137,26 @@ class TestConvolution(unittest.TestCase):
true = conv(im, msk, 49, 49)
self.assertAlmostEqualObjects(result, true)
# don't test conva, it's still not done
def dont_test_conva(self):
for im in self.all_images:
for msk in self.all_masks:
print("msk:")
msk.matrixprint()
print("im.bands = %s" % im.bands)
convolved = im.conv(msk, precision = Vips.Precision.APPROXIMATE)
result = convolved(25, 50)
true = conv(im, msk, 24, 49)
print("result = %s, true = %s" % (result, true))
self.assertLessThreshold(result, true, 5)
result = convolved(50, 50)
true = conv(im, msk, 49, 49)
print("result = %s, true = %s" % (result, true))
self.assertLessThreshold(result, true, 5)
def test_compass(self):
for im in self.all_images:
for msk in self.all_masks:

2
test/test_create.py
View File

@ -446,7 +446,7 @@ class TestCreate(unittest.TestCase):
self.assertEqual(im.width, 512)
self.assertEqual(im.height, 512)
self.assertEqual(im.bands, 1)
self.assertEqual(im.format, Vips.BandFormat.INT)
self.assertEqual(im.format, Vips.BandFormat.FLOAT)
def test_perlin(self):
im = Vips.Image.perlin(512, 512)