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

merge yafr branch back into trunk

This commit is contained in:
John Cupitt 2008-11-17 22:56:16 +00:00
parent dc68f4a3ec
commit b3dd71a55a
32 changed files with 2462 additions and 1621 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
ChangeLog
View File

@ -6,6 +6,11 @@
- added yafrsmooth interpolation
- added yafrtest
- added yafrnohalo
- ubuntu 8.10 changes
- interpolators get an output pointer, not region
- tuning for bicubic
- revised transform / clip code, we now do corner not centre
- yafr-smooth reworked along the lines of bicubic
11/9/08 started 7.16.3
- oop typo in manpage for im_project()

113
TODO
View File

@ -1,3 +1,116 @@
- make a new package for "resample"? im_shrink & friends could go in there too
- try walking the class hierarchy below VipsInterpolate and see if we can see
all the interpolators
classes need some extra fields
nickname, eg. "bicubic"
caption, eg. "catmull-rom bicubic interpolation"
- how to expose things like yafrsmooth's "sharpening" parameter to
nip2/C++/Python?
look at vips8, can we use the parametter code there?
would need to be added to VipsObject I guess
- repeat the benchmarks with the final code, how close are we? do we need to
repeat all of them because of the errors in the code we were testing?
user time, laptop, 8-bit data
bilinear 15.1 (was 16)
bicubic, fixed-point, tables 26.5 (was 29.5)
bicubic, floating-point, tables 35.8 (no change)
bicubic, floating-point, no tables 38.3 (was 44.2)
yafrsmooth, fixed-point, tables 114
yafrsmooth, from gegl 149
- we need an introspection thing to list interpolation classes
how do we list the subclasses of VipsInterpolate? it'll need some extra
stuff to present interpolations to the user (caption?)
we should add params as gobject properties
- benchmarking notes:
time vips im_affinei wtc.v wtc2.v 2 0.97 0 0 0.97 0 0 0 0 10000 10000
user time, opteron, 8-bit data
old-vips 3.0
bilinear 3.0
bicubic, fixed-point, separate tables 8.3
bicubic, fixed-point, tables 10.9
bicubic, floating-point, separate tables 12.8
bicubic, floating-point, tables 14.2
bicubic, floating-point, no tables, separate 15.0
bicubic, floating-point, no tables 17.6
yafr-smooth 36.3
user time, opteron, float data (real time is huge, thanks to very large file
read/write)
bicubic, fp, tables 8.2
bicubic, fp, tables, separate 9.7
bicubic, fp, no tables, separate 10.0
bicubic, fp, no tables 13.2
user time, opteron, 8 bit -> float, -> 8-bit
bicubic, fp, no tables 16
user time, laptop, 8-bit data
old-vips 14.7
bilinear 16.0
bicubic, fixed-point, tables 29.5
bicubic, floating-point, tables 35.8
bicubic, floating-point, no tables 44.2
user time, laptop, float data
bicubic, fp, tables 32.2
bicubic, fp, no tables 51.5
user time, laptop, 8 bit -> float, -> 8-bit
bicubic, fp, no tables 69
- still to check:
range clipping
think again about centred vs non-centred pixels
coordinate clipping
image margins
don't use tables for large expansions? we'll get banding
put bicubic's improvements into yafrsmooth / yafrnohalo
- think about position of . in int path
16 bit needs 16 for mantissa, plus 4 for adding the 16 components
together, plus a bit of overflow for safety, plus a sign bit
... so 8 bits for the fractional part sounds about right
for 8 bit data, we need 8 + 4 + 1 + 1, say 16 bits for the fractional
part
but better to keep it all with the same point position, or we need to
have two copies of the tables
- vips_object_print, set_name etc. need writing
- im_buf_t -> VipsBuf

1
include/vips/Makefile.am
View File

@ -3,6 +3,7 @@ pkginclude_HEADERS = \
VError.h \
VImage.h \
VMask.h \
bicubic.h \
vipscpp.h \
colour.h \
debug.h \

86
include/vips/bicubic.h Normal file
View File

@ -0,0 +1,86 @@
/* Bicubic (catmull-rom) interpolator.
*/
/*
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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
/*
These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
*/
#ifndef VIPS_BICUBIC_H
#define VIPS_BICUBIC_H
#ifdef __cplusplus
extern "C" {
#endif /*__cplusplus*/
#define VIPS_TYPE_INTERPOLATE_BICUBIC \
(vips_interpolate_bicubic_get_type())
#define VIPS_INTERPOLATE_BICUBIC( obj ) \
(G_TYPE_CHECK_INSTANCE_CAST( (obj), \
VIPS_TYPE_INTERPOLATE_BICUBIC, VipsInterpolateBicubic ))
#define VIPS_INTERPOLATE_BICUBIC_CLASS( klass ) \
(G_TYPE_CHECK_CLASS_CAST( (klass), \
VIPS_TYPE_INTERPOLATE_BICUBIC, VipsInterpolateBicubicClass))
#define VIPS_IS_INTERPOLATE_BICUBIC( obj ) \
(G_TYPE_CHECK_INSTANCE_TYPE( (obj), VIPS_TYPE_INTERPOLATE_BICUBIC ))
#define VIPS_IS_INTERPOLATE_BICUBIC_CLASS( klass ) \
(G_TYPE_CHECK_CLASS_TYPE( (klass), VIPS_TYPE_INTERPOLATE_BICUBIC ))
#define VIPS_INTERPOLATE_BICUBIC_GET_CLASS( obj ) \
(G_TYPE_INSTANCE_GET_CLASS( (obj), \
VIPS_TYPE_INTERPOLATE_BICUBIC, VipsInterpolateBicubicClass ))
typedef struct _VipsInterpolateBicubic {
VipsInterpolate parent_object;
} VipsInterpolateBicubic;
typedef struct _VipsInterpolateBicubicClass {
VipsInterpolateClass parent_class;
/* Precalculated interpolation matricies. int (used for pel sizes up
* to short), and double (for all others). We go to scale + 1, so
* we can round-to-nearest safely.
*/
/* We could keep a large set of 2d 4x4 matricies, but this actually
* works out slower, since for many resizes the thing will no longer
* fit in L1.
*/
int matrixi[VIPS_TRANSFORM_SCALE + 1][4];
double matrixf[VIPS_TRANSFORM_SCALE + 1][4];
} VipsInterpolateBicubicClass;
GType vips_interpolate_bicubic_get_type( void );
VipsInterpolate *vips_interpolate_bicubic_new( void );
/* Convenience: return a static default bicubic, so no need to free it.
*/
VipsInterpolate *vips_interpolate_bicubic_static( void );
#ifdef __cplusplus
}
#endif /*__cplusplus*/
#endif /*VIPS_BICUBIC_H*/

66
include/vips/interpolate.h
View File

@ -57,11 +57,11 @@ typedef struct _VipsInterpolate {
} VipsInterpolate;
/* An interpolation function. This is a class method, but we have a lookup
* function for it to speed up dispatch.
* function for it to speed up dispatch. Write to the memory at "out",
* interpolate the value at position (x, y) in "in".
*/
typedef void (*VipsInterpolateMethod)( VipsInterpolate *,
REGION *out, REGION *in,
int out_x, int out_y, double in_x, double in_y );
PEL *out, REGION *in, double x, double y );
typedef struct _VipsInterpolateClass {
VipsObjectClass parent_class;
@ -69,7 +69,7 @@ typedef struct _VipsInterpolateClass {
/* Write to pixel out(x,y), interpolating from in(x,y). The caller has
* to set the regions up.
*/
VipsInterpolateMethod interpolate;
VipsInterpolateMethod interpolate;
/* This interpolator needs a window this many pixels across and down.
*/
@ -81,8 +81,8 @@ typedef struct _VipsInterpolateClass {
} VipsInterpolateClass;
GType vips_interpolate_get_type( void );
void vips_interpolate( VipsInterpolate *interpolate, REGION *out, REGION *in,
int out_x, int out_y, double in_x, double in_y );
void vips_interpolate( VipsInterpolate *interpolate,
PEL *out, REGION *in, double x, double y );
VipsInterpolateMethod vips_interpolate_get_method( VipsInterpolate * );
int vips_interpolate_get_window_size( VipsInterpolate *interpolate );
@ -131,9 +131,11 @@ VipsInterpolate *vips_interpolate_nearest_static( void );
#define VIPS_TRANSFORM_SCALE (1 << VIPS_TRANSFORM_SHIFT)
/* How many bits of precision we keep for interpolation, ie. where the decimal
* is in the fixed-point tables.
* is in the fixed-point tables. For 16-bit pixels, we need 16 bits for the
* data, 4 bits to add 16 values together, another bit for the sign and some
* other stuff, so say 24 total. That leaves 8 bits for the fractional part.
*/
#define VIPS_INTERPOLATE_SHIFT (13)
#define VIPS_INTERPOLATE_SHIFT (8)
#define VIPS_INTERPOLATE_SCALE (1 << VIPS_INTERPOLATE_SHIFT)
#define VIPS_TYPE_INTERPOLATE_BILINEAR (vips_interpolate_bilinear_get_type())
@ -160,59 +162,21 @@ typedef struct _VipsInterpolateBilinearClass {
VipsInterpolateClass parent_class;
/* Precalculated interpolation matricies. int (used for pel sizes up
* to short), and double (for all others). We go to scale + 1, so
* we can round-to-nearest safely.
* to short), and float (for all others). We go to scale + 1, so
* we can round-to-nearest safely. Don't bother with double, since
* this is an approximation anyway.
*/
int matrixi[VIPS_TRANSFORM_SCALE + 1][VIPS_TRANSFORM_SCALE + 1][4];
double matrixd[VIPS_TRANSFORM_SCALE + 1][VIPS_TRANSFORM_SCALE + 1][4];
float matrixd[VIPS_TRANSFORM_SCALE + 1][VIPS_TRANSFORM_SCALE + 1][4];
} VipsInterpolateBilinearClass;
GType vips_interpolate_bilinear_get_type( void );
VipsInterpolate *vips_interpolate_bilinear_new( void );
/* Convenience: return a static fast bilinear, so no need to free it.
/* Convenience: return a static bilinear, so no need to free it.
*/
VipsInterpolate *vips_interpolate_bilinear_static( void );
/* Slow bilinear class starts.
*/
#define VIPS_TYPE_INTERPOLATE_BILINEAR_SLOW \
(vips_interpolate_bilinear_slow_get_type())
#define VIPS_INTERPOLATE_BILINEAR_SLOW( obj ) \
(G_TYPE_CHECK_INSTANCE_CAST( (obj), \
VIPS_TYPE_INTERPOLATE_BILINEAR_SLOW, VipsInterpolateBilinearSlow ))
#define VIPS_INTERPOLATE_BILINEAR_SLOW_CLASS( klass ) \
(G_TYPE_CHECK_CLASS_CAST( (klass), \
VIPS_TYPE_INTERPOLATE_BILINEAR_SLOW, \
VipsInterpolateBilinearSlowClass))
#define VIPS_IS_INTERPOLATE_BILINEAR_SLOW( obj ) \
(G_TYPE_CHECK_INSTANCE_TYPE( (obj), \
VIPS_TYPE_INTERPOLATE_BILINEAR_SLOW ))
#define VIPS_IS_INTERPOLATE_BILINEAR_SLOW_CLASS( klass ) \
(G_TYPE_CHECK_CLASS_TYPE( (klass), \
VIPS_TYPE_INTERPOLATE_BILINEAR_SLOW ))
#define VIPS_INTERPOLATE_BILINEAR_SLOW_GET_CLASS( obj ) \
(G_TYPE_INSTANCE_GET_CLASS( (obj), \
VIPS_TYPE_INTERPOLATE_BILINEAR_SLOW, VipsInterpolateBilinearSlowClass ))
typedef struct _VipsInterpolateBilinearSlow {
VipsInterpolate parent_object;
} VipsInterpolateBilinearSlow;
typedef struct _VipsInterpolateBilinearSlowClass {
VipsInterpolateClass parent_class;
} VipsInterpolateBilinearSlowClass;
GType vips_interpolate_bilinear_slow_get_type( void );
VipsInterpolate *vips_interpolate_bilinear_slow_new( void );
/* Convenience: return a static fast bilinear_slow, so no need to free it.
*/
VipsInterpolate *vips_interpolate_bilinear_slow_static( void );
#ifdef __cplusplus
}
#endif /*__cplusplus*/

27
include/vips/proto.h
View File

@ -52,7 +52,6 @@ extern "C" {
/* Need these for some protos.
*/
#include <stdio.h>
#include <stdarg.h>
#include <sys/types.h>
#include <glib-object.h>
@ -65,10 +64,6 @@ extern "C" {
# endif
#endif /*SWIG*/
typedef int (*im_callback_fn)( void *, void * );
typedef void *(*im_construct_fn)( void *, void *, void * );
typedef void *(*im_header_map_fn)( IMAGE *, const char *, GValue *, void * );
/* iofuncs
*/
int im_init_world( const char *argv0 );
@ -87,6 +82,7 @@ int im_header_double( IMAGE *im, const char *field, double *out );
int im_header_string( IMAGE *im, const char *field, char **out );
GType im_header_get_type( IMAGE *im, const char *field );
int im_header_get( IMAGE *im, const char *field, GValue *value_copy );
typedef void *(*im_header_map_fn)( IMAGE *, const char *, GValue *, void * );
void *im_header_map( IMAGE *im, im_header_map_fn fn, void *a );
const char *im_version_string( void );
@ -415,6 +411,8 @@ typedef enum {
IM_ARCH_MSB_FIRST
} im_arch_type;
gboolean im_isnative( im_arch_type arch );
DOUBLEMASK *im_vips2mask( IMAGE *, const char * );
int im_mask2vips( DOUBLEMASK *, IMAGE * );
int im_copy_set( IMAGE *, IMAGE *, int, float, float, int, int );
@ -637,6 +635,19 @@ int im_match_linear_search( IMAGE *ref, IMAGE *sec, IMAGE *out,
int xr2, int yr2, int xs2, int ys2,
int hwindowsize, int hsearchsize );
int im_affinei( IMAGE *in, IMAGE *out,
VipsInterpolate *interpolate,
double a, double b, double c, double d, double dx, double dy,
int ox, int oy, int ow, int oh );
int im_correl( IMAGE *ref, IMAGE *sec,
int xref, int yref, int xsec, int ysec,
int hwindowsize, int hsearchsize,
double *correlation, int *x, int *y );
int im_remosaic( IMAGE *in, IMAGE *out,
const char *old_str, const char *new_str );
/* Old stuff, for compat.
*/
int im_affine( IMAGE *in, IMAGE *out,
double a, double b, double c, double d, double dx, double dy,
int ox, int oy, int ow, int oh );
@ -645,12 +656,6 @@ int im_similarity( IMAGE *in, IMAGE *out,
int im_similarity_area( IMAGE *in, IMAGE *out,
double a, double b, double dx, double dy,
int ox, int oy, int ow, int oh );
int im_correl( IMAGE *ref, IMAGE *sec,
int xref, int yref, int xsec, int ysec,
int hwindowsize, int hsearchsize,
double *correlation, int *x, int *y );
int im_remosaic( IMAGE *in, IMAGE *out,
const char *old_str, const char *new_str );
int im_align_bands( IMAGE *in, IMAGE *out );
int im_maxpos_subpel( IMAGE *in, double *x, double *y );

9
include/vips/util.h
View File

@ -38,6 +38,8 @@
extern "C" {
#endif /*__cplusplus*/
#include <stdio.h>
/* Some platforms don't have M_PI in math.h :-(
*/
#define IM_PI (3.14159265358979323846)
@ -269,6 +271,11 @@ extern "C" {
(im_construct_fn) im_open, (im_callback_fn) im_close, \
(void *) (NAME), (void *) (MODE), NULL ))
/* Basic function types.
*/
typedef int (*im_callback_fn)( void *, void * );
typedef void *(*im_construct_fn)( void *, void *, void * );
/* strtok replacement.
*/
char *im__break_token( char *str, char *brk );
@ -322,8 +329,6 @@ char *im__file_read( FILE *fp, const char *name, unsigned int *length_out );
char *im__file_read_name( const char *name, unsigned int *length_out );
int im__file_write( void *data, size_t size, size_t nmemb, FILE *stream );
gboolean im_isnative( im_arch_type arch );
#ifdef __cplusplus
}
#endif /*__cplusplus*/

5
include/vips/vips.h
View File

@ -494,7 +494,7 @@ typedef struct {
(X) * IM_IMAGE_SIZEOF_PEL(I))
#endif /*DEBUG*/
#include <vips/proto.h>
#include <vips/util.h>
#include <vips/colour.h>
/* #include <vips/vector.h> */
#include <vips/format.h>
@ -504,10 +504,11 @@ typedef struct {
#include <vips/yafrsmooth.h>
#include <vips/yafrnohalo.h>
#include <vips/yafrtest.h>
#include <vips/bicubic.h>
#include <vips/semaphore.h>
#include <vips/threadgroup.h>
#include <vips/meta.h>
#include <vips/util.h>
#include <vips/proto.h>
#ifdef __cplusplus
}

13
include/vips/yafrsmooth.h
View File

@ -1,4 +1,4 @@
/* YAFRSMOOTH interpolator.
/* Yafrsmooth (catmull-rom) interpolator.
*/
/*
@ -88,6 +88,17 @@ typedef struct _VipsInterpolateYafrsmooth {
typedef struct _VipsInterpolateYafrsmoothClass {
VipsInterpolateClass parent_class;
/* Precalculated interpolation matricies. int (used for pel sizes up
* to short), and double (for all others). We go to scale + 1, so
* we can round-to-nearest safely.
*/
/* We could keep a large set of 2d 4x4 matricies, but this actually
* works out slower, since for many resizes the thing will no longer
* fit in L1.
*/
int matrixi[VIPS_TRANSFORM_SCALE + 1][4];
double matrixf[VIPS_TRANSFORM_SCALE + 1][4];
} VipsInterpolateYafrsmoothClass;
GType vips_interpolate_yafrsmooth_get_type( void );

7
libsrc/mosaicing/Makefile.am
View File

@ -1,9 +1,9 @@
noinst_LTLIBRARIES = libmosaicing.la
libmosaicing_la_SOURCES = \
im_affinei.c \
im_affine.c \
im_align_bands.c \
bicubic.cpp \
match.c \
mosaic1.c \
mosaicing_dispatch.c \
@ -13,7 +13,7 @@ libmosaicing_la_SOURCES = \
im_chkpair.c \
im_clinear.c \
interpolate.c \
yafrsmooth.c \
yafrsmooth.cpp \
yafrnohalo.c \
yafrtest.cpp \
im_improve.c \
@ -26,6 +26,9 @@ libmosaicing_la_SOURCES = \
im_tbmerge.c \
im_remosaic.c \
im_tbmosaic.c \
templates.h \
transform.h \
transform.c \
merge.h \
global_balance.h \
mosaic.h

459
libsrc/mosaicing/bicubic.cpp Normal file
View File

@ -0,0 +1,459 @@
/* bicubic ... catmull-rom interpolator
*/
/*
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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
/*
These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
*/
/* Bicubic (catmull-rom) interpolator derived from Nicolas Robidoux's YAFR
* resampler with permission and thanks.
*/
/*
#define DEBUG
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif /*HAVE_CONFIG_H*/
#include <vips/intl.h>
#include <stdio.h>
#include <stdlib.h>
#include <vips/vips.h>
#include <vips/internal.h>
#include "templates.h"
#ifdef WITH_DMALLOC
#include <dmalloc.h>
#endif /*WITH_DMALLOC*/
static VipsInterpolateClass *vips_interpolate_bicubic_parent_class = NULL;
/* Pointers to write to / read from, number of bands,
* how many bytes to add to move down a line.
*/
/* T is the type of pixels we are reading and writing.
*/
/* Fixed-point version, for 8 and 16-bit types.
*/
template <typename T, int min_value, int max_value>
static void inline
bicubic_int_tab( PEL *pout, const PEL *pin,
const int bands, const int lskip,
const int *cx, const int *cy )
{
T* restrict out = (T *) pout;
const T* restrict in = (T *) pin;
const int b1 = bands;
const int b2 = b1 + b1;
const int b3 = b1 + b2;
const int l1 = lskip / sizeof( T );
const int l2 = l1 + l1;
const int l3 = l1 + l2;
const int l1_plus_b1 = l1 + b1;
const int l1_plus_b2 = l1 + b2;
const int l1_plus_b3 = l1 + b3;
const int l2_plus_b1 = l2 + b1;
const int l2_plus_b2 = l2 + b2;
const int l2_plus_b3 = l2 + b3;
const int l3_plus_b1 = l3 + b1;
const int l3_plus_b2 = l3 + b2;
const int l3_plus_b3 = l3 + b3;
for( int z = 0; z < bands; z++ ) {
const T uno_one = in[0];
const T uno_two = in[b1];
const T uno_thr = in[b2];
const T uno_fou = in[b3];
const T dos_one = in[l1];
const T dos_two = in[l1_plus_b1];
const T dos_thr = in[l1_plus_b2];
const T dos_fou = in[l1_plus_b3];
const T tre_one = in[l2];
const T tre_two = in[l2_plus_b1];
const T tre_thr = in[l2_plus_b2];
const T tre_fou = in[l2_plus_b3];
const T qua_one = in[l3];
const T qua_two = in[l3_plus_b1];
const T qua_thr = in[l3_plus_b2];
const T qua_fou = in[l3_plus_b3];
int bicubic = bicubic_int<T>(
uno_one, uno_two, uno_thr, uno_fou,
dos_one, dos_two, dos_thr, dos_fou,
tre_one, tre_two, tre_thr, tre_fou,
qua_one, qua_two, qua_thr, qua_fou,
cx, cy );
if( bicubic < min_value )
bicubic = min_value;
else if( bicubic > max_value )
bicubic = max_value;
out[z] = bicubic;
in += 1;
}
}
/* Floating-point version, for int/float types.
*/
template <typename T>
static void inline
bicubic_float_tab( PEL *pout, const PEL *pin,
const int bands, const int lskip,
const double *cx, const double *cy )
{
T* restrict out = (T *) pout;
const T* restrict in = (T *) pin;
const int b1 = bands;
const int b2 = b1 + b1;
const int b3 = b1 + b2;
const int l1 = lskip / sizeof( T );
const int l2 = l1 + l1;
const int l3 = l1 + l2;
const int l1_plus_b1 = l1 + b1;
const int l1_plus_b2 = l1 + b2;
const int l1_plus_b3 = l1 + b3;
const int l2_plus_b1 = l2 + b1;
const int l2_plus_b2 = l2 + b2;
const int l2_plus_b3 = l2 + b3;
const int l3_plus_b1 = l3 + b1;
const int l3_plus_b2 = l3 + b2;
const int l3_plus_b3 = l3 + b3;
for( int z = 0; z < bands; z++ ) {
const T uno_one = in[0];
const T uno_two = in[b1];
const T uno_thr = in[b2];
const T uno_fou = in[b3];
const T dos_one = in[l1];
const T dos_two = in[l1_plus_b1];
const T dos_thr = in[l1_plus_b2];
const T dos_fou = in[l1_plus_b3];
const T tre_one = in[l2];
const T tre_two = in[l2_plus_b1];
const T tre_thr = in[l2_plus_b2];
const T tre_fou = in[l2_plus_b3];
const T qua_one = in[l3];
const T qua_two = in[l3_plus_b1];
const T qua_thr = in[l3_plus_b2];
const T qua_fou = in[l3_plus_b3];
const T bicubic = bicubic_float<T>(
uno_one, uno_two, uno_thr, uno_fou,
dos_one, dos_two, dos_thr, dos_fou,
tre_one, tre_two, tre_thr, tre_fou,
qua_one, qua_two, qua_thr, qua_fou,
cx, cy );
out[z] = bicubic;
in += 1;
}
}
/* Ultra-high-quality version for double images.
*/
template <typename T>
static void inline
bicubic_notab( PEL *pout, const PEL *pin,
const int bands, const int lskip,
double x, double y )
{
T* restrict out = (T *) pout;
const T* restrict in = (T *) pin;
const int b1 = bands;
const int b2 = b1 + b1;
const int b3 = b1 + b2;
const int l1 = lskip / sizeof( T );
const int l2 = l1 + l1;
const int l3 = l1 + l2;
const int l1_plus_b1 = l1 + b1;
const int l1_plus_b2 = l1 + b2;
const int l1_plus_b3 = l1 + b3;
const int l2_plus_b1 = l2 + b1;
const int l2_plus_b2 = l2 + b2;
const int l2_plus_b3 = l2 + b3;
const int l3_plus_b1 = l3 + b1;
const int l3_plus_b2 = l3 + b2;
const int l3_plus_b3 = l3 + b3;
double cx[4];
double cy[4];
calculate_coefficients_catmull( x, cx );
calculate_coefficients_catmull( y, cy );
for( int z = 0; z < bands; z++ ) {
const T uno_one = in[0];
const T uno_two = in[b1];
const T uno_thr = in[b2];
const T uno_fou = in[b3];
const T dos_one = in[l1];
const T dos_two = in[l1_plus_b1];
const T dos_thr = in[l1_plus_b2];
const T dos_fou = in[l1_plus_b3];
const T tre_one = in[l2];
const T tre_two = in[l2_plus_b1];
const T tre_thr = in[l2_plus_b2];
const T tre_fou = in[l2_plus_b3];
const T qua_one = in[l3];
const T qua_two = in[l3_plus_b1];
const T qua_thr = in[l3_plus_b2];
const T qua_fou = in[l3_plus_b3];
const T bicubic = bicubic_float<T>(
uno_one, uno_two, uno_thr, uno_fou,
dos_one, dos_two, dos_thr, dos_fou,
tre_one, tre_two, tre_thr, tre_fou,
qua_one, qua_two, qua_thr, qua_fou,
cx, cy );
out[z] = bicubic;
in += 1;
}
}
static void
vips_interpolate_bicubic_interpolate( VipsInterpolate *interpolate,
PEL *out, REGION *in, double x, double y )
{
VipsInterpolateBicubicClass *bicubic_class =
VIPS_INTERPOLATE_BICUBIC_GET_CLASS( interpolate );
/* Scaled int.
*/
const double sx = x * VIPS_TRANSFORM_SCALE;
const double sy = y * VIPS_TRANSFORM_SCALE;
const int sxi = FLOOR( sx );
const int syi = FLOOR( sy );
/* Get index into interpolation table and unscaled integer
* position.
*/
const int tx = sxi & (VIPS_TRANSFORM_SCALE - 1);
const int ty = syi & (VIPS_TRANSFORM_SCALE - 1);
const int ix = sxi >> VIPS_TRANSFORM_SHIFT;
const int iy = syi >> VIPS_TRANSFORM_SHIFT;
/* Look up the tables we need.
*/
const int *cxi = bicubic_class->matrixi[tx];
const int *cyi = bicubic_class->matrixi[ty];
const double *cxf = bicubic_class->matrixf[tx];
const double *cyf = bicubic_class->matrixf[ty];
/* Back and up one to get the top-left of the 4x4.
*/
const PEL *p = (PEL *) IM_REGION_ADDR( in, ix - 1, iy - 1 );
/* Pel size and line size.
*/
const int bands = in->im->Bands;
const int lskip = IM_REGION_LSKIP( in );
#ifdef DEBUG
printf( "vips_interpolate_bicubic_interpolate: %g %g\n", x, y );
printf( "\tleft=%d, top=%d, width=%d, height=%d\n",
ix - 1, iy - 1, 4, 4 );
#endif /*DEBUG*/
switch( in->im->BandFmt ) {
case IM_BANDFMT_UCHAR:
bicubic_int_tab<unsigned char, 0, UCHAR_MAX>(
out, p, bands, lskip,
cxi, cyi );
/*
Handy for benchmarking
bicubic_float_tab<unsigned char>(
out, p, bands, lskip,
cxf, cyf );
bicubic_notab<unsigned char>(
out, p, bands, lskip,
x - ix, y - iy );
*/
break;
case IM_BANDFMT_CHAR:
bicubic_int_tab<signed char, SCHAR_MIN, SCHAR_MAX>(
out, p, bands, lskip,
cxi, cyi );
break;
case IM_BANDFMT_USHORT:
bicubic_int_tab<unsigned short, 0, USHRT_MAX>(
out, p, bands, lskip,
cxi, cyi );
break;
case IM_BANDFMT_SHORT:
bicubic_int_tab<signed short, SHRT_MIN, SHRT_MAX>(
out, p, bands, lskip,
cxi, cyi );
break;
case IM_BANDFMT_UINT:
bicubic_float_tab<unsigned int>( out, p, bands, lskip,
cxf, cyf );
break;
case IM_BANDFMT_INT:
bicubic_float_tab<signed int>( out, p, bands, lskip,
cxf, cyf );
break;
case IM_BANDFMT_FLOAT:
bicubic_float_tab<float>( out, p, bands, lskip,
cxf, cyf );
break;
case IM_BANDFMT_DOUBLE:
bicubic_notab<double>( out, p, bands, lskip,
x - ix, y - iy );
break;
case IM_BANDFMT_COMPLEX:
bicubic_float_tab<float>( out, p, bands * 2, lskip,
cxf, cyf );
break;
case IM_BANDFMT_DPCOMPLEX:
bicubic_notab<double>( out, p, bands * 2, lskip,
x - ix, y - iy );
break;
default:
break;
}
}
static void
vips_interpolate_bicubic_class_init( VipsInterpolateBicubicClass *iclass )
{
VipsInterpolateClass *interpolate_class =
VIPS_INTERPOLATE_CLASS( iclass );
vips_interpolate_bicubic_parent_class =
VIPS_INTERPOLATE_CLASS( g_type_class_peek_parent( iclass ) );
interpolate_class->interpolate = vips_interpolate_bicubic_interpolate;
interpolate_class->window_size = 4;
/* Build the tables of pre-computed coefficients.
*/
for( int x = 0; x < VIPS_TRANSFORM_SCALE + 1; x++ ) {
calculate_coefficients_catmull(
(float) x / VIPS_TRANSFORM_SCALE,
iclass->matrixf[x] );
for( int i = 0; i < 4; i++ )
iclass->matrixi[x][i] =
iclass->matrixf[x][i] * VIPS_INTERPOLATE_SCALE;
}
}
static void
vips_interpolate_bicubic_init( VipsInterpolateBicubic *bicubic )
{
#ifdef DEBUG
printf( "vips_interpolate_bicubic_init: " );
vips_object_print( VIPS_OBJECT( bicubic ) );
#endif /*DEBUG*/
}
GType
vips_interpolate_bicubic_get_type()
{
static GType type = 0;
if( !type ) {
static const GTypeInfo info = {
sizeof( VipsInterpolateBicubicClass ),
NULL, /* base_init */
NULL, /* base_finalize */
(GClassInitFunc) vips_interpolate_bicubic_class_init,
NULL, /* class_finalize */
NULL, /* class_data */
sizeof( VipsInterpolateBicubic ),
32, /* n_preallocs */
(GInstanceInitFunc) vips_interpolate_bicubic_init,
};
type = g_type_register_static( VIPS_TYPE_INTERPOLATE,
"VipsInterpolateBicubic", &info,
(GTypeFlags) 0 );
}
return( type );
}
VipsInterpolate *
vips_interpolate_bicubic_new( void )
{
return( VIPS_INTERPOLATE( g_object_new(
VIPS_TYPE_INTERPOLATE_BICUBIC, NULL ) ) );
}
/* Convenience: return a static bicubic you don't need to free.
*/
VipsInterpolate *
vips_interpolate_bicubic_static( void )
{
static VipsInterpolate *interpolate = NULL;
if( !interpolate )
interpolate = vips_interpolate_bicubic_new();
return( interpolate );
}

1
libsrc/mosaicing/global_balance.c
View File

@ -104,6 +104,7 @@
#include <vips/vips.h>
#include "transform.h"
#include "merge.h"
#include "global_balance.h"

670
libsrc/mosaicing/im_affine.c
View File

@ -1,8 +1,9 @@
/* @(#) im_affine() ... affine transform, bi-linear interpolation.
/* @(#) im_affine() ... affine transform with a supplied interpolator.
* @(#)
* @(#) int im_affine(in, out, a, b, c, d, dx, dy, w, h, x, y)
* @(#) int im_affinei(in, out, interpolate, a, b, c, d, dx, dy, w, h, x, y)
* @(#)
* @(#) IMAGE *in, *out;
* @(#) VipsInterpolate *interpolate;
* @(#) double a, b, c, d, dx, dy;
* @(#) int w, h, x, y;
* @(#)
@ -73,6 +74,17 @@
* - still more tweaking, gah again
* 7/10/06
* - set THINSTRIP for no-rotate affines
* 20/10/08
* - version with interpolate parameter, from im_affine()
* 30/10/08
* - allow complex image types
* 4/11/08
* - take an interpolator as a param
* - replace im_affine with this, provide an im_affine() compat wrapper
* - break transform stuff out to transform.c
* - revise clipping / transform stuff, again
* - now do corner rather than centre: this way the identity transform
* returns the input exactly
*/
/*
@ -120,6 +132,7 @@
#include <vips/vips.h>
#include <vips/internal.h>
#include "transform.h"
#include "merge.h"
#ifdef WITH_DMALLOC
@ -130,343 +143,106 @@
*/
#define FLOOR( V ) ((V) >= 0 ? (int)(V) : (int)((V) - 1))
/* Precalculate a whole bunch of interpolation matricies. int (used for pel
* sizes up to short), and double (for all others). We go to scale + 1, so
* we can round-to-nearest safely.
FIXME ... should use seperable tables really
/* Per-call state.
*/
static int im_affine_linear_int
[TRANSFORM_SCALE + 1][TRANSFORM_SCALE + 1][4];
static double im_affine_linear_double
[TRANSFORM_SCALE + 1][TRANSFORM_SCALE + 1][4];
/* Make sure the interpolation tables are built.
*/
static void
affine_interpol_calc( void )
{
static int calced = 0;
int x, y;
if( calced )
return;
for( x = 0; x < TRANSFORM_SCALE + 1; x++ )
for( y = 0; y < TRANSFORM_SCALE + 1; y++ ) {
double X, Y, Xd, Yd;
double c1, c2, c3, c4;
/* Interpolation errors.
*/
X = (double) x / TRANSFORM_SCALE;
Y = (double) y / TRANSFORM_SCALE;
Xd = 1.0 - X;
Yd = 1.0 - Y;
/* Weights.
*/
c1 = Xd * Yd;
c2 = X * Yd;
c3 = X * Y;
c4 = Xd * Y;
im_affine_linear_double[x][y][0] = c1;
im_affine_linear_double[x][y][1] = c2;
im_affine_linear_double[x][y][2] = c3;
im_affine_linear_double[x][y][3] = c4;
im_affine_linear_int[x][y][0] = c1 * INTERPOL_SCALE;
im_affine_linear_int[x][y][1] = c2 * INTERPOL_SCALE;
im_affine_linear_int[x][y][2] = c3 * INTERPOL_SCALE;
im_affine_linear_int[x][y][3] = c4 * INTERPOL_SCALE;
}
calced = 1;
}
/* Calculate the inverse transformation.
*/
int
im__transform_calc_inverse( Transformation *trn )
{
DOUBLEMASK *msk, *msk2;
if( !(msk = im_create_dmaskv( "boink", 2, 2,
trn->a, trn->b, trn->c, trn->d )) )
return( -1 );
if( !(msk2 = im_matinv( msk, "boink2" )) ) {
(void) im_free_dmask( msk );
return( -1 );
}
trn->ia = msk2->coeff[0];
trn->ib = msk2->coeff[1];
trn->ic = msk2->coeff[2];
trn->id = msk2->coeff[3];
(void) im_free_dmask( msk );
(void) im_free_dmask( msk2 );
return( 0 );
}
/* Init a Transform.
*/
void
im__transform_init( Transformation *trn )
{
trn->oarea.left = 0;
trn->oarea.top = 0;
trn->oarea.width = -1;
trn->oarea.height = -1;
trn->iarea.left = 0;
trn->iarea.top = 0;
trn->iarea.width = -1;
trn->iarea.height = -1;
trn->a = 1.0; /* Identity transform */
trn->b = 0.0;
trn->c = 0.0;
trn->d = 1.0;
trn->dx = 0.0;
trn->dy = 0.0;
(void) im__transform_calc_inverse( trn );
}
/* Test for transform is identity function.
*/
int
im__transform_isidentity( Transformation *trn )
{
if( trn->a == 1.0 && trn->b == 0.0 && trn->c == 0.0 &&
trn->d == 1.0 && trn->dx == 0.0 && trn->dy == 0.0 )
return( 1 );
else
return( 0 );
}
/* Map a pixel coordinate through the transform.
*/
void
im__transform_forward( Transformation *trn,
double x, double y, /* In input space */
double *ox, double *oy ) /* In output space */
{
*ox = trn->a * x + trn->b * y + trn->dx;
*oy = trn->c * x + trn->d * y + trn->dy;
}
/* Map a pixel coordinate through the inverse transform.
*/
void
im__transform_inverse( Transformation *trn,
double x, double y, /* In output space */
double *ox, double *oy ) /* In input space */
{
double mx = x - trn->dx;
double my = y - trn->dy;
*ox = trn->ia * mx + trn->ib * my;
*oy = trn->ic * mx + trn->id * my;
}
/* Combine two transformations. out can be one of the ins.
*/
int
im__transform_add( Transformation *in1, Transformation *in2,
Transformation *out )
{
out->a = in1->a * in2->a + in1->c * in2->b;
out->b = in1->b * in2->a + in1->d * in2->b;
out->c = in1->a * in2->c + in1->c * in2->d;
out->d = in1->b * in2->c + in1->d * in2->d;
out->dx = in1->dx * in2->a + in1->dy * in2->b + in2->dx;
out->dy = in1->dx * in2->c + in1->dy * in2->d + in2->dy;
if( im__transform_calc_inverse( out ) )
return( -1 );
return( 0 );
}
void
im__transform_print( Transformation *trn )
{
printf( "im__transform_print:\n" );
printf( " iarea: left=%d, top=%d, width=%d, height=%d\n",
trn->iarea.left,
trn->iarea.top,
trn->iarea.width,
trn->iarea.height );
printf( " oarea: left=%d, top=%d, width=%d, height=%d\n",
trn->oarea.left,
trn->oarea.top,
trn->oarea.width,
trn->oarea.height );
printf( " mat: a=%g, b=%g, c=%g, d=%g\n",
trn->a, trn->b, trn->c, trn->d );
printf( " off: dx=%g, dy=%g\n",
trn->dx, trn->dy );
}
/* Map a point through the inverse transform. Used for clipping calculations,
* so it takes account of iarea and oarea.
*/
static void
invert_point( Transformation *trn,
double x, double y, /* In output space */
double *ox, double *oy ) /* In input space */
{
double xin = x - trn->oarea.left - trn->dx;
double yin = y - trn->oarea.top - trn->dy;
/* Find the inverse transform of current (x, y)
*/
*ox = trn->ia * xin + trn->ib * yin;
*oy = trn->ic * xin + trn->id * yin;
}
/* Given a bounding box for an area in the output image, set the bounding box
* for the corresponding pixels in the input image.
*/
static void
invert_rect( Transformation *trn,
Rect *in, /* In output space */
Rect *out ) /* In input space */
{
double x1, y1; /* Map corners */
double x2, y2;
double x3, y3;
double x4, y4;
double left, right, top, bottom;
/* Map input Rect.
*/
invert_point( trn, in->left, in->top, &x1, &y1 );
invert_point( trn, in->left, IM_RECT_BOTTOM(in), &x2, &y2 );
invert_point( trn, IM_RECT_RIGHT(in), in->top, &x3, &y3 );
invert_point( trn, IM_RECT_RIGHT(in), IM_RECT_BOTTOM(in), &x4, &y4 );
/* Find bounding box for these four corners.
*/
left = IM_MIN( x1, IM_MIN( x2, IM_MIN( x3, x4 ) ) );
right = IM_MAX( x1, IM_MAX( x2, IM_MAX( x3, x4 ) ) );
top = IM_MIN( y1, IM_MIN( y2, IM_MIN( y3, y4 ) ) );
bottom = IM_MAX( y1, IM_MAX( y2, IM_MAX( y3, y4 ) ) );
/* Set output Rect.
*/
out->left = left;
out->top = top;
out->width = right - left + 1;
out->height = bottom - top + 1;
/* Add a border for interpolation. You'd think +1 would do it, but
* we need to allow for rounding clipping as well.
FIXME ... will need adjusting when we add bicubic
*/
im_rect_marginadjust( out, 2 );
}
/* Interpolate a section ... int8/16 types.
*/
#define DO_IPEL(TYPE) { \
TYPE *tq = (TYPE *) q; \
\
int c1 = im_affine_linear_int[xi][yi][0]; \
int c2 = im_affine_linear_int[xi][yi][1]; \
int c3 = im_affine_linear_int[xi][yi][2]; \
int c4 = im_affine_linear_int[xi][yi][3]; \
\
/* p1 points to location (x_int, y_int) \
* p2 " " " (x_int+1, y_int) \
* p4 " " " (x_int+1, y_int+1) \
* p3 " " " (x_int, y_int+1) \
*/ \
PEL *p1 = (PEL *) IM_REGION_ADDR( ir, x_int, y_int ); \
PEL *p2 = p1 + ofs2; \
PEL *p3 = p1 + ofs3; \
PEL *p4 = p1 + ofs4; \
TYPE *tp1 = (TYPE *) p1; \
TYPE *tp2 = (TYPE *) p2; \
TYPE *tp3 = (TYPE *) p3; \
TYPE *tp4 = (TYPE *) p4; \
\
/* Interpolate each band. \
*/ \
for( z = 0; z < in->Bands; z++ ) \
tq[z] = (c1*tp1[z] + c2*tp2[z] + \
c3*tp3[z] + c4*tp4[z]) >> INTERPOL_SHIFT; \
}
/* Interpolate a pel ... int32 and float types.
*/
#define DO_FPEL(TYPE) { \
TYPE *tq = (TYPE *) q; \
\
double c1 = im_affine_linear_double[xi][yi][0]; \
double c2 = im_affine_linear_double[xi][yi][1]; \
double c3 = im_affine_linear_double[xi][yi][2]; \
double c4 = im_affine_linear_double[xi][yi][3]; \
\
/* p1 points to location (x_int, y_int) \
* p2 " " " (x_int+1, y_int) \
* p4 " " " (x_int+1, y_int+1) \
* p3 " " " (x_int, y_int+1) \
*/ \
PEL *p1 = (PEL *) IM_REGION_ADDR( ir, x_int, y_int ); \
PEL *p2 = p1 + ofs2; \
PEL *p3 = p1 + ofs3; \
PEL *p4 = p1 + ofs4; \
TYPE *tp1 = (TYPE *) p1; \
TYPE *tp2 = (TYPE *) p2; \
TYPE *tp3 = (TYPE *) p3; \
TYPE *tp4 = (TYPE *) p4; \
\
/* Interpolate each band. \
*/ \
for( z = 0; z < in->Bands; z++ ) \
tq[z] = c1*tp1[z] + c2*tp2[z] + \
c3*tp3[z] + c4*tp4[z]; \
}
typedef struct _Affine {
IMAGE *in;
IMAGE *out;
VipsInterpolate *interpolate;
Transformation trn;
} Affine;
static int
affine_gen( REGION *or, void *seq, void *a, void *b )
affine_free( Affine *affine )
{
IM_FREEF( g_object_unref, affine->interpolate );
return( 0 );
}
/* We have five (!!) coordinate systems. Working forward through them, there
* are:
*
* 1. The original input image
*
* 2. This is embedded in a larger image to provide borders for the
* interpolator. iarea->left/top give the offset. These are the coordinates we
* pass to IM_REGION_ADDR()/im_prepare() for the input image.
*
* 3. We need point (0, 0) in (1) to be at (0, 0) for the transformation. So
* shift everything up and left to make the displaced input image. This is the
* space that the transformation maps from, and can have negative pixels
* (up and left of the image, for interpolation).
*
* 4. Output transform space. This is the where the transform maps to. Pixels
* can be negative, since a rotated image can go up and left of the origin.
*
* 5. Output image space. This is the wh of the xywh passed to im_affine()
* below. These are the coordinates we pass to IM_REGION_ADDR() for the
* output image, and that affinei_gen() is asked for.
*/
static int
affinei_gen( REGION *or, void *seq, void *a, void *b )
{
REGION *ir = (REGION *) seq;
IMAGE *in = (IMAGE *) a;
Transformation *trn = (Transformation *) b;
const IMAGE *in = (IMAGE *) a;
const Affine *affine = (Affine *) b;
const int window_size =
vips_interpolate_get_window_size( affine->interpolate );
const int half_window_size = window_size / 2;
const VipsInterpolateMethod interpolate =
vips_interpolate_get_method( affine->interpolate );
/* Output area for this call.
/* Area we generate in the output image.
*/
Rect *r = &or->valid;
int le = r->left;
int ri = IM_RECT_RIGHT(r);
int to = r->top;
int bo = IM_RECT_BOTTOM(r);
Rect *iarea = &trn->iarea;
Rect *oarea = &trn->oarea;
const Rect *r = &or->valid;
const int le = r->left;
const int ri = IM_RECT_RIGHT( r );
const int to = r->top;
const int bo = IM_RECT_BOTTOM( r );
const Rect *iarea = &affine->trn.iarea;
const Rect *oarea = &affine->trn.oarea;
int ps = IM_IMAGE_SIZEOF_PEL( in );
int x, y, z;
/* Interpolation variables.
*/
int ofs2, ofs3, ofs4;
Rect image, want, need, clipped;
/* Clipping Rects.
#ifdef DEBUG
printf( "affine: generating left=%d, top=%d, width=%d, height=%d\n",
r->left,
r->top,
r->width,
r->height );
#endif /*DEBUG*/
/* We are generating this chunk of the transformed image.
*/
Rect image, need, clipped;
want = *r;
want.left += oarea->left;
want.top += oarea->top;
/* Find the area of the input image we need.
*/
im__transform_invert_rect( &affine->trn, &want, &need );
/* Now go to space (2) above.
*/
need.left += iarea->left;
need.top += iarea->top;
/* Add a border for interpolation.
*/
im_rect_marginadjust( &need, half_window_size );
/* Clip against the size of (2).
*/
image.left = 0;
image.top = 0;
image.width = in->Xsize;
image.height = in->Ysize;
invert_rect( trn, r, &need );
im_rect_intersectrect( &need, &image, &clipped );
/* Outside input image? All black.
@ -490,127 +266,63 @@ affine_gen( REGION *or, void *seq, void *a, void *b )
clipped.height );
#endif /*DEBUG*/
/* Calculate pel offsets.
*/
ofs2 = IM_IMAGE_SIZEOF_PEL( in );
ofs3 = ofs2 + IM_REGION_LSKIP( ir );
ofs4 = IM_REGION_LSKIP( ir );
/* Resample!
/* Resample! x/y loop over pixels in the output image (5).
*/
for( y = to; y < bo; y++ ) {
/* Continuous cods in output space.
/* Input clipping rectangle.
*/
double oy = y - oarea->top - trn->dy;
double ox;
const int ile = iarea->left;
const int ito = iarea->top;
const int iri = iarea->left + iarea->width;
const int ibo = iarea->top + iarea->height;
/* Input clipping rectangle.
*/
int ile = iarea->left;
int ito = iarea->top;
int iri = iarea->left + iarea->width;
int ibo = iarea->top + iarea->height;
/* Derivative of matrix.
*/
double dx = trn->ia;
double dy = trn->ic;
const double ddx = affine->trn.ia;
const double ddy = affine->trn.ic;
/* Continuous cods in transformed space.
*/
const double ox = le + oarea->left - affine->trn.dx;
const double oy = y + oarea->top - affine->trn.dy;
/* Continuous cods in input space.
*/
double ix, iy;
PEL *q;
q = (PEL *) IM_REGION_ADDR( or, le, y );
ox = le - oarea->left - trn->dx;
ix = trn->ia * ox + trn->ib * oy;
iy = trn->ic * ox + trn->id * oy;
/* To (3).
*/
ix = affine->trn.ia * ox + affine->trn.ib * oy;
iy = affine->trn.ic * ox + affine->trn.id * oy;
/* Offset ix/iy input by iarea.left/top ... so we skip the
* image edges we added for interpolation.
/* Now move to (2).
*/
ix += iarea->left;
iy += iarea->top;
q = (PEL *) IM_REGION_ADDR( or, le, y );
for( x = le; x < ri; x++ ) {
int fx, fy;
fx = FLOOR( ix );
fy = FLOOR( iy );
/* Clipping! Use >= for right/bottom, since IPOL needs
* to see one pixel more each way.
/* Clipping!
*/
if( fx < ile || fx >= iri || fy < ito || fy >= ibo ) {
for( z = 0; z < ps; z++ )
q[z] = 0;
}
else {
double sx, sy;
int x_int, y_int;
int xi, yi;
/* Subtract 0.5 to centre the bilinear.
FIXME ... need to adjust for bicubic.
*/
sx = ix - 0.5;
sy = iy - 0.5;
/* Now go to scaled int.
*/
sx *= TRANSFORM_SCALE;
sy *= TRANSFORM_SCALE;
x_int = FLOOR( sx );
y_int = FLOOR( sy );
/* Get index into interpolation table and
* unscaled integer position.
*/
xi = x_int & (TRANSFORM_SCALE - 1);
yi = y_int & (TRANSFORM_SCALE - 1);
x_int = x_int >> TRANSFORM_SHIFT;
y_int = y_int >> TRANSFORM_SHIFT;
/* Interpolate for each input type.
*/
switch( in->BandFmt ) {
case IM_BANDFMT_UCHAR:
DO_IPEL( unsigned char );
break;
case IM_BANDFMT_CHAR:
DO_IPEL( char );
break;
case IM_BANDFMT_USHORT:
DO_IPEL( unsigned short );
break;
case IM_BANDFMT_SHORT:
DO_IPEL( short );
break;
case IM_BANDFMT_UINT:
DO_FPEL( unsigned int );
break;
case IM_BANDFMT_INT:
DO_FPEL( int );
break;
case IM_BANDFMT_FLOAT:
DO_FPEL( float );
break;
case IM_BANDFMT_DOUBLE:
DO_FPEL( double );
break;
default:
error_exit( "im_affine: panic!");
/*NOTREACHED*/
}
interpolate( affine->interpolate,
q, ir, ix, iy );
}
ix += dx;
iy += dy;
ix += ddx;
iy += ddy;
q += ps;
}
}
@ -619,41 +331,42 @@ affine_gen( REGION *or, void *seq, void *a, void *b )
}
static int
affine( IMAGE *in, IMAGE *out, Transformation *trn )
affinei( IMAGE *in, IMAGE *out,
VipsInterpolate *interpolate, Transformation *trn )
{
Transformation *trn2;
Affine *affine;
double edge;
if( im_iscomplex( in ) ) {
im_error( "im_affine",
"%s", _( "complex input not supported" ) );
return( -1 );
}
/* We output at (0,0), so displace output by that amount -ve to get
* output at (ox,oy). Alter our copy of trn.
*/
if( !(trn2 = IM_NEW( out, Transformation )) )
return( -1 );
*trn2 = *trn;
trn2->oarea.left = -trn->oarea.left;
trn2->oarea.top = -trn->oarea.top;
if( im__transform_calc_inverse( trn2 ) )
return( -1 );
/* Make output image.
*/
if( im_piocheck( in, out ) )
return( -1 );
if( im_cp_desc( out, in ) )
return( -1 );
out->Xsize = trn2->oarea.width;
out->Ysize = trn2->oarea.height;
/* Need a copy of the params for the lifetime of out.
*/
if( !(affine = IM_NEW( out, Affine )) )
return( -1 );
affine->interpolate = NULL;
if( im_add_close_callback( out,
(im_callback_fn) affine_free, affine, NULL ) )
return( -1 );
affine->in = in;
affine->out = out;
affine->interpolate = interpolate;
g_object_ref( interpolate );
affine->trn = *trn;
if( im__transform_calc_inverse( &affine->trn ) )
return( -1 );
out->Xsize = affine->trn.oarea.width;
out->Ysize = affine->trn.oarea.height;
/* Normally SMALLTILE ... except if this is a size up/down affine.
*/
if( trn->b == 0.0 && trn->c == 0.0 ) {
if( affine->trn.b == 0.0 && affine->trn.c == 0.0 ) {
if( im_demand_hint( out, IM_FATSTRIP, in, NULL ) )
return( -1 );
}
@ -665,11 +378,11 @@ affine( IMAGE *in, IMAGE *out, Transformation *trn )
/* Check for coordinate overflow ... we want to be able to hold the
* output space inside INT_MAX / TRANSFORM_SCALE.
*/
edge = INT_MAX / TRANSFORM_SCALE;
if( trn2->oarea.left < -edge || trn2->oarea.top < -edge ||
IM_RECT_RIGHT( &trn2->oarea ) > edge ||
IM_RECT_BOTTOM( &trn2->oarea ) > edge ) {
im_error( "im_affine",
edge = INT_MAX / VIPS_TRANSFORM_SCALE;
if( affine->trn.oarea.left < -edge || affine->trn.oarea.top < -edge ||
IM_RECT_RIGHT( &affine->trn.oarea ) > edge ||
IM_RECT_BOTTOM( &affine->trn.oarea ) > edge ) {
im_error( "im_affinei",
"%s", _( "output coordinates out of range" ) );
return( -1 );
}
@ -677,7 +390,7 @@ affine( IMAGE *in, IMAGE *out, Transformation *trn )
/* Generate!
*/
if( im_generate( out,
im_start_one, affine_gen, im_stop_one, in, trn2 ) )
im_start_one, affinei_gen, im_stop_one, in, affine ) )
return( -1 );
return( 0 );
@ -685,53 +398,48 @@ affine( IMAGE *in, IMAGE *out, Transformation *trn )
/* As above, but do IM_CODING_LABQ too. And embed the input.
*/
int
im__affine( IMAGE *in, IMAGE *out, Transformation *trn )
static int
im__affinei( IMAGE *in, IMAGE *out,
VipsInterpolate *interpolate, Transformation *trn )
{
IMAGE *t3 = im_open_local( out, "im_affine:3", "p" );
const int window_size = vips_interpolate_get_window_size( interpolate );
Transformation trn2;
#ifdef DEBUG_GEOMETRY
printf( "im__affine: %s\n", in->filename );
im__transform_print( trn );
#endif /*DEBUG_GEOMETRY*/
/* Add new pixels around the input so we can interpolate at the edges.
* Bilinear needs 0.5 pixels on all edges.
FIXME ... will need to fiddle with this when we add bicubic
*/
if( !t3 ||
im_embed( in, t3, 1,
1, 1, in->Xsize + 2, in->Ysize + 2 ) )
window_size / 2, window_size / 2,
in->Xsize + window_size, in->Ysize + window_size ) )
return( -1 );
/* Set iarea so we know what part of the input we can take.
*/
trn2 = *trn;
trn2.iarea.left += 1;
trn2.iarea.top += 1;
trn2.iarea.left += window_size / 2;
trn2.iarea.top += window_size / 2;
affine_interpol_calc();
#ifdef DEBUG_GEOMETRY
printf( "im__affinei: %s\n", in->filename );
im__transform_print( &trn2 );
#endif /*DEBUG_GEOMETRY*/
if( in->Coding == IM_CODING_LABQ ) {
IMAGE *t1 = im_open_local( out, "im_affine:1", "p" );
IMAGE *t2 = im_open_local( out, "im_affine:2", "p" );
IMAGE *t[2];
if( !t1 || !t2 ||
im_LabQ2LabS( t3, t1 ) ||
affine( t1, t2, &trn2 ) ||
im_LabS2LabQ( t2, out ) )
if( im_open_local_array( out, t, 2, "im_affine:2", "p" ) ||
im_LabQ2LabS( t3, t[0] ) ||
affinei( t[0], t[1], interpolate, &trn2 ) ||
im_LabS2LabQ( t[1], out ) )
return( -1 );
}
else if( in->Coding == IM_CODING_NONE ) {
if( affine( t3, out, &trn2 ) )
if( affinei( t3, out, interpolate, &trn2 ) )
return( -1 );
}
else {
im_error( "im_affine",
"%s", _( "unknown coding type" ) );
im_error( "im_affinei", "%s", _( "unknown coding type" ) );
return( -1 );
}
@ -744,21 +452,22 @@ im__affine( IMAGE *in, IMAGE *out, Transformation *trn )
}
int
im_affine( IMAGE *in, IMAGE *out,
double a, double b, double c, double d,
double dx, double dy,
im_affinei( IMAGE *in, IMAGE *out, VipsInterpolate *interpolate,
double a, double b, double c, double d, double dx, double dy,
int ox, int oy, int ow, int oh )
{
Transformation trn;
trn.iarea.left = 0;
trn.iarea.top = 0;
trn.iarea.width = in->Xsize;
trn.iarea.height = in->Ysize;
trn.oarea.left = ox;
trn.oarea.top = oy;
trn.oarea.width = ow;
trn.oarea.height = oh;
trn.iarea.left = 0;
trn.iarea.top = 0;
trn.iarea.width = in->Xsize;
trn.iarea.height = in->Ysize;
trn.a = a;
trn.b = b;
trn.c = c;
@ -766,5 +475,26 @@ im_affine( IMAGE *in, IMAGE *out,
trn.dx = dx;
trn.dy = dy;
return( im__affine( in, out, &trn ) );
return( im__affinei( in, out, interpolate, &trn ) );
}
/* Provide the old im__affine()/im_affine() as bilinear affinei.
*/
int
im__affine( IMAGE *in, IMAGE *out, Transformation *trn )
{
return( im__affinei( in, out,
vips_interpolate_bilinear_static(), trn ) );
}
int
im_affine( IMAGE *in, IMAGE *out,
double a, double b, double c, double d, double dx, double dy,
int ox, int oy, int ow, int oh )
{
return( im_affinei( in, out,
vips_interpolate_bicubic_static(),
a, b, c, d, dx, dy,
ox, oy, ow, oh ) );
}

501
libsrc/mosaicing/im_affinei.c
View File

@ -1,501 +0,0 @@
/* @(#) im_affine() ... affine transform with a supplied interpolator.
* @(#)
* @(#) int im_affinei(in, out, interpolate, a, b, c, d, dx, dy, w, h, x, y)
* @(#)
* @(#) IMAGE *in, *out;
* @(#) VipsInterpolate *interpolate;
* @(#) double a, b, c, d, dx, dy;
* @(#) int w, h, x, y;
* @(#)
* @(#) Forward transform
* @(#) X = a * x + b * y + dx
* @(#) Y = c * x + d * y + dy
* @(#)
* @(#) x and y are the coordinates in input image.
* @(#) X and Y are the coordinates in output image.
* @(#) (0,0) is the upper left corner.
*
* Copyright N. Dessipris
* Written on: 01/11/1991
* Modified on: 12/3/92 JC
* - rounding error in interpolation routine fixed
* - test for scale=1, angle=0 case fixed
* - clipping of output removed: redundant
* - various little tidies
* - problems remain with scale>20, size<10
*
* Re-written on: 20/08/92, J.Ph Laurent
*
* 21/02/93, JC
* - speed-ups
* - simplifications
* - im_similarity now calculates a window and calls this routine
* 6/7/93 JC
* - rewritten for partials
* - ANSIfied
* - now rotates any non-complex type
* 3/6/94 JC
* - C revised in bug search
* 9/6/94 JC
* - im_prepare() was preparing too small an area! oops
* 22/5/95 JC
* - added code to detect all-black output area case - helps lazy ip
* 3/7/95 JC
* - IM_CODING_LABQ handling moved to here
* 31/7/97 JC
* - dx/dy sign reversed to be less confusing ... now follows comment at
* top ... ax - by + dx etc.
* - tiny speed up, replaced the *++ on interpolation with [z]
* - im_similarity() moved in here
* - args swapped: was whxy, now xywh
* - didn't agree with dispatch fns before :(
* 3/3/98 JC
* - im_demand_hint() added
* 20/12/99 JC
* - im_affine() made from im_similarity_area()
* - transform stuff cleaned up a bit
* 14/4/01 JC
* - oops, invert_point() had a rounding problem
* 23/2/02 JC
* - pre-calculate interpolation matricies
* - integer interpolation for int8/16 types, double for
* int32/float/double
* - faster transformation
* 15/8/02 JC
* - records Xoffset/Yoffset
* 14/4/04
* - rounding, clipping and transforming revised, now pixel-perfect (or
* better than gimp, anyway)
* 22/6/05
* - all revised again, simpler and more reliable now
* 30/3/06
* - gah, still an occasional clipping problem
* 12/7/06
* - still more tweaking, gah again
* 7/10/06
* - set THINSTRIP for no-rotate affines
* 20/10/08
* - version with interpolate parameter, from im_affine()
*/
/*
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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
/*
These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
*/
/*
#define DEBUG
#define DEBUG_GEOMETRY
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif /*HAVE_CONFIG_H*/
#include <vips/intl.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <limits.h>
#include <vips/vips.h>
#include <vips/internal.h>
#include "merge.h"
#ifdef WITH_DMALLOC
#include <dmalloc.h>
#endif /*WITH_DMALLOC*/
/* "fast" floor() ... on my laptop, anyway.
*/
#define FLOOR( V ) ((V) >= 0 ? (int)(V) : (int)((V) - 1))
/* Map a point through the inverse transform. Used for clipping calculations,
* so it takes account of iarea and oarea.
*/
static void
invert_point( Transformation *trn,
double x, double y, /* In output space */
double *ox, double *oy ) /* In input space */
{
double xin = x - trn->oarea.left - trn->dx;
double yin = y - trn->oarea.top - trn->dy;
/* Find the inverse transform of current (x, y)
*/
*ox = trn->ia * xin + trn->ib * yin;
*oy = trn->ic * xin + trn->id * yin;
}
/* Given a bounding box for an area in the output image, set the bounding box
* for the corresponding pixels in the input image.
*/
static void
invert_rect( Transformation *trn,
Rect *in, /* In output space */
Rect *out ) /* In input space */
{
double x1, y1; /* Map corners */
double x2, y2;
double x3, y3;
double x4, y4;
double left, right, top, bottom;
/* Map input Rect.
*/
invert_point( trn, in->left, in->top, &x1, &y1 );
invert_point( trn, in->left, IM_RECT_BOTTOM(in), &x2, &y2 );
invert_point( trn, IM_RECT_RIGHT(in), in->top, &x3, &y3 );
invert_point( trn, IM_RECT_RIGHT(in), IM_RECT_BOTTOM(in), &x4, &y4 );
/* Find bounding box for these four corners.
*/
left = IM_MIN( x1, IM_MIN( x2, IM_MIN( x3, x4 ) ) );
right = IM_MAX( x1, IM_MAX( x2, IM_MAX( x3, x4 ) ) );
top = IM_MIN( y1, IM_MIN( y2, IM_MIN( y3, y4 ) ) );
bottom = IM_MAX( y1, IM_MAX( y2, IM_MAX( y3, y4 ) ) );
/* Set output Rect.
*/
out->left = left;
out->top = top;
out->width = right - left + 1;
out->height = bottom - top + 1;
}
/* Per-call state.
*/
typedef struct _Affine {
IMAGE *in;
IMAGE *out;
VipsInterpolate *interpolate;
Transformation trn;
} Affine;
static int
affine_free( Affine *affine )
{
IM_FREEF( g_object_unref, affine->interpolate );
return( 0 );
}
static int
affinei_gen( REGION *or, void *seq, void *a, void *b )
{
REGION *ir = (REGION *) seq;
IMAGE *in = (IMAGE *) a;
Affine *affine = (Affine *) b;
const int window_size =
vips_interpolate_get_window_size( affine->interpolate );
const int half_window_size = window_size / 2;
VipsInterpolateMethod interpolate =
vips_interpolate_get_method( affine->interpolate );
/* Output area for this call.
*/
Rect *r = &or->valid;
int le = r->left;
int ri = IM_RECT_RIGHT(r);
int to = r->top;
int bo = IM_RECT_BOTTOM(r);
Rect *iarea = &affine->trn.iarea;
Rect *oarea = &affine->trn.oarea;
int ps = IM_IMAGE_SIZEOF_PEL( in );
int x, y, z;
/* Clipping Rects.
*/
Rect image, need, clipped;
/* Find the area of the input image we need.
*/
image.left = 0;
image.top = 0;
image.width = in->Xsize;
image.height = in->Ysize;
invert_rect( &affine->trn, r, &need );
/* Add a border for interpolation. You'd think +1 would do it, but
* we need to allow for rounding clipping as well.
*/
im_rect_marginadjust( &need, window_size );
im_rect_intersectrect( &need, &image, &clipped );
/* Outside input image? All black.
*/
if( im_rect_isempty( &clipped ) ) {
im__black_region( or );
return( 0 );
}
/* We do need some pixels from the input image to make our output -
* ask for them.
*/
if( im_prepare( ir, &clipped ) )
return( -1 );
#ifdef DEBUG
printf( "affine: preparing left=%d, top=%d, width=%d, height=%d\n",
clipped.left,
clipped.top,
clipped.width,
clipped.height );
#endif /*DEBUG*/
/* Resample!
*/
for( y = to; y < bo; y++ ) {
/* Continuous cods in output space.
*/
double oy = y - oarea->top - affine->trn.dy;
double ox;
/* Input clipping rectangle.
*/
int ile = iarea->left;
int ito = iarea->top;
int iri = iarea->left + iarea->width;
int ibo = iarea->top + iarea->height;
/* Derivative of matrix.
*/
double dx = affine->trn.ia;
double dy = affine->trn.ic;
/* Continuous cods in input space.
*/
double ix, iy;
PEL *q;
ox = le - oarea->left - affine->trn.dx;
ix = affine->trn.ia * ox + affine->trn.ib * oy;
iy = affine->trn.ic * ox + affine->trn.id * oy;
/* Offset ix/iy input by iarea.left/top ... so we skip the
* image edges we added for interpolation.
*/
ix += iarea->left;
iy += iarea->top;
q = (PEL *) IM_REGION_ADDR( or, le, y );
for( x = le; x < ri; x++ ) {
int fx, fy;
fx = FLOOR( ix );
fy = FLOOR( iy );
/* Clipping! Use >= for right/bottom, since IPOL needs
* to see one pixel more each way.
*/
if( fx <= ile - half_window_size ||
fx >= iri + half_window_size ||
fy <= ito - half_window_size ||
fy >= ibo + half_window_size ) {
for( z = 0; z < ps; z++ )
q[z] = 0;
}
else {
interpolate( affine->interpolate,
or, ir,
x, y, ix, iy );
}
ix += dx;
iy += dy;
q += ps;
}
}
return( 0 );
}
static int
affinei( IMAGE *in, IMAGE *out,
VipsInterpolate *interpolate, Transformation *trn )
{
Affine *affine;
double edge;
if( im_iscomplex( in ) ) {
im_error( "im_affinei",
"%s", _( "complex input not supported" ) );
return( -1 );
}
/* Make output image.
*/
if( im_piocheck( in, out ) )
return( -1 );
if( im_cp_desc( out, in ) )
return( -1 );
/* Need a copy of the params for the lifetime of out.
*/
if( !(affine = IM_NEW( out, Affine )) )
return( -1 );
affine->interpolate = NULL;
if( im_add_close_callback( out,
(im_callback_fn) affine_free, affine, NULL ) )
return( -1 );
affine->in = in;
affine->out = out;
affine->interpolate = interpolate;
g_object_ref( interpolate );
affine->trn = *trn;
/* We output at (0,0), so displace output by that amount -ve to get
* output at (ox,oy). Alter our copy of trn.
*/
affine->trn.oarea.left = -affine->trn.oarea.left;
affine->trn.oarea.top = -affine->trn.oarea.top;
if( im__transform_calc_inverse( &affine->trn ) )
return( -1 );
out->Xsize = affine->trn.oarea.width;
out->Ysize = affine->trn.oarea.height;
/* Normally SMALLTILE ... except if this is a size up/down affine.
*/
if( affine->trn.b == 0.0 && affine->trn.c == 0.0 ) {
if( im_demand_hint( out, IM_FATSTRIP, in, NULL ) )
return( -1 );
}
else {
if( im_demand_hint( out, IM_SMALLTILE, in, NULL ) )
return( -1 );
}
/* Check for coordinate overflow ... we want to be able to hold the
* output space inside INT_MAX / TRANSFORM_SCALE.
*/
edge = INT_MAX / VIPS_TRANSFORM_SCALE;
if( affine->trn.oarea.left < -edge || affine->trn.oarea.top < -edge ||
IM_RECT_RIGHT( &affine->trn.oarea ) > edge ||
IM_RECT_BOTTOM( &affine->trn.oarea ) > edge ) {
im_error( "im_affinei",
"%s", _( "output coordinates out of range" ) );
return( -1 );
}
/* Generate!
*/
if( im_generate( out,
im_start_one, affinei_gen, im_stop_one, in, affine ) )
return( -1 );
return( 0 );
}
/* As above, but do IM_CODING_LABQ too. And embed the input.
*/
int
im__affinei( IMAGE *in, IMAGE *out,
VipsInterpolate *interpolate, Transformation *trn )
{
IMAGE *t3 = im_open_local( out, "im_affine:3", "p" );
const int window_size = vips_interpolate_get_window_size( interpolate );
Transformation trn2;
#ifdef DEBUG_GEOMETRY
printf( "im__affinei: %s\n", in->filename );
im__transform_print( trn );
#endif /*DEBUG_GEOMETRY*/
/* Add new pixels around the input so we can interpolate at the edges.
*/
if( !t3 ||
im_embed( in, t3, 1,
window_size / 2, window_size / 2,
in->Xsize + window_size, in->Ysize + window_size ) )
return( -1 );
/* Set iarea so we know what part of the input we can take.
*/
trn2 = *trn;
trn2.iarea.left += window_size / 2;
trn2.iarea.top += window_size / 2;
if( in->Coding == IM_CODING_LABQ ) {
IMAGE *t1 = im_open_local( out, "im_affine:1", "p" );
IMAGE *t2 = im_open_local( out, "im_affine:2", "p" );
if( !t1 || !t2 ||
im_LabQ2LabS( t3, t1 ) ||
affinei( t1, t2, interpolate, &trn2 ) ||
im_LabS2LabQ( t2, out ) )
return( -1 );
}
else if( in->Coding == IM_CODING_NONE ) {
if( affinei( t3, out, interpolate, &trn2 ) )
return( -1 );
}
else {
im_error( "im_affinei", "%s", _( "unknown coding type" ) );
return( -1 );
}
/* Finally: can now set Xoffset/Yoffset.
*/
out->Xoffset = trn->dx - trn->oarea.left;
out->Yoffset = trn->dy - trn->oarea.top;
return( 0 );
}
int
im_affinei( IMAGE *in, IMAGE *out, VipsInterpolate *interpolate,
double a, double b, double c, double d,
double dx, double dy,
int ox, int oy, int ow, int oh )
{
Transformation trn;
trn.iarea.left = 0;
trn.iarea.top = 0;
trn.iarea.width = in->Xsize;
trn.iarea.height = in->Ysize;
trn.oarea.left = ox;
trn.oarea.top = oy;
trn.oarea.width = ow;
trn.oarea.height = oh;
trn.a = a;
trn.b = b;
trn.c = c;
trn.d = d;
trn.dx = dx;
trn.dy = dy;
return( im__affinei( in, out, interpolate, &trn ) );
}

1
libsrc/mosaicing/im_lrmerge.c
View File

@ -134,6 +134,7 @@
#include <vips/vips.h>
#include <vips/thread.h>
#include "transform.h"
#include "merge.h"
#ifdef WITH_DMALLOC

1
libsrc/mosaicing/im_remosaic.c
View File

@ -51,6 +51,7 @@
#include <vips/vips.h>
#include "transform.h"
#include "merge.h"
#include "global_balance.h"

1
libsrc/mosaicing/im_tbmerge.c
View File

@ -108,6 +108,7 @@
#include <vips/vips.h>
#include <vips/thread.h>
#include "transform.h"
#include "merge.h"
#ifdef WITH_DMALLOC

194
libsrc/mosaicing/interpolate.c
View File

@ -55,7 +55,6 @@
static VipsInterpolateClass *vips_interpolate_parent_class = NULL;
static VipsInterpolateClass *vips_interpolate_nearest_parent_class = NULL;
static VipsInterpolateClass *vips_interpolate_bilinear_parent_class = NULL;
static VipsInterpolateClass *vips_interpolate_bilinear_slow_parent_class = NULL;
#ifdef DEBUG
static void
@ -133,13 +132,13 @@ vips_interpolate_get_type( void )
* in_x, in_y in REGION in. Don't do this as a signal ffor speed.
*/
void
vips_interpolate( VipsInterpolate *interpolate, REGION *out, REGION *in,
int out_x, int out_y, double in_x, double in_y )
vips_interpolate( VipsInterpolate *interpolate,
PEL *out, REGION *in, double x, double y )
{
VipsInterpolateClass *class = VIPS_INTERPOLATE_GET_CLASS( interpolate );
g_assert( class->interpolate );
class->interpolate( interpolate, out, in, out_x, out_y, in_x, in_y );
class->interpolate( interpolate, out, in, x, y );
}
/* As above, but return the function pointer. Use this to cache method
@ -171,30 +170,22 @@ vips_interpolate_get_window_size( VipsInterpolate *interpolate )
static void
vips_interpolate_nearest_interpolate( VipsInterpolate *interpolate,
REGION *out, REGION *in,
int out_x, int out_y, double in_x, double in_y )
PEL *out, REGION *in, double x, double y )
{
/* Pel size and line size.
*/
const int ps = IM_IMAGE_SIZEOF_PEL( in->im );
int z;
PEL *q = (PEL *) IM_REGION_ADDR( out, out_x, out_y );
/* Subtract 0.5 to centre the nearest.
*/
const double cx = in_x - 0.5;
const double cy = in_y - 0.5;
/* Top left corner we interpolate from.
*/
const int xi = FLOOR( cx );
const int yi = FLOOR( cy );
const int xi = FLOOR( x );
const int yi = FLOOR( y );
const PEL *p = (PEL *) IM_REGION_ADDR( in, xi, yi );
for( z = 0; z < ps; z++ )
q[z] = p[z];
out[z] = p[z];
}
static void
@ -276,7 +267,7 @@ vips_interpolate_nearest_static( void )
/* Interpolate a section ... int8/16 types.
*/
#define BILINEAR_INT( TYPE ) { \
TYPE *tq = (TYPE *) q; \
TYPE *tq = (TYPE *) out; \
\
const int c1 = class->matrixi[xi][yi][0]; \
const int c2 = class->matrixi[xi][yi][1]; \
@ -296,7 +287,7 @@ vips_interpolate_nearest_static( void )
/* Interpolate a pel ... int32 and float types.
*/
#define BILINEAR_FLOAT( TYPE ) { \
TYPE *tq = (TYPE *) q; \
TYPE *tq = (TYPE *) out; \
\
const double c1 = class->matrixd[xi][yi][0]; \
const double c2 = class->matrixd[xi][yi][1]; \
@ -333,8 +324,7 @@ vips_interpolate_nearest_static( void )
static void
vips_interpolate_bilinear_interpolate( VipsInterpolate *interpolate,
REGION *out, REGION *in,
int out_x, int out_y, double in_x, double in_y )
PEL *out, REGION *in, double x, double y )
{
VipsInterpolateBilinearClass *class =
VIPS_INTERPOLATE_BILINEAR_GET_CLASS( interpolate );
@ -345,15 +335,10 @@ vips_interpolate_bilinear_interpolate( VipsInterpolate *interpolate,
const int ls = IM_REGION_LSKIP( in );
const int b = in->im->Bands;
/* Subtract 0.5 to centre the bilinear.
*/
const double cx = in_x - 0.5;
const double cy = in_y - 0.5;
/* Now go to scaled int.
*/
const double sx = cx * VIPS_TRANSFORM_SCALE;
const double sy = cy * VIPS_TRANSFORM_SCALE;
const double sx = x * VIPS_TRANSFORM_SCALE;
const double sy = y * VIPS_TRANSFORM_SCALE;
const int sxi = FLOOR( sx );
const int syi = FLOOR( sy );
@ -362,15 +347,14 @@ vips_interpolate_bilinear_interpolate( VipsInterpolate *interpolate,
*/
const int xi = sxi & (VIPS_TRANSFORM_SCALE - 1);
const int yi = syi & (VIPS_TRANSFORM_SCALE - 1);
const int in_x_int = sxi >> VIPS_TRANSFORM_SHIFT;
const int in_y_int = syi >> VIPS_TRANSFORM_SHIFT;
const int x_int = sxi >> VIPS_TRANSFORM_SHIFT;
const int y_int = syi >> VIPS_TRANSFORM_SHIFT;
const PEL *p1 = (PEL *) IM_REGION_ADDR( in, in_x_int, in_y_int );
const PEL *p1 = (PEL *) IM_REGION_ADDR( in, x_int, y_int );
const PEL *p2 = p1 + ps;
const PEL *p3 = p1 + ls;
const PEL *p4 = p3 + ps;
PEL *q = (PEL *) IM_REGION_ADDR( out, out_x, out_y );
int z;
SWITCH_INTERPOLATE( in->im->BandFmt,
@ -478,151 +462,3 @@ vips_interpolate_bilinear_static( void )
return( interpolate );
}
/* VipsInterpolateBilinearSlow class
*/
/* Slow mode is really just for testing ... it doesn't use the pre-calculated
* interpolation factors or the fixed-point arithmetic.
*/
/* in this class, name vars in the 2x2 grid as eg.
* p1 p2
* p3 p4
*/
/* Interpolate a pel ... don't use the pre-calcuated matricies.
*/
#define BILINEAR_SLOW( TYPE ) { \
TYPE *tq = (TYPE *) q; \
\
const TYPE *tp1 = (TYPE *) p1; \
const TYPE *tp2 = (TYPE *) p2; \
const TYPE *tp3 = (TYPE *) p3; \
const TYPE *tp4 = (TYPE *) p4; \
\
for( z = 0; z < b; z++ ) \
tq[z] = c1 * tp1[z] + c2 * tp2[z] + \
c3 * tp3[z] + c4 * tp4[z]; \
}
static void
vips_interpolate_bilinear_slow_interpolate( VipsInterpolate *interpolate,
REGION *out, REGION *in,
int out_x, int out_y, double in_x, double in_y )
{
/* Pel size and line size.
*/
const int ps = IM_IMAGE_SIZEOF_PEL( in->im );
const int ls = IM_REGION_LSKIP( in );
const int b = in->im->Bands;
/* Subtract 0.5 to centre the bilinear.
*/
const double cx = in_x - 0.5;
const double cy = in_y - 0.5;
/* Top left corner we interpolate from.
*/
const int xi = FLOOR( cx );
const int yi = FLOOR( cy );
/* Fractional part.
*/
const double X = cx - xi;
const double Y = cy - yi;
/* Residual.
*/
const double Xd = 1.0 - X;
const double Yd = 1.0 - Y;
/* Weights.
*/
const double c1 = Xd * Yd;
const double c2 = X * Yd;
const double c3 = Xd * Y;
const double c4 = X * Y;
const PEL *p1 = (PEL *) IM_REGION_ADDR( in, xi, yi );
const PEL *p2 = p1 + ps;
const PEL *p3 = p1 + ls;
const PEL *p4 = p3 + ps;
PEL *q = (PEL *) IM_REGION_ADDR( out, out_x, out_y );
int z;
SWITCH_INTERPOLATE( in->im->BandFmt,
BILINEAR_SLOW, BILINEAR_SLOW );
}
static void
vips_interpolate_bilinear_slow_class_init( VipsInterpolateBilinearSlowClass *class )
{
VipsInterpolateClass *interpolate_class =
(VipsInterpolateClass *) class;
vips_interpolate_bilinear_slow_parent_class =
g_type_class_peek_parent( class );
interpolate_class->interpolate =
vips_interpolate_bilinear_slow_interpolate;
interpolate_class->window_size = 2;
}
static void
vips_interpolate_bilinear_slow_init(
VipsInterpolateBilinearSlow *bilinear_slow )
{
#ifdef DEBUG
printf( "vips_interpolate_bilinear_slow_init: " );
vips_object_print( VIPS_OBJECT( bilinear_slow ) );
#endif /*DEBUG*/
}
GType
vips_interpolate_bilinear_slow_get_type( void )
{
static GType type = 0;
if( !type ) {
static const GTypeInfo info = {
sizeof( VipsInterpolateBilinearSlowClass ),
NULL, /* base_init */
NULL, /* base_finalize */
(GClassInitFunc)
vips_interpolate_bilinear_slow_class_init,
NULL, /* class_finalize */
NULL, /* class_data */
sizeof( VipsInterpolateBilinearSlow ),
32, /* n_preallocs */
(GInstanceInitFunc) vips_interpolate_bilinear_slow_init,
};
type = g_type_register_static( VIPS_TYPE_INTERPOLATE,
"VipsInterpolateBilinearSlow", &info, 0 );
}
return( type );
}
VipsInterpolate *
vips_interpolate_bilinear_slow_new( void )
{
return( VIPS_INTERPOLATE( g_object_new(
VIPS_TYPE_INTERPOLATE_BILINEAR_SLOW, NULL ) ) );
}
/* Convenience: return a static bilinear_slow you don't need to free.
*/
VipsInterpolate *
vips_interpolate_bilinear_slow_static( void )
{
static VipsInterpolate *interpolate = NULL;
if( !interpolate )
interpolate = vips_interpolate_bilinear_slow_new();
return( interpolate );
}

43
libsrc/mosaicing/merge.h
View File

@ -36,16 +36,6 @@
*/
#define BLEND_SCALE (4096)
/* How many bits of precision we keep for transformations in im_affine().
*/
#define TRANSFORM_SHIFT (5)
#define TRANSFORM_SCALE (1<<TRANSFORM_SHIFT)
/* How many bits of precision we keep for interpolation.
*/
#define INTERPOL_SHIFT (13)
#define INTERPOL_SCALE (1<<INTERPOL_SHIFT)
/* Keep state for each call in one of these.
*/
typedef struct _Overlapping {
@ -93,22 +83,6 @@ typedef struct _MergeInfo {
float *merge;
} MergeInfo;
/* Params for similarity transformation.
*/
typedef struct {
/* Geometry.
*/
Rect iarea; /* Area in input we can use */
Rect oarea; /* Area in output we generate */
/* The transform.
*/
double a, b, c, d;
double dx, dy;
double ia, ib, ic, id; /* Inverse of matrix abcd */
} Transformation;
/* Functions shared between lr and tb.
*/
extern double *im__coef1;
@ -139,20 +113,3 @@ int im__tbmerge1( IMAGE *ref, IMAGE *sec, IMAGE *out,
*/
int im__affine( IMAGE *in, IMAGE *out, Transformation *trn );
int
im_affinei( IMAGE *in, IMAGE *out, VipsInterpolate *interpolate,
double a, double b, double c, double d,
double dx, double dy,
int ox, int oy, int ow, int oh );
void im__transform_init( Transformation *trn );
int im__transform_calc_inverse( Transformation *trn );
int im__transform_isidentity( Transformation *trn );
void im__transform_forward( Transformation *trn,
double x, double y, double *ox, double *oy );
void im__transform_inverse( Transformation *trn,
double x, double y, double *ox, double *oy );
void im__transform_set_area( Transformation * );
int im__transform_add( Transformation *in1, Transformation *in2,
Transformation *out );
void im__transform_print( Transformation *trn );

9
libsrc/mosaicing/mosaic1.c
View File

@ -52,6 +52,7 @@
#include <vips/buf.h>
#include "mosaic.h"
#include "transform.h"
#include "merge.h"
#ifdef WITH_DMALLOC
@ -281,8 +282,8 @@ rotjoin_search( IMAGE *ref, IMAGE *sec, IMAGE *out, joinfn jfn,
/* Map points on sec to rotated image.
*/
im__transform_forward( &trn, xs1, ys1, &xs3, &ys3 );
im__transform_forward( &trn, xs2, ys2, &xs4, &ys4 );
im__transform_forward_point( &trn, xs1, ys1, &xs3, &ys3 );
im__transform_forward_point( &trn, xs2, ys2, &xs4, &ys4 );
/* Refine tie-points on rotated image. Remember the clip
* im__transform_set_area() has set, and move the sec tie-points
@ -315,8 +316,8 @@ rotjoin_search( IMAGE *ref, IMAGE *sec, IMAGE *out, joinfn jfn,
/* ... and now back to input space again.
*/
im__transform_inverse( &trn, xs5, ys5, &xs7, &ys7 );
im__transform_inverse( &trn, xs6, ys6, &xs8, &ys8 );
im__transform_invert_point( &trn, xs5, ys5, &xs7, &ys7 );
im__transform_invert_point( &trn, xs6, ys6, &xs8, &ys8 );
/* Recalc the transform using the refined points.
*/

13
libsrc/mosaicing/mosaicing_dispatch.c
View File

@ -39,6 +39,7 @@
#include <vips/vips.h>
#include <vips/internal.h>
#include "transform.h"
#include "merge.h"
#ifdef WITH_DMALLOC
@ -565,7 +566,7 @@ affinei_vec( im_object *argv )
break;
case 3:
interpolate = vips_interpolate_bilinear_slow_new();
interpolate = vips_interpolate_bicubic_new();
break;
case 4:
@ -573,16 +574,6 @@ affinei_vec( im_object *argv )
break;
case 5:
interpolate = vips_interpolate_yafrsmooth_new();
vips_interpolate_yafrsmooth_set_sharpening(
VIPS_INTERPOLATE_YAFRSMOOTH( interpolate ), 2.0 );
break;
case 6:
interpolate = vips_interpolate_yafr_test_new();
break;
case 7:
interpolate = vips_interpolate_yafrnohalo_new();
break;

37
libsrc/mosaicing/similarity.c
View File

@ -65,14 +65,13 @@
#include <vips/vips.h>
#include "transform.h"
#include "merge.h"
#ifdef WITH_DMALLOC
#include <dmalloc.h>
#endif /*WITH_DMALLOC*/
/* Call point from VIPS.
*/
int
im_similarity_area( IMAGE *in, IMAGE *out,
double a, double b, double dx, double dy,
@ -98,40 +97,6 @@ im_similarity_area( IMAGE *in, IMAGE *out,
return( im__affine( in, out, &trn ) );
}
/* Set output area of trn so that it just holds all of our input pels.
*/
void
im__transform_set_area( Transformation *trn )
{
double xA, xB, xC, xD;
double yA, yB, yC, yD;
int xmin, xmax, ymin, ymax;
im__transform_forward( trn,
trn->iarea.left, trn->iarea.top,
&xA, &yA );
im__transform_forward( trn,
IM_RECT_RIGHT( &trn->iarea ) - 1, trn->iarea.top,
&xB, &yB );
im__transform_forward( trn,
trn->iarea.left, IM_RECT_BOTTOM( &trn->iarea ) - 1,
&xC, &yC );
im__transform_forward( trn,
IM_RECT_RIGHT( &trn->iarea ) - 1,
IM_RECT_BOTTOM( &trn->iarea ) - 1,
&xD, &yD );
xmin = IM_MIN( xA, IM_MIN( xB, IM_MIN( xC, xD ) ) );
ymin = IM_MIN( yA, IM_MIN( yB, IM_MIN( yC, yD ) ) );
xmax = IM_MAX( xA, IM_MAX( xB, IM_MAX( xC, xD ) ) );
ymax = IM_MAX( yA, IM_MAX( yB, IM_MAX( yC, yD ) ) );
trn->oarea.left = xmin;
trn->oarea.top = ymin;
trn->oarea.width = xmax - xmin + 1;
trn->oarea.height = ymax - ymin + 1;
}
/* Output the rect holding all our input PELs.
*/
int

135
libsrc/mosaicing/templates.h Normal file
View File

@ -0,0 +1,135 @@
/* various interpolation templates
*/
/*
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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
/*
These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
*/
/* "fast" floor() ... on my laptop, anyway.
*/
#define FLOOR( V ) ((V) >= 0 ? (int)(V) : (int)((V) - 1))
#ifndef restrict
#ifdef __restrict
#define restrict __restrict
#else
#ifdef __restrict__
#define restrict __restrict__
#else
#define restrict
#endif
#endif
#endif
/* Fixed-point integer bicubic, used for 8 and 16-bit types.
*/
template <typename T> static int inline
bicubic_int(
const T uno_one, const T uno_two, const T uno_thr, const T uno_fou,
const T dos_one, const T dos_two, const T dos_thr, const T dos_fou,
const T tre_one, const T tre_two, const T tre_thr, const T tre_fou,
const T qua_one, const T qua_two, const T qua_thr, const T qua_fou,
const int *cx, const int *cy )
{
const int r0 =
(cx[0] * uno_one +
cx[1] * uno_two +
cx[2] * uno_thr +
cx[3] * uno_fou) >> VIPS_INTERPOLATE_SHIFT;
const int r1 =
(cx[0] * dos_one +
cx[1] * dos_two +
cx[2] * dos_thr +
cx[3] * dos_fou) >> VIPS_INTERPOLATE_SHIFT;
const int r2 =
(cx[0] * tre_one +
cx[1] * tre_two +
cx[2] * tre_thr +
cx[3] * tre_fou) >> VIPS_INTERPOLATE_SHIFT;
const int r3 =
(cx[0] * qua_one +
cx[1] * qua_two +
cx[2] * qua_thr +
cx[3] * qua_fou) >> VIPS_INTERPOLATE_SHIFT;
return( (cy[0] * r0 +
cy[1] * r1 +
cy[2] * r2 +
cy[3] * r3) >> VIPS_INTERPOLATE_SHIFT );
}
/* Floating-point bicubic, used for int/float/double types.
*/
template <typename T> static T inline
bicubic_float(
const T uno_one, const T uno_two, const T uno_thr, const T uno_fou,
const T dos_one, const T dos_two, const T dos_thr, const T dos_fou,
const T tre_one, const T tre_two, const T tre_thr, const T tre_fou,
const T qua_one, const T qua_two, const T qua_thr, const T qua_fou,
const double *cx, const double *cy )
{
return(
cy[0] * (cx[0] * uno_one +
cx[1] * uno_two +
cx[2] * uno_thr +
cx[3] * uno_fou) +
cy[1] * (cx[0] * dos_one +
cx[1] * dos_two +
cx[2] * dos_thr +
cx[3] * dos_fou) +
cy[2] * (cx[0] * tre_one +
cx[1] * tre_two +
cx[2] * tre_thr +
cx[3] * tre_fou) +
cy[3] * (cx[0] * qua_one +
cx[1] * qua_two +
cx[2] * qua_thr +
cx[3] * qua_fou) );
}
/* Given an offset in [0,1] (we can have x == 1 when building tables),
* calculate c0, c1, c2, c3, the catmull-rom coefficients. This is called
* from the interpolator, as well as from the table builder.
*/
static void inline
calculate_coefficients_catmull( const double x, double c[4] )
{
const double dx = 1.f - x;
const double x2 = dx * x;
const double mx2 = -.5f * x2;
g_assert( x >= 0 && x <= 1 );
c[0] = mx2 * dx;
c[1] = x2 * (-1.5f * x + 1.f) + dx;
c[2] = 1.f - (mx2 + c[1]);
c[3] = mx2 * x;
}

242
libsrc/mosaicing/transform.c Normal file
View File

@ -0,0 +1,242 @@
/* affine transforms
*/
/*
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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
/*
These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
*/
/*
*/
#define DEBUG
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif /*HAVE_CONFIG_H*/
#include <vips/intl.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <limits.h>
#include <vips/vips.h>
#include <vips/internal.h>
#include "transform.h"
#ifdef WITH_DMALLOC
#include <dmalloc.h>
#endif /*WITH_DMALLOC*/
/* Calculate the inverse transformation.
*/
int
im__transform_calc_inverse( Transformation *trn )
{
DOUBLEMASK *msk, *msk2;
if( !(msk = im_create_dmaskv( "boink", 2, 2,
trn->a, trn->b, trn->c, trn->d )) )
return( -1 );
if( !(msk2 = im_matinv( msk, "boink2" )) ) {
(void) im_free_dmask( msk );
return( -1 );
}
trn->ia = msk2->coeff[0];
trn->ib = msk2->coeff[1];
trn->ic = msk2->coeff[2];
trn->id = msk2->coeff[3];
(void) im_free_dmask( msk );
(void) im_free_dmask( msk2 );
return( 0 );
}
/* Init a Transform.
*/
void
im__transform_init( Transformation *trn )
{
trn->oarea.left = 0;
trn->oarea.top = 0;
trn->oarea.width = -1;
trn->oarea.height = -1;
trn->iarea.left = 0;
trn->iarea.top = 0;
trn->iarea.width = -1;
trn->iarea.height = -1;
trn->a = 1.0; /* Identity transform */
trn->b = 0.0;
trn->c = 0.0;
trn->d = 1.0;
trn->dx = 0.0;
trn->dy = 0.0;
(void) im__transform_calc_inverse( trn );
}
/* Test for transform is identity function.
*/
int
im__transform_isidentity( const Transformation *trn )
{
if( trn->a == 1.0 && trn->b == 0.0 && trn->c == 0.0 &&
trn->d == 1.0 && trn->dx == 0.0 && trn->dy == 0.0 )
return( 1 );
else
return( 0 );
}
/* Combine two transformations. out can be one of the ins.
*/
int
im__transform_add( const Transformation *in1, const Transformation *in2,
Transformation *out )
{
out->a = in1->a * in2->a + in1->c * in2->b;
out->b = in1->b * in2->a + in1->d * in2->b;
out->c = in1->a * in2->c + in1->c * in2->d;
out->d = in1->b * in2->c + in1->d * in2->d;
out->dx = in1->dx * in2->a + in1->dy * in2->b + in2->dx;
out->dy = in1->dx * in2->c + in1->dy * in2->d + in2->dy;
if( im__transform_calc_inverse( out ) )
return( -1 );
return( 0 );
}
void
im__transform_print( const Transformation *trn )
{
printf( "im__transform_print:\n" );
printf( " iarea: left=%d, top=%d, width=%d, height=%d\n",
trn->iarea.left,
trn->iarea.top,
trn->iarea.width,
trn->iarea.height );
printf( " oarea: left=%d, top=%d, width=%d, height=%d\n",
trn->oarea.left,
trn->oarea.top,
trn->oarea.width,
trn->oarea.height );
printf( " mat: a=%g, b=%g, c=%g, d=%g\n",
trn->a, trn->b, trn->c, trn->d );
printf( " off: dx=%g, dy=%g\n",
trn->dx, trn->dy );
}
/* Map a pixel coordinate through the transform.
*/
void
im__transform_forward_point( const Transformation *trn,
const double x, const double y, /* In input space */
double *ox, double *oy ) /* In output space */
{
*ox = trn->a * x + trn->b * y + trn->dx;
*oy = trn->c * x + trn->d * y + trn->dy;
}
/* Map a pixel coordinate through the inverse transform.
*/
void
im__transform_invert_point( const Transformation *trn,
const double x, const double y, /* In output space */
double *ox, double *oy ) /* In input space */
{
double mx = x - trn->dx;
double my = y - trn->dy;
*ox = trn->ia * mx + trn->ib * my;
*oy = trn->ic * mx + trn->id * my;
}
typedef void (*transform_fn)( const Transformation *,
const double, const double, double*, double* );
/* Transform a rect using a point transformer.
*/
static void
transform_rect( const Transformation *trn, transform_fn transform,
const Rect *in, /* In input space */
Rect *out ) /* In output space */
{
double x1, y1; /* Map corners */
double x2, y2;
double x3, y3;
double x4, y4;
double left, right, top, bottom;
/* Map input Rect.
*/
transform( trn, in->left, in->top, &x1, &y1 );
transform( trn, in->left, IM_RECT_BOTTOM( in ), &x3, &y3 );
transform( trn, IM_RECT_RIGHT( in ), in->top, &x2, &y2 );
transform( trn, IM_RECT_RIGHT( in ), IM_RECT_BOTTOM( in ), &x4, &y4 );
/* Find bounding box for these four corners.
*/
left = IM_MIN( x1, IM_MIN( x2, IM_MIN( x3, x4 ) ) );
right = IM_MAX( x1, IM_MAX( x2, IM_MAX( x3, x4 ) ) );
top = IM_MIN( y1, IM_MIN( y2, IM_MIN( y3, y4 ) ) );
bottom = IM_MAX( y1, IM_MAX( y2, IM_MAX( y3, y4 ) ) );
out->left = floor( left );
out->top = floor( top );
out->width = ceil( right ) - out->left;
out->height = ceil( bottom ) - out->top;
}
/* Given an area in the input image, calculate the bounding box for those
* pixels in the output image.
*/
void
im__transform_forward_rect( const Transformation *trn,
const Rect *in, /* In input space */
Rect *out ) /* In output space */
{
transform_rect( trn, im__transform_forward_point, in, out );
}
/* Given an area in the output image, calculate the bounding box for the
* corresponding pixels in the input image.
*/
void
im__transform_invert_rect( const Transformation *trn,
const Rect *in, /* In output space */
Rect *out ) /* In input space */
{
transform_rect( trn, im__transform_invert_point, in, out );
}
/* Set output area of trn so that it just holds all of our input pels.
*/
void
im__transform_set_area( Transformation *trn )
{
im__transform_forward_rect( trn, &trn->iarea, &trn->oarea );
}

67
libsrc/mosaicing/transform.h Normal file
View File

@ -0,0 +1,67 @@
/* Affine transforms.
*/
/*
Copyright (C) 1991-2003 The National Gallery
This program 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
/*
These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
*/
/* Params for an affine transformation.
*/
typedef struct {
/* Area of input we can use. This can be smaller than the real input
* image: we expand the input to add extra pixels for interpolation.
*/
Rect iarea;
/* The area of the output we've been asked to generate. left/top can
* be negative.
*/
Rect oarea;
/* The transform.
*/
double a, b, c, d;
double dx, dy;
double ia, ib, ic, id; /* Inverse of matrix abcd */
} Transformation;
void im__transform_init( Transformation *trn );
int im__transform_calc_inverse( Transformation *trn );
int im__transform_isidentity( const Transformation *trn );
int im__transform_add( const Transformation *in1, const Transformation *in2,
Transformation *out );
void im__transform_print( const Transformation *trn );
void im__transform_forward_point( const Transformation *trn,
const double x, const double y, double *ox, double *oy );
void im__transform_invert_point( const Transformation *trn,
const double x, const double y, double *ox, double *oy );
void im__transform_forward_rect( const Transformation *trn,
const Rect *in, Rect *out );
void im__transform_invert_rect( const Transformation *trn,
const Rect *in, Rect *out );
void im__transform_set_area( Transformation * );

6
libsrc/mosaicing/yafrnohalo.c
View File

@ -686,8 +686,7 @@ bilinear_yafrnohalo (float* restrict out, const float* restrict in,
static void
vips_interpolate_yafrnohalo_interpolate( VipsInterpolate *interpolate,
REGION *out, REGION *in,
int out_x, int out_y, double x, double y )
PEL *out, REGION *in, double x, double y )
{
VipsInterpolateYafrnohalo *yafrnohalo =
VIPS_INTERPOLATE_YAFRNOHALO( interpolate );
@ -977,11 +976,10 @@ vips_interpolate_yafrnohalo_interpolate( VipsInterpolate *interpolate,
/* Where we write the result.
*/
PEL *q = (PEL *) IM_REGION_ADDR( out, out_x, out_y );
int z;
for( z = 0; z < channels; z++ )
bilinear_yafrnohalo ((float *) q + z, (float *) p + z,
bilinear_yafrnohalo ((float *) out + z, (float *) p + z,
channels, pixels_per_buffer_row,
yafrnohalo->sharpening,
c_horizo_left___up,

759
libsrc/mosaicing/yafrsmooth.cpp Normal file
View File

@ -0,0 +1,759 @@
/* yafrsmooth interpolator
*/
/*
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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
/*
These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
*/
/*
* 2008 (c) Nicolas Robidoux (developer of Yet Another Fast
* Resampler).
*
* Acknowledgement: N. Robidoux's research on YAFRSMOOTH funded in part by
* an NSERC (National Science and Engineering Research Council of
* Canada) Discovery Grant.
*/
/* Hacked for vips by J. Cupitt, 12/11/08.
*
* Bicubic component replaced with the one from bicubbic.cpp.
*/
/*
* YAFRSMOOTH = Yet Another Fast Resampler
*
* Yet Another Fast Resampler is a nonlinear resampler which consists
* of a linear scheme (in this version, Catmull-Rom) plus a nonlinear
* sharpening correction the purpose of which is the straightening of
* diagonal interfaces between flat colour areas.
*
* Key properties:
*
* YAFRSMOOTH (smooth) is interpolatory:
*
* If asked for the value at the center of an input pixel, it will
* return the corresponding value, unchanged.
*
* YAFRSMOOTH (smooth) preserves local averages:
*
* The average of the reconstructed intensity surface over any region
* is the same as the average of the piecewise constant surface with
* values over pixel areas equal to the input pixel values (the
* "nearest neighbour" surface), except for a small amount of blur at
* the boundary of the region. More precicely: YAFRSMOOTH (smooth) is a box
* filtered exact area method.
*
* Main weaknesses of YAFRSMOOTH (smooth):
*
* Weakness 1: YAFRSMOOTH (smooth) improves on Catmull-Rom only for images
* with at least a little bit of smoothness.
*
* Weakness 2: Catmull-Rom introduces a lot of haloing. YAFRSMOOTH (smooth)
* is based on Catmull-Rom, and consequently it too introduces a lot
* of haloing.
*
* More details regarding Weakness 1:
*
* If a portion of the image is such that every pixel has immediate
* neighbours in the horizontal and vertical directions which have
* exactly the same pixel value, then YAFRSMOOTH (smooth) boils down to
* Catmull-Rom, and the computation of the correction is a waste.
* Extreme case: If all the pixels are either pure black or pure white
* in some region, as in some text images (more generally, if the
* region is "bichromatic"), then the YAFRSMOOTH (smooth) correction is 0 in
* the interior of the bichromatic region.
*/
/*
#define DEBUG
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif /*HAVE_CONFIG_H*/
#include <vips/intl.h>
#include <stdio.h>
#include <stdlib.h>
#include <vips/vips.h>
#include <vips/internal.h>
#include "templates.h"
#ifdef WITH_DMALLOC
#include <dmalloc.h>
#endif /*WITH_DMALLOC*/
/* "fast" floor() ... on my laptop, anyway.
*/
#define FLOOR( V ) ((V) >= 0 ? (int)(V) : (int)((V) - 1))
#ifndef restrict
#ifdef __restrict
#define restrict __restrict
#else
#ifdef __restrict__
#define restrict __restrict__
#else
#define restrict
#endif
#endif
#endif
/* Scale sharpening by this to normalise.
*/
#define SMOOTH_SHARPENING_SCALE (0.453125f)
static VipsInterpolateClass *vips_interpolate_yafrsmooth_parent_class = NULL;
/* T is the type of pixels we are computing, D is a type large enough to hold
* (Ta - Tb) ** 2.
*/
/* The 16 values for this interpolation, four constants for this
* interpolation position.
*/
template <typename T, typename D> static float inline
yafrsmooth(
const T uno_one, const T uno_two, const T uno_thr, const T uno_fou,
const T dos_one, const T dos_two, const T dos_thr, const T dos_fou,
const T tre_one, const T tre_two, const T tre_thr, const T tre_fou,
const T qua_one, const T qua_two, const T qua_thr, const T qua_fou,
const double *c )
{
/*
* Computation of the YAFRSMOOTH correction:
*
* Basically, if two consecutive pixel value differences have the
* same sign, the smallest one (in absolute value) is taken to be
* the corresponding slope. If they don't have the same sign, the
* corresponding slope is set to 0.
*
* Four such pairs (vertical and horizontal) of slopes need to be
* computed, one pair for each of the pixels which potentially
* overlap the unit area centered at the interpolation point.
*/
/*
* Beginning of the computation of the "up" horizontal slopes:
*/
const D prem__up = dos_two - dos_one;
const D deux__up = dos_thr - dos_two;
const D troi__up = dos_fou - dos_thr;
/*
* "down" horizontal slopes:
*/
const D prem_dow = tre_two - tre_one;
const D deux_dow = tre_thr - tre_two;
const D troi_dow = tre_fou - tre_thr;
/*
* "left" vertical slopes:
*/
const D prem_left = dos_two - uno_two;
const D deux_left = tre_two - dos_two;
const D troi_left = qua_two - tre_two;
/*
* "right" vertical slopes:
*/
const D prem_rite = dos_thr - uno_thr;
const D deux_rite = tre_thr - dos_thr;
const D troi_rite = qua_thr - tre_thr;
/*
* Back to "up":
*/
const D prem__up_squared = prem__up * prem__up;
const D deux__up_squared = deux__up * deux__up;
const D troi__up_squared = troi__up * troi__up;
/*
* Back to "down":
*/
const D prem_dow_squared = prem_dow * prem_dow;
const D deux_dow_squared = deux_dow * deux_dow;
const D troi_dow_squared = troi_dow * troi_dow;
/*
* Back to "left":
*/
const D prem_left_squared = prem_left * prem_left;
const D deux_left_squared = deux_left * deux_left;
const D troi_left_squared = troi_left * troi_left;
/*
* Back to "right":
*/
const D prem_rite_squared = prem_rite * prem_rite;
const D deux_rite_squared = deux_rite * deux_rite;
const D troi_rite_squared = troi_rite * troi_rite;
/*
* "up":
*/
const D prem__up_times_deux__up = prem__up * deux__up;
const D deux__up_times_troi__up = deux__up * troi__up;
/*
* "down":
*/
const D prem_dow_times_deux_dow = prem_dow * deux_dow;
const D deux_dow_times_troi_dow = deux_dow * troi_dow;
/*
* "left":
*/
const D prem_left_times_deux_left = prem_left * deux_left;
const D deux_left_times_troi_left = deux_left * troi_left;
/*
* "right":
*/
const D prem_rite_times_deux_rite = prem_rite * deux_rite;
const D deux_rite_times_troi_rite = deux_rite * troi_rite;
/*
* Branching parts of the computation of the YAFRSMOOTH correction
* (could be unbranched using arithmetic branching and C99 math
* intrinsics, although the compiler may be smart enough to remove
* the branching on its own):
*/
/*
* "up":
*/
const D prem__up_vs_deux__up =
prem__up_squared < deux__up_squared ? prem__up : deux__up;
const D deux__up_vs_troi__up =
deux__up_squared < troi__up_squared ? deux__up : troi__up;
/*
* "down":
*/
const D prem_dow_vs_deux_dow =
prem_dow_squared < deux_dow_squared ? prem_dow : deux_dow;
const D deux_dow_vs_troi_dow =
deux_dow_squared < troi_dow_squared ? deux_dow : troi_dow;
/*
* "left":
*/
const D prem_left_vs_deux_left =
prem_left_squared < deux_left_squared ? prem_left : deux_left;
const D deux_left_vs_troi_left =
deux_left_squared < troi_left_squared ? deux_left : troi_left;
/*
* "right":
*/
const D prem_rite_vs_deux_rite =
prem_rite_squared < deux_rite_squared ? prem_rite : deux_rite;
const D deux_rite_vs_troi_rite =
deux_rite_squared < troi_rite_squared ? deux_rite : troi_rite;
/*
* Computation of the YAFRSMOOTH slopes.
*/
/*
* "up":
*/
const D mx_left__up =
prem__up_times_deux__up < 0.f ? 0.f : prem__up_vs_deux__up;
const D mx_rite__up =
deux__up_times_troi__up < 0.f ? 0.f : deux__up_vs_troi__up;
/*
* "down":
*/
const D mx_left_dow =
prem_dow_times_deux_dow < 0.f ? 0.f : prem_dow_vs_deux_dow;
const D mx_rite_dow =
deux_dow_times_troi_dow < 0.f ? 0.f : deux_dow_vs_troi_dow;
/*
* "left":
*/
const D my_left__up =
prem_left_times_deux_left < 0.f ? 0.f : prem_left_vs_deux_left;
const D my_left_dow =
deux_left_times_troi_left < 0.f ? 0.f : deux_left_vs_troi_left;
/*
* "right":
*/
const D my_rite__up =
prem_rite_times_deux_rite < 0.f ? 0.f : prem_rite_vs_deux_rite;
const D my_rite_dow =
deux_rite_times_troi_rite < 0.f ? 0.f : deux_rite_vs_troi_rite;
/*
* Assemble the unweighted YAFRSMOOTH correction:
*/
const float yafr =
c[0] * (mx_left__up - mx_rite__up) +
c[1] * (mx_left_dow - mx_rite_dow) +
c[2] * (my_left__up - my_left_dow) +
c[3] * (my_rite__up - my_rite_dow);
return( yafr );
}
/* Pointers to write to / read from, number of bands,
* how many bytes to add to move down a line.
*/
/* T is the type of pixels we are reading and writing, D is a type large
* enough to hold (T1 - T2) ** 2.
*/
/* Fixed-point version for 8/16 bit ints.
*/
template <typename T, typename D, int min_value, int max_value>
static void inline
yafrsmooth_int_tab( PEL *pout, const PEL *pin,
const int bands, const int lskip,
const double sharpening,
const int *cx, const int *cy, const double *cs )
{
T* restrict out = (T *) pout;
const T* restrict in = (T *) pin;
const int b1 = bands;
const int b2 = 2 * bands;
const int b3 = 3 * bands;
const int l1 = lskip / sizeof( T );
const int l2 = 2 * lskip / sizeof( T );
const int l3 = 3 * lskip / sizeof( T );
for( int z = 0; z < bands; z++ ) {
const T uno_one = in[0];
const T uno_two = in[b1];
const T uno_thr = in[b2];
const T uno_fou = in[b3];
const T dos_one = in[l1];
const T dos_two = in[b1 + l1];
const T dos_thr = in[b2 + l1];
const T dos_fou = in[b3 + l1];
const T tre_one = in[l2];
const T tre_two = in[b1 + l2];
const T tre_thr = in[b2 + l2];
const T tre_fou = in[b3 + l2];
const T qua_one = in[l3];
const T qua_two = in[b1 + l3];
const T qua_thr = in[b2 + l3];
const T qua_fou = in[b3 + l3];
const int bicubic = bicubic_int<T>(
uno_one, uno_two, uno_thr, uno_fou,
dos_one, dos_two, dos_thr, dos_fou,
tre_one, tre_two, tre_thr, tre_fou,
qua_one, qua_two, qua_thr, qua_fou,
cx, cy );
const float yafr = yafrsmooth<T, D>(
uno_one, uno_two, uno_thr, uno_fou,
dos_one, dos_two, dos_thr, dos_fou,
tre_one, tre_two, tre_thr, tre_fou,
qua_one, qua_two, qua_thr, qua_fou,
cs );
int result = bicubic +
sharpening * SMOOTH_SHARPENING_SCALE * yafr;
if( result < min_value )
result = min_value;
else if( result > max_value )
result = max_value;
*out = result;
in += 1;
out += 1;
}
}
/* Float version for int/float types.
*/
template <typename T, typename D> static void inline
yafrsmooth_float_tab( PEL *pout, const PEL *pin,
const int bands, const int lskip,
const double sharpening,
const double *cx, const double *cy, const double *cs )
{
T* restrict out = (T *) pout;
const T* restrict in = (T *) pin;
const int b1 = bands;
const int b2 = 2 * bands;
const int b3 = 3 * bands;
const int l1 = lskip / sizeof( T );
const int l2 = 2 * lskip / sizeof( T );
const int l3 = 3 * lskip / sizeof( T );
for( int z = 0; z < bands; z++ ) {
const T uno_one = in[0];
const T uno_two = in[b1];
const T uno_thr = in[b2];
const T uno_fou = in[b3];
const T dos_one = in[l1];
const T dos_two = in[b1 + l1];
const T dos_thr = in[b2 + l1];
const T dos_fou = in[b3 + l1];
const T tre_one = in[l2];
const T tre_two = in[b1 + l2];
const T tre_thr = in[b2 + l2];
const T tre_fou = in[b3 + l2];
const T qua_one = in[l3];
const T qua_two = in[b1 + l3];
const T qua_thr = in[b2 + l3];
const T qua_fou = in[b3 + l3];
const T bicubic = bicubic_float<T>(
uno_one, uno_two, uno_thr, uno_fou,
dos_one, dos_two, dos_thr, dos_fou,
tre_one, tre_two, tre_thr, tre_fou,
qua_one, qua_two, qua_thr, qua_fou,
cx, cy );
const float yafr = yafrsmooth<T, D>(
uno_one, uno_two, uno_thr, uno_fou,
dos_one, dos_two, dos_thr, dos_fou,
tre_one, tre_two, tre_thr, tre_fou,
qua_one, qua_two, qua_thr, qua_fou,
cs );
*out = bicubic + sharpening * SMOOTH_SHARPENING_SCALE * yafr;
in += 1;
out += 1;
}
}
/* Given an offset in [0,1], calculate c0, c1, c2, c3, the yafr-smooth pixel
* weights.
*/
static void inline
calculate_coefficients_smooth( const double x, const double y, double c[4] )
{
const double dx = 1.f - x;
const double dy = 1.f - y;
g_assert( x >= 0 && x < 1 );
g_assert( y >= 0 && y < 1 );
c[0] = dx * x * dy;
c[1] = dx * x * y;
c[2] = dy * y * dx;
c[3] = dy * y * x;
}
/* High-quality double-only version.
*/
static void inline
yafrsmooth_notab( PEL *pout, const PEL *pin,
const int bands, const int lskip,
const double sharpening,
double x, double y )
{
double * restrict out = (double *) pout;
const double * restrict in = (double *) pin;
const int b1 = bands;
const int b2 = 2 * bands;
const int b3 = 3 * bands;
const int l1 = lskip / sizeof( double );
const int l2 = 2 * lskip / sizeof( double );
const int l3 = 3 * lskip / sizeof( double );
double cx[4];
double cy[4];
calculate_coefficients_catmull( x, cx );
calculate_coefficients_catmull( y, cy );
double cs[4];
calculate_coefficients_smooth( x, y, cs );
for( int z = 0; z < bands; z++ ) {
const double uno_one = in[0];
const double uno_two = in[b1];
const double uno_thr = in[b2];
const double uno_fou = in[b3];
const double dos_one = in[l1];
const double dos_two = in[b1 + l1];
const double dos_thr = in[b2 + l1];
const double dos_fou = in[b3 + l1];
const double tre_one = in[l2];
const double tre_two = in[b1 + l2];
const double tre_thr = in[b2 + l2];
const double tre_fou = in[b3 + l2];
const double qua_one = in[l3];
const double qua_two = in[b1 + l3];
const double qua_thr = in[b2 + l3];
const double qua_fou = in[b3 + l3];
const double bicubic = bicubic_float<double>(
uno_one, uno_two, uno_thr, uno_fou,
dos_one, dos_two, dos_thr, dos_fou,
tre_one, tre_two, tre_thr, tre_fou,
qua_one, qua_two, qua_thr, qua_fou,
cx, cy );
const double yafr = yafrsmooth<double, double>(
uno_one, uno_two, uno_thr, uno_fou,
dos_one, dos_two, dos_thr, dos_fou,
tre_one, tre_two, tre_thr, tre_fou,
qua_one, qua_two, qua_thr, qua_fou,
cs );
*out = bicubic + sharpening * SMOOTH_SHARPENING_SCALE * yafr;
in += 1;
out += 1;
}
}
static void
vips_interpolate_yafrsmooth_interpolate( VipsInterpolate *interpolate,
PEL *out, REGION *in, double x, double y )
{
VipsInterpolateYafrsmoothClass *yafrsmooth_class =
VIPS_INTERPOLATE_YAFRSMOOTH_GET_CLASS( interpolate );
VipsInterpolateYafrsmooth *yafrsmooth =
VIPS_INTERPOLATE_YAFRSMOOTH( interpolate );
/* Scaled int.
*/
const double sx = x * VIPS_TRANSFORM_SCALE;
const double sy = y * VIPS_TRANSFORM_SCALE;
const int sxi = FLOOR( sx );
const int syi = FLOOR( sy );
/* Get index into interpolation table and unscaled integer
* position.
*/
const int tx = sxi & (VIPS_TRANSFORM_SCALE - 1);
const int ty = syi & (VIPS_TRANSFORM_SCALE - 1);
const int xi = sxi >> VIPS_TRANSFORM_SHIFT;
const int yi = syi >> VIPS_TRANSFORM_SHIFT;
/* Look up the tables we need.
*/
const int *cxi = yafrsmooth_class->matrixi[tx];
const int *cyi = yafrsmooth_class->matrixi[ty];
const double *cxf = yafrsmooth_class->matrixf[tx];
const double *cyf = yafrsmooth_class->matrixf[ty];
/* Position weights for yafrsmooth.
*/
double cs[4];
calculate_coefficients_smooth( x - xi, y - yi, cs );
/* Back and up one to get the top-left of the 4x4.
*/
const PEL *p = (PEL *) IM_REGION_ADDR( in, xi - 1, yi - 1 );
/* Pel size and line size.
*/
const int bands = in->im->Bands;
const int lskip = IM_REGION_LSKIP( in );
#ifdef DEBUG
printf( "vips_interpolate_yafrsmooth_interpolate: %g %g\n", x, y );
printf( "\tleft=%d, top=%d, width=%d, height=%d\n",
xi - 1, yi - 1, 4, 4 );
#endif /*DEBUG*/
switch( in->im->BandFmt ) {
case IM_BANDFMT_UCHAR:
yafrsmooth_int_tab<unsigned char, int, 0, UCHAR_MAX>(
out, p, bands, lskip,
yafrsmooth->sharpening,
cxi, cyi, cs );
break;
case IM_BANDFMT_CHAR:
yafrsmooth_int_tab<signed char, int, SCHAR_MIN, SCHAR_MAX>(
out, p, bands, lskip,
yafrsmooth->sharpening,
cxi, cyi, cs );
break;
case IM_BANDFMT_USHORT:
yafrsmooth_int_tab<unsigned short, int, 0, USHRT_MAX>(
out, p, bands, lskip,
yafrsmooth->sharpening,
cxi, cyi, cs );
break;
case IM_BANDFMT_SHORT:
yafrsmooth_int_tab<signed short, int, SHRT_MIN, SHRT_MAX>(
out, p, bands, lskip,
yafrsmooth->sharpening,
cxi, cyi, cs );
break;
case IM_BANDFMT_UINT:
yafrsmooth_float_tab<unsigned int, float>(
out, p, bands, lskip,
yafrsmooth->sharpening,
cxf, cyf, cs );
break;
case IM_BANDFMT_INT:
yafrsmooth_float_tab<signed int, float>(
out, p, bands, lskip,
yafrsmooth->sharpening,
cxf, cyf, cs );
break;
case IM_BANDFMT_FLOAT:
yafrsmooth_float_tab<float, float>(
out, p, bands, lskip,
yafrsmooth->sharpening,
cxf, cyf, cs );
break;
case IM_BANDFMT_DOUBLE:
yafrsmooth_notab(
out, p, bands, lskip,
yafrsmooth->sharpening,
x - xi, y - yi );
break;
case IM_BANDFMT_COMPLEX:
yafrsmooth_float_tab<float, float>(
out, p, bands * 2, lskip,
yafrsmooth->sharpening,
cxf, cyf, cs );
break;
case IM_BANDFMT_DPCOMPLEX:
yafrsmooth_notab(
out, p, bands * 2, lskip,
yafrsmooth->sharpening,
x - xi, y - yi );
break;
default:
break;
}
}
static void
vips_interpolate_yafrsmooth_class_init( VipsInterpolateYafrsmoothClass *iclass )
{
VipsInterpolateClass *interpolate_class =
VIPS_INTERPOLATE_CLASS( iclass );
vips_interpolate_yafrsmooth_parent_class =
VIPS_INTERPOLATE_CLASS( g_type_class_peek_parent( iclass ) );
interpolate_class->interpolate =
vips_interpolate_yafrsmooth_interpolate;
interpolate_class->window_size = 4;
/* Build the tables of pre-computed coefficients.
*/
for( int x = 0; x < VIPS_TRANSFORM_SCALE + 1; x++ ) {
calculate_coefficients_catmull(
(float) x / VIPS_TRANSFORM_SCALE,
iclass->matrixf[x] );
for( int i = 0; i < 4; i++ )
iclass->matrixi[x][i] =
iclass->matrixf[x][i] * VIPS_INTERPOLATE_SCALE;
}
}
static void
vips_interpolate_yafrsmooth_init( VipsInterpolateYafrsmooth *yafrsmooth )
{
#ifdef DEBUG
printf( "vips_interpolate_yafrsmooth_init: " );
vips_object_print( VIPS_OBJECT( yafrsmooth ) );
#endif /*DEBUG*/
yafrsmooth->sharpening = 1.0;
}
GType
vips_interpolate_yafrsmooth_get_type()
{
static GType type = 0;
if( !type ) {
static const GTypeInfo info = {
sizeof( VipsInterpolateYafrsmoothClass ),
NULL, /* base_init */
NULL, /* base_finalize */
(GClassInitFunc) vips_interpolate_yafrsmooth_class_init,
NULL, /* class_finalize */
NULL, /* class_data */
sizeof( VipsInterpolateYafrsmooth ),
32, /* n_preallocs */
(GInstanceInitFunc) vips_interpolate_yafrsmooth_init,
};
type = g_type_register_static( VIPS_TYPE_INTERPOLATE,
"VipsInterpolateYafrsmooth", &info,
(GTypeFlags) 0 );
}
return( type );
}
VipsInterpolate *
vips_interpolate_yafrsmooth_new( void )
{
return( VIPS_INTERPOLATE( g_object_new(
VIPS_TYPE_INTERPOLATE_YAFRSMOOTH, NULL ) ) );
}
void
vips_interpolate_yafrsmooth_set_sharpening(
VipsInterpolateYafrsmooth *yafrsmooth,
double sharpening )
{
yafrsmooth->sharpening = sharpening;
}
/* Convenience: return a static yafrsmooth you don't need to free.
*/
VipsInterpolate *
vips_interpolate_yafrsmooth_static( void )
{
static VipsInterpolate *interpolate = NULL;
if( !interpolate )
interpolate = vips_interpolate_yafrsmooth_new();
return( interpolate );
}

15
libsrc/mosaicing/yafrsmooth.c → libsrc/mosaicing/yafrsmooth.old
View File

@ -461,8 +461,7 @@ catrom_yafrsmooth (float* restrict out, const float* restrict in,
static void
vips_interpolate_yafrsmooth_interpolate( VipsInterpolate *interpolate,
REGION *out, REGION *in,
int out_x, int out_y, double x, double y )
PEL *out, REGION *in, double x, double y )
{
VipsInterpolateYafrsmooth *yafrsmooth =
VIPS_INTERPOLATE_YAFRSMOOTH( interpolate );
@ -580,11 +579,10 @@ vips_interpolate_yafrsmooth_interpolate( VipsInterpolate *interpolate,
/* Where we write the result.
*/
PEL *q = (PEL *) IM_REGION_ADDR( out, out_x, out_y );
int z;
for( z = 0; z < channels; z++ )
catrom_yafrsmooth ((float *) q + z, (float *) p + z,
catrom_yafrsmooth ((float *) out + z, (float *) p + z,
channels, pixels_per_buffer_row,
yafrsmooth->sharpening,
cardinal_one,
@ -602,13 +600,13 @@ vips_interpolate_yafrsmooth_interpolate( VipsInterpolate *interpolate,
}
static void
vips_interpolate_yafrsmooth_class_init( VipsInterpolateYafrsmoothClass *class )
vips_interpolate_yafrsmooth_class_init( VipsInterpolateYafrsmoothClass *iclass )
{
VipsInterpolateClass *interpolate_class =
VIPS_INTERPOLATE_CLASS( class );
VIPS_INTERPOLATE_CLASS( iclass );
vips_interpolate_yafrsmooth_parent_class =
g_type_class_peek_parent( class );
VIPS_INTERPOLATE_CLASS( g_type_class_peek_parent( iclass ) );
interpolate_class->interpolate =
vips_interpolate_yafrsmooth_interpolate;
@ -645,7 +643,8 @@ vips_interpolate_yafrsmooth_get_type( void )
};
type = g_type_register_static( VIPS_TYPE_INTERPOLATE,
"VipsInterpolateYafrsmooth", &info, 0 );
"VipsInterpolateYafrsmooth", &info,
(GTypeFlags) 0 );
}
return( type );

9
libsrc/mosaicing/yafrtest.cpp
View File

@ -435,8 +435,7 @@ catrom_yafr_test(
static void
vips_interpolate_yafr_test_interpolate( VipsInterpolate *interpolate,
REGION *out, REGION *in,
int out_x, int out_y, double x, double y )
PEL *out, REGION *in, double x, double y )
{
VipsInterpolateYafrTest *yafr_test =
VIPS_INTERPOLATE_YAFR_TEST( interpolate );
@ -552,14 +551,10 @@ vips_interpolate_yafr_test_interpolate( VipsInterpolate *interpolate,
const int lskip = IM_REGION_LSKIP( in );
const int esize = IM_IMAGE_SIZEOF_ELEMENT( in->im );
/* Where we write the result.
*/
PEL *q = (PEL *) IM_REGION_ADDR( out, out_x, out_y );
/* Put this in a macro to save some typing.
*/
#define CALL(T, D) \
catrom_yafr_test<T, D>(q + z * esize, p + z * esize, \
catrom_yafr_test<T, D>(out + z * esize, p + z * esize, \
channels, lskip, \
yafr_test->sharpening, \
cardinal_one, \

0
po/LINGUAS Normal file
View File

BIN
po/en_GB.gmo
View File

Binary file not shown.

588
po/en_GB.po
View File

File diff suppressed because it is too large Load Diff