1 /* horizontal reduce by a float factor with a kernel
2 *
3 * 29/1/16
4 * - from shrinkh.c
5 * 10/3/16
6 * - add other kernels
7 * 15/8/16
8 * - rename xshrink as hshrink for consistency
9 * 9/9/16
10 * - add @centre option
11 * 6/6/20 kleisauke
12 * - deprecate @centre option, it's now always on
13 * - fix pixel shift
14 */
15
16 /*
17
18 This file is part of VIPS.
19
20 VIPS is free software; you can redistribute it and/or modify
21 it under the terms of the GNU Lesser General Public License as published by
22 the Free Software Foundation; either version 2 of the License, or
23 (at your option) any later version.
24
25 This program is distributed in the hope that it will be useful,
26 but WITHOUT ANY WARRANTY; without even the implied warranty of
27 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
28 GNU Lesser General Public License for more details.
29
30 You should have received a copy of the GNU Lesser General Public License
31 along with this program; if not, write to the Free Software
32 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
33 02110-1301 USA
34
35 */
36
37 /*
38
39 These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
40
41 */
42
43 /*
44 #define DEBUG
45 */
46
47 #ifdef HAVE_CONFIG_H
48 #include <config.h>
49 #endif /*HAVE_CONFIG_H*/
50 #include <vips/intl.h>
51
52 #include <stdio.h>
53 #include <stdlib.h>
54 #include <math.h>
55
56 #include <vips/vips.h>
57 #include <vips/debug.h>
58 #include <vips/internal.h>
59
60 #include "presample.h"
61 #include "templates.h"
62
63 typedef struct _VipsReduceh {
64 VipsResample parent_instance;
65
66 double hshrink; /* Reduce factor */
67
68 /* The thing we use to make the kernel.
69 */
70 VipsKernel kernel;
71
72 /* Number of points in kernel.
73 */
74 int n_point;
75
76 /* Horizontal displacement.
77 */
78 double hoffset;
79
80 /* Precalculated interpolation matrices. int (used for pel
81 * sizes up to short), and double (for all others). We go to
82 * scale + 1 so we can round-to-nearest safely.
83 */
84 int *matrixi[VIPS_TRANSFORM_SCALE + 1];
85 double *matrixf[VIPS_TRANSFORM_SCALE + 1];
86
87 /* Deprecated.
88 */
89 gboolean centre;
90
91 } VipsReduceh;
92
93 typedef VipsResampleClass VipsReducehClass;
94
95 /* We need C linkage for this.
96 */
97 extern "C" {
98 G_DEFINE_TYPE( VipsReduceh, vips_reduceh, VIPS_TYPE_RESAMPLE );
99 }
100
101 /* Get n points. @shrink is the shrink factor, so 2 for a 50% reduction.
102 */
103 int
vips_reduce_get_points(VipsKernel kernel,double shrink)104 vips_reduce_get_points( VipsKernel kernel, double shrink )
105 {
106 switch( kernel ) {
107 case VIPS_KERNEL_NEAREST:
108 return( 1 );
109
110 case VIPS_KERNEL_LINEAR:
111 return( 2 * rint( shrink ) + 1 );
112
113 case VIPS_KERNEL_CUBIC:
114 case VIPS_KERNEL_MITCHELL:
115 return( 2 * rint( 2 * shrink ) + 1 );
116
117 case VIPS_KERNEL_LANCZOS2:
118 /* Needs to be in sync with calculate_coefficients_lanczos().
119 */
120 return( 2 * rint( 2 * shrink ) + 1 );
121
122 case VIPS_KERNEL_LANCZOS3:
123 return( 2 * rint( 3 * shrink ) + 1 );
124
125 default:
126 g_assert_not_reached();
127 return( 0 );
128 }
129 }
130
131 /* Calculate a mask element.
132 */
133 void
vips_reduce_make_mask(double * c,VipsKernel kernel,double shrink,double x)134 vips_reduce_make_mask( double *c, VipsKernel kernel, double shrink, double x )
135 {
136 switch( kernel ) {
137 case VIPS_KERNEL_NEAREST:
138 c[0] = 1.0;
139 break;
140
141 case VIPS_KERNEL_LINEAR:
142 calculate_coefficients_triangle( c, shrink, x );
143 break;
144
145 case VIPS_KERNEL_CUBIC:
146 /* Catmull-Rom.
147 */
148 calculate_coefficients_cubic( c, shrink, x, 0.0, 0.5 );
149 break;
150
151 case VIPS_KERNEL_MITCHELL:
152 calculate_coefficients_cubic( c, shrink, x,
153 1.0 / 3.0, 1.0 / 3.0 );
154 break;
155
156 case VIPS_KERNEL_LANCZOS2:
157 calculate_coefficients_lanczos( c, 2, shrink, x );
158 break;
159
160 case VIPS_KERNEL_LANCZOS3:
161 calculate_coefficients_lanczos( c, 3, shrink, x );
162 break;
163
164 default:
165 g_assert_not_reached();
166 break;
167 }
168 }
169
170 template <typename T, int max_value>
171 static void inline
reduceh_unsigned_int_tab(VipsReduceh * reduceh,VipsPel * pout,const VipsPel * pin,const int bands,const int * restrict cx)172 reduceh_unsigned_int_tab( VipsReduceh *reduceh,
173 VipsPel *pout, const VipsPel *pin,
174 const int bands, const int * restrict cx )
175 {
176 T* restrict out = (T *) pout;
177 const T* restrict in = (T *) pin;
178 const int n = reduceh->n_point;
179
180 for( int z = 0; z < bands; z++ ) {
181 int sum;
182
183 sum = reduce_sum<T, int>( in + z, bands, cx, n );
184 sum = unsigned_fixed_round( sum );
185 sum = VIPS_CLIP( 0, sum, max_value );
186
187 out[z] = sum;
188 }
189 }
190
191 template <typename T, int min_value, int max_value>
192 static void inline
reduceh_signed_int_tab(VipsReduceh * reduceh,VipsPel * pout,const VipsPel * pin,const int bands,const int * restrict cx)193 reduceh_signed_int_tab( VipsReduceh *reduceh,
194 VipsPel *pout, const VipsPel *pin,
195 const int bands, const int * restrict cx )
196 {
197 T* restrict out = (T *) pout;
198 const T* restrict in = (T *) pin;
199 const int n = reduceh->n_point;
200
201 for( int z = 0; z < bands; z++ ) {
202 int sum;
203
204 sum = reduce_sum<T, int>( in + z, bands, cx, n );
205 sum = signed_fixed_round( sum );
206 sum = VIPS_CLIP( min_value, sum, max_value );
207
208 out[z] = sum;
209 }
210 }
211
212 /* Floating-point version.
213 */
214 template <typename T>
215 static void inline
reduceh_float_tab(VipsReduceh * reduceh,VipsPel * pout,const VipsPel * pin,const int bands,const double * cx)216 reduceh_float_tab( VipsReduceh *reduceh,
217 VipsPel *pout, const VipsPel *pin,
218 const int bands, const double *cx )
219 {
220 T* restrict out = (T *) pout;
221 const T* restrict in = (T *) pin;
222 const int n = reduceh->n_point;
223
224 for( int z = 0; z < bands; z++ )
225 out[z] = reduce_sum<T, double>( in + z, bands, cx, n );
226 }
227
228 /* 32-bit int output needs a double intermediate.
229 */
230
231 template <typename T, int max_value>
232 static void inline
reduceh_unsigned_int32_tab(VipsReduceh * reduceh,VipsPel * pout,const VipsPel * pin,const int bands,const double * restrict cx)233 reduceh_unsigned_int32_tab( VipsReduceh *reduceh,
234 VipsPel *pout, const VipsPel *pin,
235 const int bands, const double * restrict cx )
236 {
237 T* restrict out = (T *) pout;
238 const T* restrict in = (T *) pin;
239 const int n = reduceh->n_point;
240
241 for( int z = 0; z < bands; z++ ) {
242 double sum;
243
244 sum = reduce_sum<T, double>( in + z, bands, cx, n );
245 out[z] = VIPS_CLIP( 0, sum, max_value );
246 }
247 }
248
249 template <typename T, int min_value, int max_value>
250 static void inline
reduceh_signed_int32_tab(VipsReduceh * reduceh,VipsPel * pout,const VipsPel * pin,const int bands,const double * restrict cx)251 reduceh_signed_int32_tab( VipsReduceh *reduceh,
252 VipsPel *pout, const VipsPel *pin,
253 const int bands, const double * restrict cx )
254 {
255 T* restrict out = (T *) pout;
256 const T* restrict in = (T *) pin;
257 const int n = reduceh->n_point;
258
259 for( int z = 0; z < bands; z++ ) {
260 double sum;
261
262 sum = reduce_sum<T, double>( in + z, bands, cx, n );
263 sum = VIPS_CLIP( min_value, sum, max_value );
264 out[z] = sum;
265 }
266 }
267
268 /* Ultra-high-quality version for double images.
269 */
270 template <typename T>
271 static void inline
reduceh_notab(VipsReduceh * reduceh,VipsPel * pout,const VipsPel * pin,const int bands,double x)272 reduceh_notab( VipsReduceh *reduceh,
273 VipsPel *pout, const VipsPel *pin,
274 const int bands, double x )
275 {
276 T* restrict out = (T *) pout;
277 const T* restrict in = (T *) pin;
278 const int n = reduceh->n_point;
279
280 double cx[MAX_POINT];
281
282 vips_reduce_make_mask( cx, reduceh->kernel, reduceh->hshrink, x );
283
284 for( int z = 0; z < bands; z++ ) {
285 double sum;
286 sum = reduce_sum<T, double>( in + z, bands, cx, n );
287
288 out[z] = VIPS_ROUND_UINT( sum );
289 }
290 }
291
292 /* Tried a vector path (see reducev) but it was slower. The vectors for
293 * horizontal reduce are just too small to get a useful speedup.
294 */
295
296 static int
vips_reduceh_gen(VipsRegion * out_region,void * seq,void * a,void * b,gboolean * stop)297 vips_reduceh_gen( VipsRegion *out_region, void *seq,
298 void *a, void *b, gboolean *stop )
299 {
300 VipsImage *in = (VipsImage *) a;
301 VipsReduceh *reduceh = (VipsReduceh *) b;
302 const int ps = VIPS_IMAGE_SIZEOF_PEL( in );
303 VipsRegion *ir = (VipsRegion *) seq;
304 VipsRect *r = &out_region->valid;
305
306 /* Double bands for complex.
307 */
308 const int bands = in->Bands *
309 (vips_band_format_iscomplex( in->BandFmt ) ? 2 : 1);
310
311 VipsRect s;
312
313 #ifdef DEBUG
314 printf( "vips_reduceh_gen: generating %d x %d at %d x %d\n",
315 r->width, r->height, r->left, r->top );
316 #endif /*DEBUG*/
317
318 s.left = r->left * reduceh->hshrink - reduceh->hoffset;
319 s.top = r->top;
320 s.width = r->width * reduceh->hshrink + reduceh->n_point;
321 s.height = r->height;
322 if( vips_region_prepare( ir, &s ) )
323 return( -1 );
324
325 VIPS_GATE_START( "vips_reduceh_gen: work" );
326
327 for( int y = 0; y < r->height; y++ ) {
328 VipsPel *p0;
329 VipsPel *q;
330
331 double X;
332
333 q = VIPS_REGION_ADDR( out_region, r->left, r->top + y );
334
335 X = (r->left + 0.5) * reduceh->hshrink - 0.5 -
336 reduceh->hoffset;
337
338 /* We want p0 to be the start (ie. x == 0) of the input
339 * scanline we are reading from. We can then calculate the p we
340 * need for each pixel with a single mul and avoid calling ADDR
341 * for each pixel.
342 *
343 * We can't get p0 directly with ADDR since it could be outside
344 * valid, so get the leftmost pixel in valid and subtract a
345 * bit.
346 */
347 p0 = VIPS_REGION_ADDR( ir, ir->valid.left, r->top + y ) -
348 ir->valid.left * ps;
349
350 for( int x = 0; x < r->width; x++ ) {
351 const int ix = (int) X;
352 VipsPel *p = p0 + ix * ps;
353 const int sx = X * VIPS_TRANSFORM_SCALE * 2;
354 const int six = sx & (VIPS_TRANSFORM_SCALE * 2 - 1);
355 const int tx = (six + 1) >> 1;
356 const int *cxi = reduceh->matrixi[tx];
357 const double *cxf = reduceh->matrixf[tx];
358
359 switch( in->BandFmt ) {
360 case VIPS_FORMAT_UCHAR:
361 reduceh_unsigned_int_tab
362 <unsigned char, UCHAR_MAX>(
363 reduceh,
364 q, p, bands, cxi );
365 break;
366
367 case VIPS_FORMAT_CHAR:
368 reduceh_signed_int_tab
369 <signed char, SCHAR_MIN, SCHAR_MAX>(
370 reduceh,
371 q, p, bands, cxi );
372 break;
373
374 case VIPS_FORMAT_USHORT:
375 reduceh_unsigned_int_tab
376 <unsigned short, USHRT_MAX>(
377 reduceh,
378 q, p, bands, cxi );
379 break;
380
381 case VIPS_FORMAT_SHORT:
382 reduceh_signed_int_tab
383 <signed short, SHRT_MIN, SHRT_MAX>(
384 reduceh,
385 q, p, bands, cxi );
386 break;
387
388 case VIPS_FORMAT_UINT:
389 reduceh_unsigned_int32_tab
390 <unsigned int, INT_MAX>(
391 reduceh,
392 q, p, bands, cxf );
393 break;
394
395 case VIPS_FORMAT_INT:
396 reduceh_signed_int32_tab
397 <signed int, INT_MIN, INT_MAX>(
398 reduceh,
399 q, p, bands, cxf );
400 break;
401
402 case VIPS_FORMAT_FLOAT:
403 case VIPS_FORMAT_COMPLEX:
404 reduceh_float_tab<float>( reduceh,
405 q, p, bands, cxf );
406 break;
407
408 case VIPS_FORMAT_DOUBLE:
409 case VIPS_FORMAT_DPCOMPLEX:
410 reduceh_notab<double>( reduceh,
411 q, p, bands, X - ix );
412 break;
413
414 default:
415 g_assert_not_reached();
416 break;
417 }
418
419 X += reduceh->hshrink;
420 q += ps;
421 }
422 }
423
424 VIPS_GATE_STOP( "vips_reduceh_gen: work" );
425
426 VIPS_COUNT_PIXELS( out_region, "vips_reduceh_gen" );
427
428 return( 0 );
429 }
430
431 static int
vips_reduceh_build(VipsObject * object)432 vips_reduceh_build( VipsObject *object )
433 {
434 VipsObjectClass *object_class = VIPS_OBJECT_GET_CLASS( object );
435 VipsResample *resample = VIPS_RESAMPLE( object );
436 VipsReduceh *reduceh = (VipsReduceh *) object;
437 VipsImage **t = (VipsImage **)
438 vips_object_local_array( object, 2 );
439
440 VipsImage *in;
441 double width, extra_pixels;
442
443 if( VIPS_OBJECT_CLASS( vips_reduceh_parent_class )->build( object ) )
444 return( -1 );
445
446 in = resample->in;
447
448 if( reduceh->hshrink < 1 ) {
449 vips_error( object_class->nickname,
450 "%s", _( "reduce factors should be >= 1" ) );
451 return( -1 );
452 }
453
454 if( reduceh->hshrink == 1 )
455 return( vips_image_write( in, resample->out ) );
456
457 reduceh->n_point =
458 vips_reduce_get_points( reduceh->kernel, reduceh->hshrink );
459 g_info( "reduceh: %d point mask", reduceh->n_point );
460 if( reduceh->n_point > MAX_POINT ) {
461 vips_error( object_class->nickname,
462 "%s", _( "reduce factor too large" ) );
463 return( -1 );
464 }
465
466 /* Output size. We need to always round to nearest, so round(), not
467 * rint().
468 */
469 width = VIPS_ROUND_UINT(
470 (double) resample->in->Xsize / reduceh->hshrink );
471
472 /* How many pixels we are inventing in the input, -ve for
473 * discarding.
474 */
475 extra_pixels =
476 width * reduceh->hshrink - resample->in->Xsize;
477
478 /* If we are rounding down, we are not using some input
479 * pixels. We need to move the origin *inside* the input image
480 * by half that distance so that we discard pixels equally
481 * from left and right.
482 */
483 reduceh->hoffset = (1 + extra_pixels) / 2.0 - 1;
484
485 /* Build the tables of pre-computed coefficients.
486 */
487 for( int x = 0; x < VIPS_TRANSFORM_SCALE + 1; x++ ) {
488 reduceh->matrixf[x] =
489 VIPS_ARRAY( object, reduceh->n_point, double );
490 reduceh->matrixi[x] =
491 VIPS_ARRAY( object, reduceh->n_point, int );
492 if( !reduceh->matrixf[x] ||
493 !reduceh->matrixi[x] )
494 return( -1 );
495
496 vips_reduce_make_mask( reduceh->matrixf[x],
497 reduceh->kernel, reduceh->hshrink,
498 (float) x / VIPS_TRANSFORM_SCALE );
499
500 for( int i = 0; i < reduceh->n_point; i++ )
501 reduceh->matrixi[x][i] = reduceh->matrixf[x][i] *
502 VIPS_INTERPOLATE_SCALE;
503
504 #ifdef DEBUG
505 printf( "vips_reduceh_build: mask %d\n ", x );
506 for( int i = 0; i < reduceh->n_point; i++ )
507 printf( "%d ", reduceh->matrixi[x][i] );
508 printf( "\n" );
509 #endif /*DEBUG*/
510 }
511
512 /* Unpack for processing.
513 */
514 if( vips_image_decode( in, &t[0] ) )
515 return( -1 );
516 in = t[0];
517
518 /* Add new pixels around the input so we can interpolate at the edges.
519 */
520 if( vips_embed( in, &t[1],
521 VIPS_CEIL( reduceh->n_point / 2.0 ) - 1, 0,
522 in->Xsize + reduceh->n_point, in->Ysize,
523 "extend", VIPS_EXTEND_COPY,
524 (void *) NULL ) )
525 return( -1 );
526 in = t[1];
527
528 if( vips_image_pipelinev( resample->out,
529 VIPS_DEMAND_STYLE_THINSTRIP, in, (void *) NULL ) )
530 return( -1 );
531
532 /* Size output. We need to always round to nearest, so round(), not
533 * rint().
534 *
535 * Don't change xres/yres, leave that to the application layer. For
536 * example, vipsthumbnail knows the true reduce factor (including the
537 * fractional part), we just see the integer part here.
538 */
539 resample->out->Xsize = width;
540 if( resample->out->Xsize <= 0 ) {
541 vips_error( object_class->nickname,
542 "%s", _( "image has shrunk to nothing" ) );
543 return( -1 );
544 }
545
546 #ifdef DEBUG
547 printf( "vips_reduceh_build: reducing %d x %d image to %d x %d\n",
548 in->Xsize, in->Ysize,
549 resample->out->Xsize, resample->out->Ysize );
550 #endif /*DEBUG*/
551
552 if( vips_image_generate( resample->out,
553 vips_start_one, vips_reduceh_gen, vips_stop_one,
554 in, reduceh ) )
555 return( -1 );
556
557 vips_reorder_margin_hint( resample->out, reduceh->n_point );
558
559 return( 0 );
560 }
561
562 static void
vips_reduceh_class_init(VipsReducehClass * reduceh_class)563 vips_reduceh_class_init( VipsReducehClass *reduceh_class )
564 {
565 GObjectClass *gobject_class = G_OBJECT_CLASS( reduceh_class );
566 VipsObjectClass *vobject_class = VIPS_OBJECT_CLASS( reduceh_class );
567 VipsOperationClass *operation_class =
568 VIPS_OPERATION_CLASS( reduceh_class );
569
570 VIPS_DEBUG_MSG( "vips_reduceh_class_init\n" );
571
572 gobject_class->set_property = vips_object_set_property;
573 gobject_class->get_property = vips_object_get_property;
574
575 vobject_class->nickname = "reduceh";
576 vobject_class->description = _( "shrink an image horizontally" );
577 vobject_class->build = vips_reduceh_build;
578
579 operation_class->flags = VIPS_OPERATION_SEQUENTIAL;
580
581 VIPS_ARG_DOUBLE( reduceh_class, "hshrink", 3,
582 _( "Hshrink" ),
583 _( "Horizontal shrink factor" ),
584 VIPS_ARGUMENT_REQUIRED_INPUT,
585 G_STRUCT_OFFSET( VipsReduceh, hshrink ),
586 1, 1000000, 1 );
587
588 VIPS_ARG_ENUM( reduceh_class, "kernel", 4,
589 _( "Kernel" ),
590 _( "Resampling kernel" ),
591 VIPS_ARGUMENT_OPTIONAL_INPUT,
592 G_STRUCT_OFFSET( VipsReduceh, kernel ),
593 VIPS_TYPE_KERNEL, VIPS_KERNEL_LANCZOS3 );
594
595 /* Old name.
596 */
597 VIPS_ARG_DOUBLE( reduceh_class, "xshrink", 3,
598 _( "Xshrink" ),
599 _( "Horizontal shrink factor" ),
600 VIPS_ARGUMENT_REQUIRED_INPUT | VIPS_ARGUMENT_DEPRECATED,
601 G_STRUCT_OFFSET( VipsReduceh, hshrink ),
602 1, 1000000, 1 );
603
604 /* We used to let people pick centre or corner, but it's automatic now.
605 */
606 VIPS_ARG_BOOL( reduceh_class, "centre", 7,
607 _( "Centre" ),
608 _( "Use centre sampling convention" ),
609 VIPS_ARGUMENT_OPTIONAL_INPUT | VIPS_ARGUMENT_DEPRECATED,
610 G_STRUCT_OFFSET( VipsReduceh, centre ),
611 FALSE );
612
613 }
614
615 static void
vips_reduceh_init(VipsReduceh * reduceh)616 vips_reduceh_init( VipsReduceh *reduceh )
617 {
618 reduceh->kernel = VIPS_KERNEL_LANCZOS3;
619 }
620
621 /* See reduce.c for the doc comment.
622 */
623
624 int
vips_reduceh(VipsImage * in,VipsImage ** out,double hshrink,...)625 vips_reduceh( VipsImage *in, VipsImage **out, double hshrink, ... )
626 {
627 va_list ap;
628 int result;
629
630 va_start( ap, hshrink );
631 result = vips_call_split( "reduceh", ap, in, out, hshrink );
632 va_end( ap );
633
634 return( result );
635 }
636