forked from libvips/libvips
/
vsqbs.cpp
415 lines (366 loc) · 14 KB
/
vsqbs.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
/* vertex-split subdivision followed by quadratic b-spline smoothing
*
* C. Racette 23-28/05/2010 based on code by N. Robidoux and J. Cupitt
*
* N. Robidoux 29-30/05/2010
*/
/*
This file is part of VIPS.
VIPS is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as
published by the Free Software Foundation; either version 2 of the
License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
02110-1301 USA
*/
/*
These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
*/
/*
* 2010 (c) Chantal Racette, Nicolas Robidoux, John Cupitt.
*
* Nicolas Robidoux thanks Adam Turcotte, Geert Jordaens, Ralf Meyer,
* Øyvind Kolås, Minglun Gong and Eric Daoust for useful comments and
* code.
*
* Chantal Racette's image resampling research and programming funded
* in part by a NSERC Discovery Grant awarded to Julien Dompierre
* (20-61098).
*/
/*
* Vertex-Split Quadratic B-Splines (VSQBS) is a brand new method
* which consists of vertex-split subdivision, a subdivision method
* with the (as yet unknown?) property that data which is (locally)
* constant on diagonals is subdivided into data which is (locally)
* constant on diagonals, followed by quadratic B-Spline smoothing.
* Because both methods are linear, their combination can be
* implemented as if there is no subdivision.
*
* At high enlargement ratios, VSQBS is very effective at "masking"
* that the original has pixels uniformly distributed on a grid. In
* particular, VSQBS produces resamples with only very mild
* staircasing. Like cubic B-Spline smoothing, however, VSQBS is not
* an interpolatory method. For example, using VSQBS to perform the
* identity geometric transformation (enlargement by a scaling factor
* equal to 1) on an image does not return the original: VSQBS
* effectively smooths out the image with the convolution mask
*
* 1/8
* 1/8 1/2 1/8
* 1/8
*
* which is a fairly moderate blur (although the checkerboard mode is
* in its nullspace).
*
* By blending VSQBS with an interpolatory method (bilinear, say) in a
* transformation adaptive environment (current GEGL, for example), it
* is quite easy to restore that resampling for identity geometric
* transformation is equivalent to the identity, and rotations are not
* affected by the above, implicit, blur. Contact N. Robidoux for
* details.
*
* An article on VSQBS is forthcoming.
*/
#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"
#define VIPS_TYPE_INTERPOLATE_VSQBS \
(vips_interpolate_vsqbs_get_type())
#define VIPS_INTERPOLATE_VSQBS( obj ) \
(G_TYPE_CHECK_INSTANCE_CAST( (obj), \
VIPS_TYPE_INTERPOLATE_VSQBS, VipsInterpolateVsqbs ))
#define VIPS_INTERPOLATE_VSQBS_CLASS( klass ) \
(G_TYPE_CHECK_CLASS_CAST( (klass), \
VIPS_TYPE_INTERPOLATE_VSQBS, VipsInterpolateVsqbsClass))
#define VIPS_IS_INTERPOLATE_VSQBS( obj ) \
(G_TYPE_CHECK_INSTANCE_TYPE( (obj), VIPS_TYPE_INTERPOLATE_VSQBS ))
#define VIPS_IS_INTERPOLATE_VSQBS_CLASS( klass ) \
(G_TYPE_CHECK_CLASS_TYPE( (klass), VIPS_TYPE_INTERPOLATE_VSQBS ))
#define VIPS_INTERPOLATE_VSQBS_GET_CLASS( obj ) \
(G_TYPE_INSTANCE_GET_CLASS( (obj), \
VIPS_TYPE_INTERPOLATE_VSQBS, VipsInterpolateVsqbsClass ))
typedef struct _VipsInterpolateVsqbs {
VipsInterpolate parent_object;
} VipsInterpolateVsqbs;
typedef struct _VipsInterpolateVsqbsClass {
VipsInterpolateClass parent_class;
} VipsInterpolateVsqbsClass;
/*
* THE STENCIL OF INPUT VALUES:
*
* Pointer arithmetic is used to implicitly reflect the input stencil
* about dos_two---assumed closer to the sampling location than other
* pixels (ties are OK)---in such a way that after reflection the
* sampling point is to the bottom right of dos_two.
*
* The following code and picture assumes that the stencil reflexion
* has already been performed. (X is the sampling location.)
*
*
* (ix,iy-1) (ix+1,iy-1)
* = uno_two = uno_thr
*
*
*
* (ix-1,iy) (ix,iy) (ix+1,iy)
* = dos_one = dos_two = dos_thr
* X
*
*
* (ix-1,iy+1) (ix,iy+1) (ix+1,iy+1)
* = tre_one = tre_two = tre_thr
*
*
* The above input pixel values are the ones needed in order to
* IMPLICITLY make available the following values, needed by quadratic
* B-Splines, which is performed on (shifted) double density data:
*
*
* uno_one_1 = uno_two_1 = uno_thr_1 =
* (ix-1/4,iy-1/4) (ix+1/4,iy-1/4) (ix+3/4,iy-1/4)
*
*
*
* X or X
* dos_one_1 = dos_two_1 = dos_thr_1 =
* (ix-1/4,iy+1/4) (ix+1/4,iy+1/4) (ix+3/4,iy+1/4)
* or X or X
*
*
*
* tre_one_1 = tre_two_1 = tre_thr_1 =
* (ix-1/4,iy+3/4) (ix+1/4,iy+3/4) (ix+3/4,iy+3/4)
*
*
* In the coefficient computations, we fix things so that coordinates
* are relative to dos_two_1, and so that distances are rescaled so
* that double density pixel locations are at a distance of 1.
*/
/*
* Call vertex-split + quadratic B-splines with a careful type
* conversion as a parameter. (It would be nice to do this with
* templates somehow---for one thing this would allow code
* comments---but we can't figure a clean way to do it.)
*/
#define VSQBS_CONVERSION( conversion ) \
template <typename T> static void inline \
vsqbs_ ## conversion( void* restrict pout, \
const VipsPel* restrict pin, \
const int bands, \
const int lskip, \
const double x_0, \
const double y_0 ) \
{ \
T* restrict out = (T *) pout; \
\
const T* restrict in = (T *) pin; \
\
const int sign_of_x_0 = 2 * ( x_0 >= 0. ) - 1; \
const int sign_of_y_0 = 2 * ( y_0 >= 0. ) - 1; \
\
const int shift_forw_1_pix = sign_of_x_0 * bands; \
const int shift_forw_1_row = sign_of_y_0 * lskip; \
\
const int shift_back_1_pix = -shift_forw_1_pix; \
const int shift_back_1_row = -shift_forw_1_row; \
\
const int uno_two_shift = shift_back_1_row; \
const int uno_thr_shift = shift_forw_1_pix + shift_back_1_row; \
\
const int dos_one_shift = shift_back_1_pix; \
const int dos_two_shift = 0; \
const int dos_thr_shift = shift_forw_1_pix; \
\
const int tre_one_shift = shift_back_1_pix + shift_forw_1_row; \
const int tre_two_shift = shift_forw_1_row; \
const int tre_thr_shift = shift_forw_1_pix + shift_forw_1_row; \
\
\
const double twice_abs_x_0 = ( 2 * sign_of_x_0 ) * x_0; \
const double twice_abs_y_0 = ( 2 * sign_of_y_0 ) * y_0; \
const double x = twice_abs_x_0 + -0.5; \
const double y = twice_abs_y_0 + -0.5; \
const double cent = 0.75 - x * x; \
const double mid = 0.75 - y * y; \
const double left = -0.5 * ( x + cent ) + 0.5; \
const double top = -0.5 * ( y + mid ) + 0.5; \
const double left_p_cent = left + cent; \
const double top_p_mid = top + mid; \
const double cent_p_rite = 1.0 - left; \
const double mid_p_bot = 1.0 - top; \
const double rite = 1.0 - left_p_cent; \
const double bot = 1.0 - top_p_mid; \
\
const double four_c_uno_two = left_p_cent * top; \
const double four_c_dos_one = left * top_p_mid; \
const double four_c_dos_two = left_p_cent + top_p_mid; \
const double four_c_dos_thr = cent_p_rite * top_p_mid + rite; \
const double four_c_tre_two = mid_p_bot * left_p_cent + bot; \
const double four_c_tre_thr = mid_p_bot * rite + cent_p_rite * bot; \
const double four_c_uno_thr = top - four_c_uno_two; \
const double four_c_tre_one = left - four_c_dos_one; \
\
\
int band = bands; \
\
do \
{ \
const double double_result = \
( \
( \
( \
four_c_uno_two * in[uno_two_shift] \
+ \
four_c_dos_one * in[dos_one_shift] \
) \
+ \
( \
four_c_dos_two * in[dos_two_shift] \
+ \
four_c_dos_thr * in[dos_thr_shift] \
) \
) \
+ \
( \
( \
four_c_tre_two * in[tre_two_shift] \
+ \
four_c_tre_thr * in[tre_thr_shift] \
) \
+ \
( \
four_c_uno_thr * in[uno_thr_shift] \
+ \
four_c_tre_one * in[tre_one_shift] \
) \
) \
) * 0.25; \
\
const T result = to_ ## conversion<T>( double_result ); \
in++; \
*out++ = result; \
\
} while (--band); \
}
VSQBS_CONVERSION( fptypes )
VSQBS_CONVERSION( withsign )
VSQBS_CONVERSION( nosign )
#define CALL( T, conversion ) \
vsqbs_ ## conversion<T>( out, \
p, \
bands, \
lskip, \
relative_x, \
relative_y );
/*
* We need C linkage:
*/
extern "C" {
G_DEFINE_TYPE( VipsInterpolateVsqbs, vips_interpolate_vsqbs,
VIPS_TYPE_INTERPOLATE );
}
static void
vips_interpolate_vsqbs_interpolate( VipsInterpolate* restrict interpolate,
void* restrict out,
VipsRegion* restrict in,
double absolute_x,
double absolute_y )
{
/* absolute_x and absolute_y are always >= 1.0 (see double-check assert
* below), so we don't need floor().
*
* It's 1 not 0 since we ask for a window_offset of 1 at the bottom.
*/
const int ix = (int) (absolute_x + 0.5);
const int iy = (int) (absolute_y + 0.5);
/*
* Move the pointer to (the first band of) the top/left pixel of the
* 2x2 group of pixel centers which contains the sampling location
* in its convex hull:
*/
const VipsPel* restrict p = VIPS_REGION_ADDR( in, ix, iy );
const double relative_x = absolute_x - ix;
const double relative_y = absolute_y - iy;
/*
* VIPS versions of Nicolas's pixel addressing values.
*/
const int lskip = VIPS_REGION_LSKIP( in ) /
VIPS_IMAGE_SIZEOF_ELEMENT( in->im );
/*
* Double the bands for complex images to account for the real and
* imaginary parts being computed independently:
*/
const int actual_bands = in->im->Bands;
const int bands =
vips_band_format_iscomplex( in->im->BandFmt ) ?
2 * actual_bands : actual_bands;
g_assert( ix - 1 >= in->valid.left );
g_assert( iy - 1 >= in->valid.top );
g_assert( ix + 1 <= VIPS_RECT_RIGHT( &in->valid ) );
g_assert( iy + 1 <= VIPS_RECT_BOTTOM( &in->valid ) );
/* Confirm that absolute_x and absolute_y are >= 1, see above.
*/
g_assert( absolute_x >= 1.0 );
g_assert( absolute_y >= 1.0 );
switch( in->im->BandFmt ) {
case VIPS_FORMAT_UCHAR:
CALL( unsigned char, nosign );
break;
case VIPS_FORMAT_CHAR:
CALL( signed char, withsign );
break;
case VIPS_FORMAT_USHORT:
CALL( unsigned short, nosign );
break;
case VIPS_FORMAT_SHORT:
CALL( signed short, withsign );
break;
case VIPS_FORMAT_UINT:
CALL( unsigned int, nosign );
break;
case VIPS_FORMAT_INT:
CALL( signed int, withsign );
break;
/*
* Complex images are handled by doubling bands:
*/
case VIPS_FORMAT_FLOAT:
case VIPS_FORMAT_COMPLEX:
CALL( float, fptypes );
break;
case VIPS_FORMAT_DOUBLE:
case VIPS_FORMAT_DPCOMPLEX:
CALL( double, fptypes );
break;
default:
g_assert( 0 );
break;
}
}
static void
vips_interpolate_vsqbs_class_init( VipsInterpolateVsqbsClass *klass )
{
VipsObjectClass *object_class = VIPS_OBJECT_CLASS( klass );
VipsInterpolateClass *interpolate_class = VIPS_INTERPOLATE_CLASS( klass );
object_class->nickname = "vsqbs";
object_class->description = _( "B-Splines with antialiasing smoothing" );
interpolate_class->interpolate = vips_interpolate_vsqbs_interpolate;
interpolate_class->window_size = 4;
interpolate_class->window_offset = 1;
}
static void
vips_interpolate_vsqbs_init( VipsInterpolateVsqbs *vsqbs )
{
}