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

VipsLinear allows complex constants

Browse Source
This commit is contained in:
John Cupitt 2014-02-22 16:08:46 +00:00
parent 3bceb5286b
commit a562f46e69
5 changed files with 363 additions and 125 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
ChangeLog
View File

@ -8,6 +8,7 @@
im_draw_smudge(), im_label_regions() as classes
- better rounding in vips_flatten()
- VipsStatistic operations are sequential
- VipsLinear has optional complex constants, added vips_linear_complex() etc.
13/2/14 started 7.38.4
- --sharpen=none option to vipsthumbnail was broken, thanks ferryfax

38
TODO
View File

@ -1,26 +1,11 @@
- ink to vec etc. should work for complex .. output or accept a double-length
vector
- ink to vec etc must have a way to give a complex constant
needs to work for eg. vips_linear() as well
eg. drawink needs a --ink_imag option with the imaginary components of the
ink
if vector.len is == bands, expand to complex by setting imaginary == 0
if vector.len == bands * 2, set imaginary from constant
if vector.len == 1, bandup with imaginary == 0
if vector.len == 2, bandup setting imaginary from vector
what about bands == 2? ambigious
we must insist vector.length == 2 or bands * 2 for complex, or have the
current behaviour of imaginary == 0 always
no, vips_linear() is better the way it is now, if we had to have different
constants for complex images, we'd need to add an extra path to all uses of
vips_linear() ... better to make a vips_linearc() (or whatever) that takes
an explicitly complex constant
maybe an optional arg to linear with the imaginary components of the
constant? that would work for draw as well
look for uses of vips__vector_to_ink() and add extra params to other places
too, eg. vips_embed(), vips_insert() etc.
@ -64,19 +49,6 @@
- use g_log() instead of vips_info()?
- object construction is threadsafe, but class construction is not
https://github.com/jcupitt/libvips/issues/64
we worked around this by adding vips_class_ping_all() to dsave build, but
this is not a good fix
find out why vips class construct fails, test on seurat, kirk's 24-core
monster
- quadratic doesn't work for order 3

439
libvips/arithmetic/linear.c
View File

@ -44,6 +44,8 @@
* - try an ORC path with the band loop unrolled
* 14/1/14
* - add uchar output option
* 21/2/14
* - add imaginary components
*/
/*
@ -84,6 +86,7 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <vips/vips.h>
@ -98,6 +101,11 @@ typedef struct _VipsLinear {
VipsArea *a;
VipsArea *b;
/* Optional imaginary part. Zero if not set.
*/
VipsArea *a_imag;
VipsArea *b_imag;
/* uchar output.
*/
gboolean uchar;
@ -107,6 +115,8 @@ typedef struct _VipsLinear {
int n;
double *a_ready;
double *b_ready;
double *a_imag_ready;
double *b_imag_ready;
} VipsLinear;
@ -123,46 +133,81 @@ vips_linear_build( VipsObject *object )
VipsLinear *linear = (VipsLinear *) object;
int bands;
VipsBandFormat format;
int i;
/* How many bands will our input image have after decoding?
/* How many bands will our input image have after decoding? Need
* format too.
*/
switch( unary->in->Coding ) {
case VIPS_CODING_RAD:
bands = 3;
format = VIPS_FORMAT_FLOAT;
break;
case VIPS_CODING_LABQ:
bands = 3;
format = VIPS_FORMAT_SHORT;
break;
default:
bands = unary->in->Bands;
format = unary->in->BandFmt;
break;
}
/* If we have a three-element vector we need to bandup the image to
/* If we have a many-element vector, we need to bandup the image to
* match.
*/
linear->n = 1;
if( linear->a )
linear->n = VIPS_MAX( linear->n, linear->b->n );
linear->n = VIPS_MAX( linear->n, linear->a->n );
if( linear->b )
linear->n = VIPS_MAX( linear->n, linear->b->n );
if( linear->a_imag )
linear->n = VIPS_MAX( linear->n, linear->a_imag->n );
if( linear->b_imag )
linear->n = VIPS_MAX( linear->n, linear->b_imag->n );
if( unary->in )
linear->n = VIPS_MAX( linear->n, bands );
arithmetic->base_bands = linear->n;
if( unary->in && linear->a && linear->b ) {
if( vips_check_vector( class->nickname,
linear->a->n, unary->in ) ||
vips_check_vector( class->nickname,
linear->b->n, unary->in ) )
if( unary->in &&
linear->a &&
vips_check_vector( class->nickname, linear->a->n, unary->in ) )
return( -1 );
if( linear->b &&
linear->a &&
vips_check_vector_length( class->nickname,
linear->b->n, linear->a->n ) )
return( -1 );
if( linear->a_imag &&
linear->a &&
vips_check_vector_length( class->nickname,
linear->a_imag->n, linear->a->n ) )
return( -1 );
if( linear->b_imag &&
linear->a &&
vips_check_vector_length( class->nickname,
linear->b_imag->n, linear->a->n ) )
return( -1 );
}
/* Make up-banded versions of our constants.
*/
linear->a_ready = VIPS_ARRAY( linear, linear->n, double );
linear->b_ready = VIPS_ARRAY( linear, linear->n, double );
/* Either complex constant can be missing, we need to default to zero.
*/
if( linear->a_imag ||
linear->b_imag ) {
linear->a_imag_ready = VIPS_ARRAY( linear, linear->n, double );
linear->b_imag_ready = VIPS_ARRAY( linear, linear->n, double );
memset( linear->a_imag_ready, 0, linear->n * sizeof( double ) );
memset( linear->b_imag_ready, 0, linear->n * sizeof( double ) );
}
for( i = 0; i < linear->n; i++ ) {
if( linear->a ) {
double *ary = (double *) linear->a->data;
@ -177,10 +222,31 @@ vips_linear_build( VipsObject *object )
linear->b_ready[i] = ary[j];
}
if( linear->a_imag ) {
double *ary = (double *) linear->a_imag->data;
int j = VIPS_MIN( i, linear->a_imag->n - 1 );
linear->a_imag_ready[i] = ary[j];
}
if( linear->b_imag ) {
double *ary = (double *) linear->b_imag->data;
int j = VIPS_MIN( i, linear->b_imag->n - 1 );
linear->b_imag_ready[i] = ary[j];
}
}
if( linear->uchar )
arithmetic->format = VIPS_FORMAT_UCHAR;
else if( linear->a_imag ||
linear->b_imag ) {
if( format == VIPS_FORMAT_DOUBLE )
arithmetic->format = VIPS_FORMAT_DPCOMPLEX;
else
arithmetic->format = VIPS_FORMAT_COMPLEX;
}
if( VIPS_OBJECT_CLASS( vips_linear_parent_class )->build( object ) )
return( -1 );
@ -188,7 +254,7 @@ vips_linear_build( VipsObject *object )
return( 0 );
}
/* Non-complex input, any output, all bands of the constant equal.
/* Non-complex input, non-complex constant, all bands of the constant equal.
*/
#define LOOP1( IN, OUT ) { \
IN * restrict p = (IN *) in[0]; \
@ -201,7 +267,7 @@ vips_linear_build( VipsObject *object )
q[x] = a1 * (OUT) p[x] + b1; \
}
/* Non-complex input, any output.
/* Non-complex input, non-complex constant, many-band constant.
*/
#define LOOPN( IN, OUT ) { \
IN * restrict p = (IN *) in[0]; \
@ -212,8 +278,25 @@ vips_linear_build( VipsObject *object )
q[i] = a[k] * (OUT) p[i] + b[k]; \
}
/* Non-complex input, complex constant, many-band constant.
*/
#define LOOPNC( IN, OUT ) { \
IN * restrict p = (IN *) in[0]; \
OUT * restrict q = (OUT *) out; \
\
for( i = 0, x = 0; x < width; x++ ) \
for( k = 0; k < nb; k++, i++ ) { \
q[0] = p[i] * a[k] + b[k]; \
q[1] = p[i] * a_imag[k] + b_imag[k]; \
q += 2; \
} \
}
#define LOOP( IN, OUT ) { \
if( linear->a->n == 1 && linear->b->n == 1 ) { \
if( linear->a_imag_ready ) { \
LOOPNC( IN, OUT ); \
} \
else if( linear->a->n == 1 && linear->b->n == 1 ) { \
LOOP1( IN, OUT ); \
} \
else { \
@ -221,7 +304,7 @@ vips_linear_build( VipsObject *object )
} \
}
/* Complex input, complex output.
/* Complex input, non-complex constant.
*/
#define LOOPCMPLXN( IN, OUT ) { \
IN * restrict p = (IN *) in[0]; \
@ -236,10 +319,42 @@ vips_linear_build( VipsObject *object )
} \
}
/* Non-complex input, any output, all bands of the constant equal, uchar
* output.
/* Complex input, complex constant.
*/
#define LOOP1uc( IN ) { \
#define LOOPCMPLXNC( IN, OUT ) { \
IN * restrict p = (IN *) in[0]; \
OUT * restrict q = (OUT *) out; \
\
for( x = 0; x < width; x++ ) \
for( k = 0; k < nb; k++ ) { \
double x1 = p[0]; \
double y1 = p[1]; \
double x2 = a[k]; \
double y2 = a_imag[k]; \
\
q[0] = x1 * x2 - y1 * y2 + b[k]; \
q[1] = x1 * y2 + x2 * y1 + b_imag[k]; \
\
q += 2; \
p += 2; \
} \
}
#define LOOPCMPLX( IN, OUT ) { \
if( linear->a_imag_ready ) { \
LOOPCMPLXNC( IN, OUT ); \
} \
else { \
LOOPCMPLXN( IN, OUT ); \
} \
}
/* Non-complex input, all bands of the constant equal, uchar output. Since we
* don't look at the imaginary component of the constant since we don't
* generate the imaginary component of the output, we work for a complex
* constant too.
*/
#define LOOP1uc( IN, DUMMY ) { \
IN * restrict p = (IN *) in[0]; \
VipsPel * restrict q = (VipsPel *) out; \
float a1 = a[0]; \
@ -253,9 +368,10 @@ vips_linear_build( VipsObject *object )
} \
}
/* Non-complex input, uchar output.
/* Non-complex input, non-complex constant, uchar output. Since we are
* outputting non-complex, we will work for a complex constant too.
*/
#define LOOPNuc( IN ) { \
#define LOOPNuc( IN, DUMMY ) { \
IN * restrict p = (IN *) in[0]; \
VipsPel * restrict q = (VipsPel *) out; \
\
@ -267,18 +383,18 @@ vips_linear_build( VipsObject *object )
} \
}
#define LOOPuc( IN ) { \
#define LOOPuc( IN, DUMMY ) { \
if( linear->a->n == 1 && linear->b->n == 1 ) { \
LOOP1uc( IN ); \
LOOP1uc( IN, DUMMY ); \
} \
else { \
LOOPNuc( IN ); \
LOOPNuc( IN, DUMMY ); \
} \
}
/* Complex input, uchar output.
/* Complex input, non-complex constant, uchar output.
*/
#define LOOPCMPLXNuc( IN ) { \
#define LOOPCMPLXNuc( IN, DUMMY ) { \
IN * restrict p = (IN *) in[0]; \
VipsPel * restrict q = (VipsPel *) out; \
\
@ -291,8 +407,62 @@ vips_linear_build( VipsObject *object )
} \
}
/* Lintra a buffer, n set of scale/offset.
/* Complex input, complex constant, uchar output.
*/
#define LOOPCMPLXNCuc( IN, DUMMY ) { \
IN * restrict p = (IN *) in[0]; \
VipsPel * restrict q = (VipsPel *) out; \
\
for( i = 0, x = 0; x < width; x++ ) \
for( k = 0; k < nb; k++, i++ ) { \
double x1 = p[0]; \
double y1 = p[1]; \
double x2 = a[k]; \
double y2 = a_imag[k]; \
double t = x1 * x2 - y1 * y2 + b[k]; \
\
q[i] = VIPS_CLIP( 0, t, 255 ); \
p += 2; \
} \
}
#define LOOPCMPLXuc( IN, OUT ) { \
if( linear->a_imag_ready ) { \
LOOPCMPLXNCuc( IN, OUT ); \
} \
else { \
LOOPCMPLXNuc( IN, OUT ); \
} \
}
#define SWITCH( REAL, CMPLX ) { \
switch( vips_image_get_format( im ) ) { \
case VIPS_FORMAT_UCHAR: \
REAL( unsigned char, float ); break; \
case VIPS_FORMAT_CHAR: \
REAL( signed char, float ); break; \
case VIPS_FORMAT_USHORT: \
REAL( unsigned short, float ); break; \
case VIPS_FORMAT_SHORT: \
REAL( signed short, float ); break; \
case VIPS_FORMAT_UINT: \
REAL( unsigned int, float ); break; \
case VIPS_FORMAT_INT: \
REAL( signed int, float ); break; \
case VIPS_FORMAT_FLOAT: \
REAL( float, float ); break; \
case VIPS_FORMAT_DOUBLE: \
REAL( double, double ); break; \
case VIPS_FORMAT_COMPLEX: \
CMPLX( float, float ); break; \
case VIPS_FORMAT_DPCOMPLEX: \
CMPLX( double, double ); break; \
\
default: \
g_assert( 0 ); \
} \
}
static void
vips_linear_buffer( VipsArithmetic *arithmetic,
VipsPel *out, VipsPel **in, int width )
@ -301,67 +471,20 @@ vips_linear_buffer( VipsArithmetic *arithmetic,
VipsLinear *linear = (VipsLinear *) arithmetic;
double * restrict a = linear->a_ready;
double * restrict b = linear->b_ready;
double * restrict a_imag = linear->a_imag_ready;
double * restrict b_imag = linear->b_imag_ready;
int nb = im->Bands;
int i, x, k;
if( linear->uchar )
switch( vips_image_get_format( im ) ) {
case VIPS_FORMAT_UCHAR:
LOOPuc( unsigned char ); break;
case VIPS_FORMAT_CHAR:
LOOPuc( signed char ); break;
case VIPS_FORMAT_USHORT:
LOOPuc( unsigned short ); break;
case VIPS_FORMAT_SHORT:
LOOPuc( signed short ); break;
case VIPS_FORMAT_UINT:
LOOPuc( unsigned int ); break;
case VIPS_FORMAT_INT:
LOOPuc( signed int ); break;
case VIPS_FORMAT_FLOAT:
LOOPuc( float ); break;
case VIPS_FORMAT_DOUBLE:
LOOPuc( double ); break;
case VIPS_FORMAT_COMPLEX:
LOOPCMPLXNuc( float ); break;
case VIPS_FORMAT_DPCOMPLEX:
LOOPCMPLXNuc( double ); break;
default:
g_assert( 0 );
if( linear->uchar ) {
SWITCH( LOOPuc, LOOPCMPLXuc );
}
else {
SWITCH( LOOP, LOOPCMPLX );
}
else
switch( vips_image_get_format( im ) ) {
case VIPS_FORMAT_UCHAR:
LOOP( unsigned char, float ); break;
case VIPS_FORMAT_CHAR:
LOOP( signed char, float ); break;
case VIPS_FORMAT_USHORT:
LOOP( unsigned short, float ); break;
case VIPS_FORMAT_SHORT:
LOOP( signed short, float ); break;
case VIPS_FORMAT_UINT:
LOOP( unsigned int, float ); break;
case VIPS_FORMAT_INT:
LOOP( signed int, float ); break;
case VIPS_FORMAT_FLOAT:
LOOP( float, float ); break;
case VIPS_FORMAT_DOUBLE:
LOOP( double, double ); break;
case VIPS_FORMAT_COMPLEX:
LOOPCMPLXN( float, float ); break;
case VIPS_FORMAT_DPCOMPLEX:
LOOPCMPLXN( double, 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
@ -373,8 +496,6 @@ vips_linear_buffer( VipsArithmetic *arithmetic,
#define D VIPS_FORMAT_DOUBLE
#define DX VIPS_FORMAT_DPCOMPLEX
/* Format doesn't change with linear.
*/
static const VipsBandFormat vips_linear_format_table[10] = {
/* UC C US S UI I F X D DX */
F, F, F, F, F, F, F, X, D, DX
@ -412,7 +533,21 @@ vips_linear_class_init( VipsLinearClass *class )
G_STRUCT_OFFSET( VipsLinear, b ),
VIPS_TYPE_ARRAY_DOUBLE );
VIPS_ARG_BOOL( class, "uchar", 112,
VIPS_ARG_BOXED( class, "a_imag", 112,
_( "a_imag" ),
_( "Multiply by this (imaginary component)" ),
VIPS_ARGUMENT_OPTIONAL_INPUT,
G_STRUCT_OFFSET( VipsLinear, a_imag ),
VIPS_TYPE_ARRAY_DOUBLE );
VIPS_ARG_BOXED( class, "b_imag", 113,
_( "b_imag" ),
_( "Add this (imaginary component)" ),
VIPS_ARGUMENT_OPTIONAL_INPUT,
G_STRUCT_OFFSET( VipsLinear, b_imag ),
VIPS_TYPE_ARRAY_DOUBLE );
VIPS_ARG_BOOL( class, "uchar", 114,
_( "uchar" ),
_( "Output should be uchar" ),
VIPS_ARGUMENT_OPTIONAL_INPUT,
@ -428,21 +563,64 @@ vips_linear_init( VipsLinear *linear )
static int
vips_linearv( VipsImage *in, VipsImage **out,
double *a, double *b, int n, va_list ap )
double *a, double *a_imag, double *b, double *b_imag, int n,
va_list ap )
{
VipsOperation *operation;
VipsArea *area_a;
VipsArea *area_b;
int result;
if( !(operation = vips_operation_new( "linear" )) )
return( -1 );
area_a = (VipsArea *) vips_array_double_new( a, n );
area_b = (VipsArea *) vips_array_double_new( b, n );
result = vips_call_split( "linear", ap, in, out, area_a, area_b );
g_object_set( operation,
"in", in,
"a", area_a,
"b", area_b,
NULL );
vips_area_unref( area_a );
vips_area_unref( area_b );
return( result );
if( a_imag ) {
VipsArea *area_a_imag;
area_a_imag = (VipsArea *) vips_array_double_new( a_imag, n );
g_object_set( operation,
"a_imag", area_a_imag,
NULL );
vips_area_unref( area_a_imag );
}
if( b_imag ) {
VipsArea *area_b_imag;
area_b_imag = (VipsArea *) vips_array_double_new( b_imag, n );
g_object_set( operation,
"b_imag", area_b_imag,
NULL );
vips_area_unref( area_b_imag );
}
(void) vips_object_set_valist( VIPS_OBJECT( operation ), ap );
if( vips_cache_operation_buildp( &operation ) ) {
vips_object_unref_outputs( VIPS_OBJECT( operation ) );
g_object_unref( operation );
return( -1 );
}
g_object_get( operation,
"out", out,
NULL );
g_object_unref( operation );
return( 0 );
}
/**
@ -457,11 +635,16 @@ vips_linearv( VipsImage *in, VipsImage **out,
* Optional arguments:
*
* @uchar: output uchar pixels
* @a_imag: #VipsArrayDouble of imaginary constants for multiplication
* @b_imag: #VipsArrayDouble of imaginary constants for addition
*
* Pass an image through a linear transform, ie. (@out = @in * @a + @b). Output
* is float for integer input, double for double input, complex for
* complex input and double complex for double complex input. Set @uchar to
* output uchar pixels.
* complex input and double complex for double complex input. If complex
* constants are specified, the output is complex, see below.
*
* Set @uchar to output uchar pixels. This is much faster than vips_linear()
* followed by vips_cast().
*
* If the arrays of constants have just one element, that constant is used for
* all image bands. If the arrays have more than one element and they have
@ -470,7 +653,12 @@ vips_linearv( VipsImage *in, VipsImage **out,
* element and the image only has a single band, the result is a many-band
* image where each band corresponds to one array element.
*
* See also: vips_linear1(), vips_add().
* Set @a_imag and @b_imag to set imagiary constants for multiplication and
* addition. If imaginary components are specified, the output is complex for
* non-double-complex inputs and double-complex for double-complex inputs.
*
* See also: vips_linear1(), vips_linear_complex(), vips_linear_complex1(),
* vips_add().
*
* Returns: 0 on success, -1 on error
*/
@ -481,7 +669,7 @@ vips_linear( VipsImage *in, VipsImage **out, double *a, double *b, int n, ... )
int result;
va_start( ap, n );
result = vips_linearv( in, out, a, b, n, ap );
result = vips_linearv( in, out, a, NULL, b, NULL, n, ap );
va_end( ap );
return( result );
@ -498,6 +686,8 @@ vips_linear( VipsImage *in, VipsImage **out, double *a, double *b, int n, ... )
* Optional arguments:
*
* @uchar: output uchar pixels
* @a_imag: #VipsArrayDouble of imaginary constants for multiplication
* @b_imag: #VipsArrayDouble of imaginary constants for addition
*
* Run vips_linear() with a single constant.
*
@ -512,7 +702,76 @@ vips_linear1( VipsImage *in, VipsImage **out, double a, double b, ... )
int result;
va_start( ap, b );
result = vips_linearv( in, out, &a, &b, 1, ap );
result = vips_linearv( in, out, &a, NULL, &b, NULL, 1, ap );
va_end( ap );
return( result );
}
/**
* vips_linear_complex:
* @in: image to transform
* @out: output image
* @a: (array length=n): array of real constants for multiplication
* @a_imag: (array length=n): array of imaginary constants for multiplication
* @b: (array length=n): array of real constants for addition
* @b_imag: (array length=n): array of imaginary constants for addition
* @n: length of constant arrays
* @...: %NULL-terminated list of optional named arguments
*
* Optional arguments:
*
* @uchar: output uchar pixels
*
* Run vips_linear() with a set of complex constants.
*
* See also: vips_linear().
*
* Returns: 0 on success, -1 on error
*/
int
vips_linear_complex( VipsImage *in, VipsImage **out,
double *a, double *a_imag, double *b, double *b_imag, int n, ... )
{
va_list ap;
int result;
va_start( ap, n );
result = vips_linearv( in, out, a, a_imag, b, b_imag, n, ap );
va_end( ap );
return( result );
}
/**
* vips_linear_complex1:
* @in: image to transform
* @out: output image
* @a: real constant for multiplication
* @a_imag: imaginary constant for multiplication
* @b: real constant for addition
* @b_imag: imaginary constant for addition
* @...: %NULL-terminated list of optional named arguments
*
* Optional arguments:
*
* @uchar: output uchar pixels
*
* Run vips_linear() with a single complex constant.
*
* See also: vips_linear().
*
* Returns: 0 on success, -1 on error
*/
int
vips_linear_complex1( VipsImage *in, VipsImage **out,
double a, double a_imag, double b, double b_imag, ... )
{
va_list ap;
int result;
va_start( ap, b_imag );
result = vips_linearv( in, out, &a, &a_imag, &b, &b_imag, 1, ap );
va_end( ap );
return( result );

6
libvips/include/vips/arithmetic.h
View File

@ -185,6 +185,12 @@ int vips_linear( VipsImage *in, VipsImage **out,
__attribute__((sentinel));
int vips_linear1( VipsImage *in, VipsImage **out, double a, double b, ... )
__attribute__((sentinel));
int vips_linear_complex( VipsImage *in, VipsImage **out,
double *a, double *a_imag, double *b, double *b_imag, int n, ... )
__attribute__((sentinel));
int vips_linear_complex1( VipsImage *in, VipsImage **out,
double a, double a_imag, double b, double b_imag, ... )
__attribute__((sentinel));
int vips_remainder( VipsImage *left, VipsImage *right, VipsImage **out, ... )
__attribute__((sentinel));
int vips_remainder_const( VipsImage *in, VipsImage **out,