1 /* complex.c ... various complex operations
2 *
3 * Copyright: 1990, N. Dessipris.
4 *
5 * Author: Nicos Dessipris
6 * Written on: 12/02/1990
7 * Modified on : 09/05/1990
8 * 15/6/93 JC
9 * - stupid stupid includes and externs fixed
10 * - I have been editing for 1 1/2 hours and I'm still drowning in
11 * rubbish extetrnshh
12 * 13/12/94 JC
13 * - modernised
14 * 9/7/02 JC
15 * - degree output, for consistency
16 * - slightly better behaviour in edge cases
17 * 27/1/10
18 * - modernised
19 * - gtk-doc
20 * 19/11/11
21 * - redo as a class
22 * 21/11/11
23 * - add vips_complexget()
24 * 29/9/15
25 * - return 0 for cross-product where one arg is zero
26 */
27
28 /*
29
30 Copyright (C) 1991-2005 The National Gallery
31
32 This library is free software; you can redistribute it and/or
33 modify it under the terms of the GNU Lesser General Public
34 License as published by the Free Software Foundation; either
35 version 2.1 of the License, or (at your option) any later version.
36
37 This library is distributed in the hope that it will be useful,
38 but WITHOUT ANY WARRANTY; without even the implied warranty of
39 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
40 Lesser General Public License for more details.
41
42 You should have received a copy of the GNU Lesser General Public
43 License along with this library; if not, write to the Free Software
44 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
45 02110-1301 USA
46
47 */
48
49 /*
50
51 These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
52
53 */
54
55 /*
56 #define DEBUG
57 */
58
59 #ifdef HAVE_CONFIG_H
60 #include <config.h>
61 #endif /*HAVE_CONFIG_H*/
62 #include <vips/intl.h>
63
64 #include <stdio.h>
65 #include <stdlib.h>
66 #include <math.h>
67
68 #include <vips/vips.h>
69
70 #include "unary.h"
71 #include "binary.h"
72
73 typedef struct _VipsComplex {
74 VipsUnary parent_instance;
75
76 VipsOperationComplex cmplx;
77
78 } VipsComplex;
79
80 typedef VipsUnaryClass VipsComplexClass;
81
82 G_DEFINE_TYPE( VipsComplex, vips_complex, VIPS_TYPE_UNARY );
83
84 #define LOOP( IN, OUT, OP ) { \
85 IN * restrict p = (IN *) in[0]; \
86 OUT * restrict q = (OUT *) out; \
87 \
88 for( x = 0; x < sz; x++ ) { \
89 OP( q, p[x], 0.0 ); \
90 \
91 q += 2; \
92 } \
93 }
94
95 #define CLOOP( IN, OUT, OP ) { \
96 IN * restrict p = (IN *) in[0]; \
97 OUT * restrict q = (OUT *) out; \
98 \
99 for( x = 0; x < sz; x++ ) { \
100 OP( q, p[0], p[1] ); \
101 \
102 p += 2; \
103 q += 2; \
104 } \
105 }
106
107 #define SWITCH( OP ) \
108 switch( vips_image_get_format( im ) ) { \
109 case VIPS_FORMAT_UCHAR: \
110 LOOP( unsigned char, float, OP ); break; \
111 case VIPS_FORMAT_CHAR: \
112 LOOP( signed char, float, OP ); break; \
113 case VIPS_FORMAT_USHORT: \
114 LOOP( unsigned short, float, OP ); break; \
115 case VIPS_FORMAT_SHORT: \
116 LOOP( signed short, float, OP ); break; \
117 case VIPS_FORMAT_UINT: \
118 LOOP( unsigned int, float, OP ); break; \
119 case VIPS_FORMAT_INT: \
120 LOOP( signed int, float, OP ); break; \
121 case VIPS_FORMAT_FLOAT: \
122 LOOP( float, float, OP ); break; \
123 case VIPS_FORMAT_DOUBLE: \
124 LOOP( double, double, OP ); break;\
125 case VIPS_FORMAT_COMPLEX: \
126 CLOOP( float, float, OP ); break; \
127 case VIPS_FORMAT_DPCOMPLEX: \
128 CLOOP( double, double, OP ); break;\
129 \
130 default: \
131 g_assert_not_reached(); \
132 }
133
134 static double
vips_complex_hypot(double a,double b)135 vips_complex_hypot( double a, double b )
136 {
137 double d;
138
139 /* hypot() is less sensitive to overflow. Use it if we can.
140 */
141 #ifdef HAVE_HYPOT
142 d = hypot( a, b );
143 #else
144 d = sqrt( a * a + b * b );
145 #endif
146
147 return( d );
148 }
149
150 static double
vips_complex_atan2(double a,double b)151 vips_complex_atan2( double a, double b )
152 {
153 double h;
154
155 /* atan2() is very slow, but is better behaved when a is near 0. Use
156 * it in preference when we can.
157 */
158 #ifdef HAVE_ATAN2
159 h = VIPS_DEG( atan2( b, a ) );
160 if( h < 0.0 )
161 h += 360;
162 #else
163 h = vips_col_ab2h( a, b );
164 #endif
165
166 return( h );
167 }
168
169 #define POLAR( Q, X, Y ) { \
170 double re = (X); \
171 double im = (Y); \
172 double am, ph; \
173 \
174 am = vips_complex_hypot( re, im ); \
175 ph = vips_complex_atan2( re, im ); \
176 \
177 Q[0] = am; \
178 Q[1] = ph; \
179 }
180
181 #define RECT( Q, X, Y ) { \
182 double am = (X); \
183 double ph = (Y); \
184 double re, im; \
185 \
186 re = am * cos( VIPS_RAD( ph ) ); \
187 im = am * sin( VIPS_RAD( ph ) ); \
188 \
189 Q[0] = re; \
190 Q[1] = im; \
191 }
192
193 #define CONJ( Q, X, Y ) { \
194 double re = (X); \
195 double im = (Y); \
196 \
197 im *= -1; \
198 \
199 Q[0] = re; \
200 Q[1] = im; \
201 }
202
203 static void
vips_complex_buffer(VipsArithmetic * arithmetic,VipsPel * out,VipsPel ** in,int width)204 vips_complex_buffer( VipsArithmetic *arithmetic,
205 VipsPel *out, VipsPel **in, int width )
206 {
207 VipsComplex *cmplx = (VipsComplex *) arithmetic;
208 VipsImage *im = arithmetic->ready[0];
209 const int sz = width * vips_image_get_bands( im );
210
211 int x;
212
213 switch( cmplx->cmplx ) {
214 case VIPS_OPERATION_COMPLEX_POLAR: SWITCH( POLAR ); break;
215 case VIPS_OPERATION_COMPLEX_RECT: SWITCH( RECT ); break;
216 case VIPS_OPERATION_COMPLEX_CONJ: SWITCH( CONJ ); break;
217
218 default:
219 g_assert_not_reached();
220 }
221 }
222
223 /* Save a bit of typing.
224 */
225 #define UC VIPS_FORMAT_UCHAR
226 #define C VIPS_FORMAT_CHAR
227 #define US VIPS_FORMAT_USHORT
228 #define S VIPS_FORMAT_SHORT
229 #define UI VIPS_FORMAT_UINT
230 #define I VIPS_FORMAT_INT
231 #define F VIPS_FORMAT_FLOAT
232 #define X VIPS_FORMAT_COMPLEX
233 #define D VIPS_FORMAT_DOUBLE
234 #define DX VIPS_FORMAT_DPCOMPLEX
235
236 static const VipsBandFormat vips_complex_format_table[10] = {
237 /* UC C US S UI I F X D DX */
238 X, X, X, X, X, X, X, X, DX, DX
239 };
240
241 static void
vips_complex_class_init(VipsComplexClass * class)242 vips_complex_class_init( VipsComplexClass *class )
243 {
244 GObjectClass *gobject_class = G_OBJECT_CLASS( class );
245 VipsObjectClass *object_class = (VipsObjectClass *) class;
246 VipsArithmeticClass *aclass = VIPS_ARITHMETIC_CLASS( class );
247
248 gobject_class->set_property = vips_object_set_property;
249 gobject_class->get_property = vips_object_get_property;
250
251 object_class->nickname = "complex";
252 object_class->description =
253 _( "perform a complex operation on an image" );
254
255 aclass->process_line = vips_complex_buffer;
256
257 vips_arithmetic_set_format_table( aclass, vips_complex_format_table );
258
259 VIPS_ARG_ENUM( class, "cmplx", 200,
260 _( "Operation" ),
261 _( "complex to perform" ),
262 VIPS_ARGUMENT_REQUIRED_INPUT,
263 G_STRUCT_OFFSET( VipsComplex, cmplx ),
264 VIPS_TYPE_OPERATION_COMPLEX, VIPS_OPERATION_COMPLEX_POLAR );
265 }
266
267 static void
vips_complex_init(VipsComplex * cmplx)268 vips_complex_init( VipsComplex *cmplx )
269 {
270 }
271
272 static int
vips_complexv(VipsImage * in,VipsImage ** out,VipsOperationComplex cmplx,va_list ap)273 vips_complexv( VipsImage *in, VipsImage **out,
274 VipsOperationComplex cmplx, va_list ap )
275 {
276 return( vips_call_split( "complex", ap, in, out, cmplx ) );
277 }
278
279 /**
280 * vips_complex: (method)
281 * @in: input #VipsImage
282 * @out: (out): output #VipsImage
283 * @cmplx: complex operation to perform
284 * @...: %NULL-terminated list of optional named arguments
285 *
286 * Perform various operations on complex images.
287 *
288 * Angles are expressed in degrees. The output type is complex unless the
289 * input is double or dpcomplex, in which case the output is dpcomplex.
290 *
291 * Returns: 0 on success, -1 on error
292 */
293 int
vips_complex(VipsImage * in,VipsImage ** out,VipsOperationComplex cmplx,...)294 vips_complex( VipsImage *in, VipsImage **out, VipsOperationComplex cmplx, ... )
295 {
296 va_list ap;
297 int result;
298
299 va_start( ap, cmplx );
300 result = vips_complexv( in, out, cmplx, ap );
301 va_end( ap );
302
303 return( result );
304 }
305
306 /**
307 * vips_polar: (method)
308 * @in: input #VipsImage
309 * @out: (out): output #VipsImage
310 * @...: %NULL-terminated list of optional named arguments
311 *
312 * Perform #VIPS_OPERATION_COMPLEX_POLAR on an image. See vips_complex().
313 *
314 * Returns: 0 on success, -1 on error
315 */
316 int
vips_polar(VipsImage * in,VipsImage ** out,...)317 vips_polar( VipsImage *in, VipsImage **out, ... )
318 {
319 va_list ap;
320 int result;
321
322 va_start( ap, out );
323 result = vips_complexv( in, out, VIPS_OPERATION_COMPLEX_POLAR, ap );
324 va_end( ap );
325
326 return( result );
327 }
328
329 /**
330 * vips_rect: (method)
331 * @in: input #VipsImage
332 * @out: (out): output #VipsImage
333 * @...: %NULL-terminated list of optional named arguments
334 *
335 * Perform #VIPS_OPERATION_COMPLEX_RECT on an image. See vips_complex().
336 *
337 * Returns: 0 on success, -1 on error
338 */
339 int
vips_rect(VipsImage * in,VipsImage ** out,...)340 vips_rect( VipsImage *in, VipsImage **out, ... )
341 {
342 va_list ap;
343 int result;
344
345 va_start( ap, out );
346 result = vips_complexv( in, out, VIPS_OPERATION_COMPLEX_RECT, ap );
347 va_end( ap );
348
349 return( result );
350 }
351
352 /**
353 * vips_conj: (method)
354 * @in: input #VipsImage
355 * @out: (out): output #VipsImage
356 * @...: %NULL-terminated list of optional named arguments
357 *
358 * Perform #VIPS_OPERATION_COMPLEX_CONJ on an image. See vips_complex().
359 *
360 * Returns: 0 on success, -1 on error
361 */
362 int
vips_conj(VipsImage * in,VipsImage ** out,...)363 vips_conj( VipsImage *in, VipsImage **out, ... )
364 {
365 va_list ap;
366 int result;
367
368 va_start( ap, out );
369 result = vips_complexv( in, out, VIPS_OPERATION_COMPLEX_CONJ, ap );
370 va_end( ap );
371
372 return( result );
373 }
374
375 typedef struct _VipsComplex2 {
376 VipsBinary parent_instance;
377
378 VipsOperationComplex2 cmplx;
379
380 } VipsComplex2;
381
382 typedef VipsBinaryClass VipsComplex2Class;
383
384 G_DEFINE_TYPE( VipsComplex2, vips_complex2, VIPS_TYPE_BINARY );
385
386 #define LOOP2( IN, OUT, OP ) { \
387 IN *p1 = (IN *) in[0]; \
388 IN *p2 = (IN *) in[1]; \
389 OUT *q = (OUT *) out; \
390 \
391 for( x = 0; x < sz; x++ ) { \
392 OP( q, p1[x], 0.0, p2[x], 0.0 ); \
393 \
394 q += 2; \
395 } \
396 }
397
398 #define CLOOP2( IN, OUT, OP ) { \
399 IN *p1 = (IN *) in[0]; \
400 IN *p2 = (IN *) in[1]; \
401 OUT *q = (OUT *) out; \
402 \
403 for( x = 0; x < sz; x++ ) { \
404 OP( q, p1[0], p1[1], p2[0], p2[1] ); \
405 \
406 p1 += 2; \
407 p2 += 2; \
408 q += 2; \
409 } \
410 }
411
412 #define SWITCH2( OP ) \
413 switch( vips_image_get_format( im ) ) { \
414 case VIPS_FORMAT_UCHAR: \
415 LOOP2( unsigned char, float, OP ); break; \
416 case VIPS_FORMAT_CHAR: \
417 LOOP2( signed char, float, OP ); break; \
418 case VIPS_FORMAT_USHORT: \
419 LOOP2( unsigned short, float, OP ); break; \
420 case VIPS_FORMAT_SHORT: \
421 LOOP2( signed short, float, OP ); break; \
422 case VIPS_FORMAT_UINT: \
423 LOOP2( unsigned int, float, OP ); break; \
424 case VIPS_FORMAT_INT: \
425 LOOP2( signed int, float, OP ); break; \
426 case VIPS_FORMAT_FLOAT: \
427 LOOP2( float, float, OP ); break; \
428 case VIPS_FORMAT_DOUBLE: \
429 LOOP2( double, double, OP ); break;\
430 case VIPS_FORMAT_COMPLEX: \
431 CLOOP2( float, float, OP ); break; \
432 case VIPS_FORMAT_DPCOMPLEX: \
433 CLOOP2( double, double, OP ); break;\
434 \
435 default: \
436 g_assert_not_reached(); \
437 }
438
439 /* There doesn't seem to be much difference in speed between these two methods
440 * (on an Athlon64), so I use the modulus argument version, since atan2() is
441 * in c89 but hypot() is c99.
442 *
443 * If you think that it might be faster on your platform, uncomment the
444 * following:
445 */
446 #define USE_MODARG_DIV
447
448 #ifdef USE_MODARG_DIV
449
450 #define CROSS( Q, X1, Y1, X2, Y2 ) { \
451 if( ((X1) == 0.0 && (Y1) == 0.0) || \
452 ((X2) == 0.0 && (Y2) == 0.0) ) { \
453 Q[0] = 0.0; \
454 Q[1] = 0.0; \
455 } \
456 else { \
457 double arg = atan2( X2, X1 ) - atan2( Y2, Y1 ); \
458 \
459 Q[0] = cos( arg ); \
460 Q[1] = sin( arg ); \
461 } \
462 }
463
464 #else /* USE_MODARG_DIV */
465
466 #define CROSS( Q, X1, Y1, X2, Y2 ) { \
467 if( ((X1) == 0.0 && (Y1) == 0.0) || \
468 ((X2) == 0.0 && (Y2) == 0.0) ) { \
469 Q[0] = 0.0; \
470 Q[1] = 0.0; \
471 } \
472 else if( ABS( Y1 ) > ABS( Y2 ) ) { \
473 double a = Y2 / Y1; \
474 double b = Y1 + Y2 * a; \
475 double re = (X1 + X2 * a) / b; \
476 double im = (X2 - X1 * a) / b; \
477 double mod = vips__hypot( re, im ); \
478 \
479 Q[0] = re / mod; \
480 Q[1] = im / mod; \
481 } \
482 else { \
483 double a = Y1 / Y2; \
484 double b = Y2 + Y1 * a; \
485 double re = (X1 * a + X2) / b; \
486 double im = (X2 * a - X1) / b; \
487 double mod = vips__hypot( re, im ); \
488 \
489 Q[0] = re / mod; \
490 Q[1] = im / mod; \
491 } \
492 }
493
494 #endif /* USE_MODARG_DIV */
495
496 static void
vips_complex2_buffer(VipsArithmetic * arithmetic,VipsPel * out,VipsPel ** in,int width)497 vips_complex2_buffer( VipsArithmetic *arithmetic,
498 VipsPel *out, VipsPel **in, int width )
499 {
500 VipsComplex2 *cmplx = (VipsComplex2 *) arithmetic;
501 VipsImage *im = arithmetic->ready[0];
502 const int sz = width * vips_image_get_bands( im );
503
504 int x;
505
506 switch( cmplx->cmplx ) {
507 case VIPS_OPERATION_COMPLEX2_CROSS_PHASE:
508 SWITCH2( CROSS );
509 break;
510
511 default:
512 g_assert_not_reached();
513 }
514 }
515
516 /* Save a bit of typing.
517 */
518 #define UC VIPS_FORMAT_UCHAR
519 #define C VIPS_FORMAT_CHAR
520 #define US VIPS_FORMAT_USHORT
521 #define S VIPS_FORMAT_SHORT
522 #define UI VIPS_FORMAT_UINT
523 #define I VIPS_FORMAT_INT
524 #define F VIPS_FORMAT_FLOAT
525 #define X VIPS_FORMAT_COMPLEX
526 #define D VIPS_FORMAT_DOUBLE
527 #define DX VIPS_FORMAT_DPCOMPLEX
528
529 static const VipsBandFormat vips_complex2_format_table[10] = {
530 /* UC C US S UI I F X D DX */
531 X, X, X, X, X, X, X, X, DX, DX
532 };
533
534 static void
vips_complex2_class_init(VipsComplex2Class * class)535 vips_complex2_class_init( VipsComplex2Class *class )
536 {
537 GObjectClass *gobject_class = G_OBJECT_CLASS( class );
538 VipsObjectClass *object_class = (VipsObjectClass *) class;
539 VipsArithmeticClass *aclass = VIPS_ARITHMETIC_CLASS( class );
540
541 gobject_class->set_property = vips_object_set_property;
542 gobject_class->get_property = vips_object_get_property;
543
544 object_class->nickname = "complex2";
545 object_class->description =
546 _( "complex binary operations on two images" );
547
548 aclass->process_line = vips_complex2_buffer;
549
550 vips_arithmetic_set_format_table( aclass, vips_complex2_format_table );
551
552 VIPS_ARG_ENUM( class, "cmplx", 200,
553 _( "Operation" ),
554 _( "binary complex operation to perform" ),
555 VIPS_ARGUMENT_REQUIRED_INPUT,
556 G_STRUCT_OFFSET( VipsComplex2, cmplx ),
557 VIPS_TYPE_OPERATION_COMPLEX2,
558 VIPS_OPERATION_COMPLEX2_CROSS_PHASE );
559 }
560
561 static void
vips_complex2_init(VipsComplex2 * cmplx)562 vips_complex2_init( VipsComplex2 *cmplx )
563 {
564 }
565
566 static int
vips_complex2v(VipsImage * left,VipsImage * right,VipsImage ** out,VipsOperationComplex2 cmplx,va_list ap)567 vips_complex2v( VipsImage *left, VipsImage *right, VipsImage **out,
568 VipsOperationComplex2 cmplx, va_list ap )
569 {
570 return( vips_call_split( "complex2", ap, left, right, out, cmplx ) );
571 }
572
573 /**
574 * vips_complex2:
575 * @left: input #VipsImage
576 * @right: input #VipsImage
577 * @out: (out): output #VipsImage
578 * @cmplx: complex2 operation to perform
579 * @...: %NULL-terminated list of optional named arguments
580 *
581 * Perform various binary operations on complex images.
582 *
583 * Angles are expressed in degrees. The output type is complex unless the
584 * input is double or dpcomplex, in which case the output is dpcomplex.
585 *
586 * Returns: 0 on success, -1 on error
587 */
588 int
vips_complex2(VipsImage * left,VipsImage * right,VipsImage ** out,VipsOperationComplex2 cmplx,...)589 vips_complex2( VipsImage *left, VipsImage *right, VipsImage **out,
590 VipsOperationComplex2 cmplx, ... )
591 {
592 va_list ap;
593 int result;
594
595 va_start( ap, cmplx );
596 result = vips_complex2v( left, right, out, cmplx, ap );
597 va_end( ap );
598
599 return( result );
600 }
601
602 /**
603 * vips_cross_phase:
604 * @left: input #VipsImage
605 * @right: input #VipsImage
606 * @out: (out): output #VipsImage
607 * @...: %NULL-terminated list of optional named arguments
608 *
609 * Perform #VIPS_OPERATION_COMPLEX2_CROSS_PHASE on an image.
610 * See vips_complex2().
611 *
612 * Returns: 0 on success, -1 on error
613 */
614 int
vips_cross_phase(VipsImage * left,VipsImage * right,VipsImage ** out,...)615 vips_cross_phase( VipsImage *left, VipsImage *right, VipsImage **out, ... )
616 {
617 va_list ap;
618 int result;
619
620 va_start( ap, out );
621 result = vips_complex2v( left, right, out,
622 VIPS_OPERATION_COMPLEX2_CROSS_PHASE, ap );
623 va_end( ap );
624
625 return( result );
626 }
627
628 typedef struct _VipsComplexget {
629 VipsUnary parent_instance;
630
631 VipsOperationComplexget get;
632
633 } VipsComplexget;
634
635 typedef VipsUnaryClass VipsComplexgetClass;
636
637 G_DEFINE_TYPE( VipsComplexget, vips_complexget, VIPS_TYPE_UNARY );
638
639 static int
vips_complexget_build(VipsObject * object)640 vips_complexget_build( VipsObject *object )
641 {
642 VipsUnary *unary = (VipsUnary *) object;
643 VipsComplexget *complexget = (VipsComplexget *) object;
644
645 if( unary->in &&
646 !vips_band_format_iscomplex( unary->in->BandFmt ) &&
647 complexget->get == VIPS_OPERATION_COMPLEXGET_REAL )
648 return( vips_unary_copy( unary ) );
649
650 if( VIPS_OBJECT_CLASS( vips_complexget_parent_class )->build( object ) )
651 return( -1 );
652
653 return( 0 );
654 }
655
656 #define GETLOOP( TYPE, OP ) { \
657 TYPE *p __attribute__ ((unused)) = (TYPE *) in[0]; \
658 TYPE *q = (TYPE *) out; \
659 \
660 for( x = 0; x < sz; x++ ) { \
661 OP( q[x], p[x], 0.0 ); \
662 } \
663 }
664
665 #define CGETLOOP( TYPE, OP ) { \
666 TYPE *p __attribute__ ((unused)) = (TYPE *) in[0]; \
667 TYPE *q = (TYPE *) out; \
668 \
669 for( x = 0; x < sz; x++ ) { \
670 OP( q[x], p[0], p[1] ); \
671 \
672 p += 2; \
673 } \
674 }
675
676 #define GETSWITCH( OP ) \
677 switch( vips_image_get_format( im ) ) { \
678 case VIPS_FORMAT_UCHAR: \
679 GETLOOP( unsigned char, OP ); break; \
680 case VIPS_FORMAT_CHAR: \
681 GETLOOP( signed char, OP ); break; \
682 case VIPS_FORMAT_USHORT: \
683 GETLOOP( unsigned short, OP ); break; \
684 case VIPS_FORMAT_SHORT: \
685 GETLOOP( signed short, OP ); break; \
686 case VIPS_FORMAT_UINT: \
687 GETLOOP( unsigned int, OP ); break; \
688 case VIPS_FORMAT_INT: \
689 GETLOOP( signed int, OP ); break; \
690 case VIPS_FORMAT_FLOAT: \
691 GETLOOP( float, OP ); break; \
692 case VIPS_FORMAT_DOUBLE: \
693 GETLOOP( double, OP ); break;\
694 case VIPS_FORMAT_COMPLEX: \
695 CGETLOOP( float, OP ); break; \
696 case VIPS_FORMAT_DPCOMPLEX: \
697 CGETLOOP( double, OP ); break;\
698 \
699 default: \
700 g_assert_not_reached(); \
701 }
702
703 #define REAL( Q, X, Y ) { \
704 Q = X; \
705 }
706
707 #define IMAG( Q, X, Y ) { \
708 Q = Y; \
709 }
710
711 static void
vips_complexget_buffer(VipsArithmetic * arithmetic,VipsPel * out,VipsPel ** in,int width)712 vips_complexget_buffer( VipsArithmetic *arithmetic,
713 VipsPel *out, VipsPel **in, int width )
714 {
715 VipsComplexget *complexget = (VipsComplexget *) arithmetic;
716 VipsImage *im = arithmetic->ready[0];
717 const int sz = width * vips_image_get_bands( im );
718
719 int x;
720
721 switch( complexget->get ) {
722 case VIPS_OPERATION_COMPLEXGET_REAL: GETSWITCH( REAL ); break;
723 case VIPS_OPERATION_COMPLEXGET_IMAG: GETSWITCH( IMAG ); break;
724
725 default:
726 g_assert_not_reached();
727 }
728 }
729
730 /* Save a bit of typing.
731 */
732 #define UC VIPS_FORMAT_UCHAR
733 #define C VIPS_FORMAT_CHAR
734 #define US VIPS_FORMAT_USHORT
735 #define S VIPS_FORMAT_SHORT
736 #define UI VIPS_FORMAT_UINT
737 #define I VIPS_FORMAT_INT
738 #define F VIPS_FORMAT_FLOAT
739 #define X VIPS_FORMAT_COMPLEX
740 #define D VIPS_FORMAT_DOUBLE
741 #define DX VIPS_FORMAT_DPCOMPLEX
742
743 static const VipsBandFormat vips_complexget_format_table[10] = {
744 /* UC C US S UI I F X D DX */
745 UC, C, US, S, UI, I, F, F, D, D
746 };
747
748 static void
vips_complexget_class_init(VipsComplexgetClass * class)749 vips_complexget_class_init( VipsComplexgetClass *class )
750 {
751 GObjectClass *gobject_class = G_OBJECT_CLASS( class );
752 VipsObjectClass *object_class = (VipsObjectClass *) class;
753 VipsArithmeticClass *aclass = VIPS_ARITHMETIC_CLASS( class );
754
755 gobject_class->set_property = vips_object_set_property;
756 gobject_class->get_property = vips_object_get_property;
757
758 object_class->nickname = "complexget";
759 object_class->description = _( "get a component from a complex image" );
760 object_class->build = vips_complexget_build;
761
762 aclass->process_line = vips_complexget_buffer;
763
764 vips_arithmetic_set_format_table( aclass,
765 vips_complexget_format_table );
766
767 VIPS_ARG_ENUM( class, "get", 200,
768 _( "Operation" ),
769 _( "complex to perform" ),
770 VIPS_ARGUMENT_REQUIRED_INPUT,
771 G_STRUCT_OFFSET( VipsComplexget, get ),
772 VIPS_TYPE_OPERATION_COMPLEXGET,
773 VIPS_OPERATION_COMPLEXGET_REAL );
774 }
775
776 static void
vips_complexget_init(VipsComplexget * complexget)777 vips_complexget_init( VipsComplexget *complexget )
778 {
779 }
780
781 static int
vips_complexgetv(VipsImage * in,VipsImage ** out,VipsOperationComplexget get,va_list ap)782 vips_complexgetv( VipsImage *in, VipsImage **out,
783 VipsOperationComplexget get, va_list ap )
784 {
785 return( vips_call_split( "complexget", ap, in, out, get ) );
786 }
787
788 /**
789 * vips_complexget: (method)
790 * @in: input #VipsImage
791 * @out: (out): output #VipsImage
792 * @get: complex operation to perform
793 * @...: %NULL-terminated list of optional named arguments
794 *
795 * Get components of complex images.
796 *
797 * The output type is the same as the input type, except #VIPS_FORMAT_COMPLEX
798 * becomes #VIPS_FORMAT_FLOAT and #VIPS_FORMAT_DPCOMPLEX becomes
799 * #VIPS_FORMAT_DOUBLE.
800 *
801 * Returns: 0 on success, -1 on error
802 */
803 int
vips_complexget(VipsImage * in,VipsImage ** out,VipsOperationComplexget get,...)804 vips_complexget( VipsImage *in, VipsImage **out,
805 VipsOperationComplexget get, ... )
806 {
807 va_list ap;
808 int result;
809
810 va_start( ap, get );
811 result = vips_complexgetv( in, out, get, ap );
812 va_end( ap );
813
814 return( result );
815 }
816
817 /**
818 * vips_real: (method)
819 * @in: input #VipsImage
820 * @out: (out): output #VipsImage
821 * @...: %NULL-terminated list of optional named arguments
822 *
823 * Perform #VIPS_OPERATION_COMPLEXGET_REAL on an image. See vips_complexget().
824 *
825 * Returns: 0 on success, -1 on error
826 */
827 int
vips_real(VipsImage * in,VipsImage ** out,...)828 vips_real( VipsImage *in, VipsImage **out, ... )
829 {
830 va_list ap;
831 int result;
832
833 va_start( ap, out );
834 result = vips_complexgetv( in, out,
835 VIPS_OPERATION_COMPLEXGET_REAL, ap );
836 va_end( ap );
837
838 return( result );
839 }
840
841 /**
842 * vips_imag: (method)
843 * @in: input #VipsImage
844 * @out: (out): output #VipsImage
845 * @...: %NULL-terminated list of optional named arguments
846 *
847 * Perform #VIPS_OPERATION_COMPLEXGET_IMAG on an image. See vips_complexget().
848 *
849 * Returns: 0 on success, -1 on error
850 */
851 int
vips_imag(VipsImage * in,VipsImage ** out,...)852 vips_imag( VipsImage *in, VipsImage **out, ... )
853 {
854 va_list ap;
855 int result;
856
857 va_start( ap, out );
858 result = vips_complexgetv( in, out,
859 VIPS_OPERATION_COMPLEXGET_IMAG, ap );
860 va_end( ap );
861
862 return( result );
863 }
864
865 typedef VipsBinary VipsComplexform;
866 typedef VipsBinaryClass VipsComplexformClass;
867
868 G_DEFINE_TYPE( VipsComplexform, vips_complexform, VIPS_TYPE_BINARY );
869
870 static int
vips_complexform_build(VipsObject * object)871 vips_complexform_build( VipsObject *object )
872 {
873 VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( object );
874 VipsBinary *binary = (VipsBinary *) object;
875
876 if( binary->left &&
877 vips_check_noncomplex( class->nickname, binary->left ) )
878 return( -1 );
879 if( binary->right &&
880 vips_check_noncomplex( class->nickname, binary->right ) )
881 return( -1 );
882
883 if( VIPS_OBJECT_CLASS( vips_complexform_parent_class )->
884 build( object ) )
885 return( -1 );
886
887 return( 0 );
888 }
889
890 #define CFORM( IN, OUT ) { \
891 IN *left = (IN *) in[0]; \
892 IN *right = (IN *) in[1]; \
893 OUT *q = (OUT *) out; \
894 \
895 for( x = 0; x < sz; x++ ) { \
896 q[0] = left[x]; \
897 q[1] = right[x]; \
898 \
899 q += 2; \
900 } \
901 }
902
903 static void
vips_complexform_buffer(VipsArithmetic * arithmetic,VipsPel * out,VipsPel ** in,int width)904 vips_complexform_buffer( VipsArithmetic *arithmetic,
905 VipsPel *out, VipsPel **in, int width )
906 {
907 VipsImage *im = arithmetic->ready[0];
908 const int sz = width * vips_image_get_bands( im );
909
910 int x;
911
912 /* Keep types here in sync with bandfmt_complexform[]
913 * below.
914 */
915 switch( vips_image_get_format( im ) ) {
916 case VIPS_FORMAT_UCHAR:
917 CFORM( unsigned char, float ); break;
918 case VIPS_FORMAT_CHAR:
919 CFORM( signed char, float ); break;
920 case VIPS_FORMAT_USHORT:
921 CFORM( unsigned short, float ); break;
922 case VIPS_FORMAT_SHORT:
923 CFORM( signed short, float ); break;
924 case VIPS_FORMAT_UINT:
925 CFORM( unsigned int, float ); break;
926 case VIPS_FORMAT_INT:
927 CFORM( signed int, float ); break;
928 case VIPS_FORMAT_FLOAT:
929 CFORM( float, float ); break;
930 case VIPS_FORMAT_DOUBLE:
931 CFORM( double, double ); break;
932
933 default:
934 g_assert_not_reached();
935 }
936 }
937
938 /* Save a bit of typing.
939 */
940 #define UC VIPS_FORMAT_UCHAR
941 #define C VIPS_FORMAT_CHAR
942 #define US VIPS_FORMAT_USHORT
943 #define S VIPS_FORMAT_SHORT
944 #define UI VIPS_FORMAT_UINT
945 #define I VIPS_FORMAT_INT
946 #define F VIPS_FORMAT_FLOAT
947 #define X VIPS_FORMAT_COMPLEX
948 #define D VIPS_FORMAT_DOUBLE
949 #define DX VIPS_FORMAT_DPCOMPLEX
950
951 /* Type promotion for form complex. Sign and value preserving. Make sure
952 * these match the case statement in complexform_buffer() above.
953 */
954 static int vips_complexform_format_table[10] = {
955 /* UC C US S UI I F X D DX */
956 X, X, X, X, X, X, X, X, DX,DX
957 };
958
959 static void
vips_complexform_class_init(VipsComplexformClass * class)960 vips_complexform_class_init( VipsComplexformClass *class )
961 {
962 VipsObjectClass *object_class = (VipsObjectClass *) class;
963 VipsArithmeticClass *aclass = VIPS_ARITHMETIC_CLASS( class );
964
965 object_class->nickname = "complexform";
966 object_class->description =
967 _( "form a complex image from two real images" );
968 object_class->build = vips_complexform_build;
969
970 aclass->process_line = vips_complexform_buffer;
971
972 vips_arithmetic_set_format_table( aclass,
973 vips_complexform_format_table );
974 }
975
976 static void
vips_complexform_init(VipsComplexform * complexform)977 vips_complexform_init( VipsComplexform *complexform )
978 {
979 }
980
981 /**
982 * vips_complexform:
983 * @left: input image
984 * @right: input image
985 * @out: (out): output image
986 * @...: %NULL-terminated list of optional named arguments
987 *
988 * Compose two real images to make a complex image. If either @left or @right
989 * are #VIPS_FORMAT_DOUBLE, @out is #VIPS_FORMAT_DPCOMPLEX. Otherwise @out
990 * is #VIPS_FORMAT_COMPLEX. @left becomes the real component of @out and
991 * @right the imaginary.
992 *
993 * If the number of bands differs, one of the images
994 * must have one band. In this case, an n-band image is formed from the
995 * one-band image by joining n copies of the one-band image together, and then
996 * the two n-band images are operated upon.
997 *
998 * See also: vips_complexget().
999 *
1000 * Returns: 0 on success, -1 on error
1001 */
1002 int
vips_complexform(VipsImage * left,VipsImage * right,VipsImage ** out,...)1003 vips_complexform( VipsImage *left, VipsImage *right, VipsImage **out, ... )
1004 {
1005 va_list ap;
1006 int result;
1007
1008 va_start( ap, out );
1009 result = vips_call_split( "complexform", ap, left, right, out );
1010 va_end( ap );
1011
1012 return( result );
1013 }
1014
1015
1016
1017
1018
1019