1 /* Divide two images
2 *
3 * Copyright: 1990, N. Dessipris.
4 *
5 * Author: Nicos Dessipris
6 * Written on: 02/05/1990
7 * Modified on:
8 * 29/4/93 JC
9 * - now works for partial images
10 * 1/7/93 JC
11 * - adapted for partial v2
12 * - ANSIfied
13 * 19/10/93 JC
14 * - coredump-inducing bug in complex*complex fixed
15 * 13/12/93
16 * - char * short bug fixed
17 * 12/6/95 JC
18 * - new im_multiply adapted to make new im_divide
19 * 27/9/04
20 * - updated for 1 band $op n band image -> n band image case
21 * 8/12/06
22 * - add liboil support
23 * 18/8/08
24 * - revise upcasting system
25 * - add gtkdoc comments
26 * 31/7/10
27 * - remove liboil support
28 * - avoid /0
29 * 6/11/11
30 * - rewrite as a class
31 * 22/2/12
32 * - avoid /0 for complex as well
33 * 6/4/12
34 * - fixed switch cases
35 * - fixed int operands with <1 result
36 */
37
38 /*
39
40 Copyright (C) 1991-2005 The National Gallery
41
42 This library is free software; you can redistribute it and/or
43 modify it under the terms of the GNU Lesser General Public
44 License as published by the Free Software Foundation; either
45 version 2.1 of the License, or (at your option) any later version.
46
47 This library is distributed in the hope that it will be useful,
48 but WITHOUT ANY WARRANTY; without even the implied warranty of
49 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
50 Lesser General Public License for more details.
51
52 You should have received a copy of the GNU Lesser General Public
53 License along with this library; if not, write to the Free Software
54 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
55 02110-1301 USA
56
57 */
58
59 /*
60
61 These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
62
63 */
64
65 /*
66 #define DEBUG
67 */
68
69 #ifdef HAVE_CONFIG_H
70 #include <config.h>
71 #endif /*HAVE_CONFIG_H*/
72 #include <vips/intl.h>
73
74 #include <stdio.h>
75 #include <stdlib.h>
76 #include <math.h>
77
78 #include <vips/vips.h>
79
80 #include "binary.h"
81
82 typedef VipsBinary VipsDivide;
83 typedef VipsBinaryClass VipsDivideClass;
84
85 G_DEFINE_TYPE( VipsDivide, vips_divide, VIPS_TYPE_BINARY );
86
87 /* Complex divide.
88 */
89 #ifdef USE_MODARG_DIV
90
91 /* This is going to be much slower */
92 #define CLOOP( TYPE ) { \
93 TYPE * restrict left = (TYPE *) in[0]; \
94 TYPE * restrict right = (TYPE *) in[1]; \
95 TYPE * restrict q = (TYPE *) out; \
96 int i; \
97 \
98 for( i = 0; i < sz; i++ ) { \
99 if( right[0] == 0.0 && \
100 right[1] == 0.0 ) { \
101 q[0] = 0.0; \
102 q[1] = 0.0; \
103 } \
104 else { \
105 double arg = atan2( left[1], left[0] ) - \
106 atan2( right[1], right[0] ); \
107 double mod = hypot( left[1], left[0] ) / \
108 hypot( right[1], right[0] ); \
109 \
110 q[0] = mod * cos( arg ); \
111 q[1] = mod * sin( arg ); \
112 } \
113 \
114 left += 2; \
115 right += 2; \
116 q += 2; \
117 } \
118 }
119
120 #else /* USE_MODARG_DIV */
121
122 #define CLOOP( TYPE ) { \
123 TYPE * restrict left = (TYPE *) in[0]; \
124 TYPE * restrict right = (TYPE *) in[1]; \
125 TYPE * restrict q = (TYPE *) out; \
126 int i; \
127 \
128 for( i = 0; i < sz; i++ ) { \
129 if( right[0] == 0.0 && \
130 right[1] == 0.0 ) { \
131 q[0] = 0.0; \
132 q[1] = 0.0; \
133 } \
134 else if( VIPS_FABS( right[0] ) > VIPS_FABS( right[1] ) ) { \
135 double a = right[1] / right[0]; \
136 double b = right[0] + right[1] * a; \
137 \
138 q[0] = (left[0] + left[1] * a) / b; \
139 q[1] = (left[1] - left[0] * a) / b; \
140 } \
141 else { \
142 double a = right[0] / right[1]; \
143 double b = right[1] + right[0] * a; \
144 \
145 q[0] = (left[0] * a + left[1]) / b; \
146 q[1] = (left[1] * a - left[0]) / b; \
147 } \
148 \
149 left += 2; \
150 right += 2; \
151 q += 2; \
152 } \
153 }
154
155 #endif /* USE_MODARG_DIV */
156
157 /* Real divide. Cast in to OUT before divide so we work for float output.
158 */
159 #define RLOOP( IN, OUT ) { \
160 IN * restrict left = (IN *) in[0]; \
161 IN * restrict right = (IN *) in[1]; \
162 OUT * restrict q = (OUT *) out; \
163 \
164 for( x = 0; x < sz; x++ ) \
165 q[x] = right[x] == 0 ? 0 : (OUT) left[x] / (OUT) right[x]; \
166 }
167
168 static void
vips_divide_buffer(VipsArithmetic * arithmetic,VipsPel * out,VipsPel ** in,int width)169 vips_divide_buffer( VipsArithmetic *arithmetic,
170 VipsPel *out, VipsPel **in, int width )
171 {
172 VipsImage *im = arithmetic->ready[0];
173 const int sz = width * vips_image_get_bands( im );
174
175 int x;
176
177 /* Keep types here in sync with vips_divide_format_table[]
178 * below.
179 */
180 switch( vips_image_get_format( im ) ) {
181 case VIPS_FORMAT_CHAR: RLOOP( signed char, float ); break;
182 case VIPS_FORMAT_UCHAR: RLOOP( unsigned char, float ); break;
183 case VIPS_FORMAT_SHORT: RLOOP( signed short, float ); break;
184 case VIPS_FORMAT_USHORT: RLOOP( unsigned short, float ); break;
185 case VIPS_FORMAT_INT: RLOOP( signed int, float ); break;
186 case VIPS_FORMAT_UINT: RLOOP( unsigned int, float ); break;
187 case VIPS_FORMAT_FLOAT: RLOOP( float, float ); break;
188 case VIPS_FORMAT_DOUBLE: RLOOP( double, double ); break;
189 case VIPS_FORMAT_COMPLEX: CLOOP( float ); break;
190 case VIPS_FORMAT_DPCOMPLEX: CLOOP( double ); break;
191
192 default:
193 g_assert_not_reached();
194 }
195 }
196
197 /* Save a bit of typing.
198 */
199 #define UC VIPS_FORMAT_UCHAR
200 #define C VIPS_FORMAT_CHAR
201 #define US VIPS_FORMAT_USHORT
202 #define S VIPS_FORMAT_SHORT
203 #define UI VIPS_FORMAT_UINT
204 #define I VIPS_FORMAT_INT
205 #define F VIPS_FORMAT_FLOAT
206 #define X VIPS_FORMAT_COMPLEX
207 #define D VIPS_FORMAT_DOUBLE
208 #define DX VIPS_FORMAT_DPCOMPLEX
209
210 /* Type promotion for division. Sign and value preserving. Make sure
211 * these match the case statement in divide_buffer() above.
212 */
213 static int vips_divide_format_table[10] = {
214 /* UC C US S UI I F X D DX */
215 F, F, F, F, F, F, F, X, D, DX
216 };
217
218 static void
vips_divide_class_init(VipsDivideClass * class)219 vips_divide_class_init( VipsDivideClass *class )
220 {
221 VipsObjectClass *object_class = (VipsObjectClass *) class;
222 VipsArithmeticClass *aclass = VIPS_ARITHMETIC_CLASS( class );
223
224 object_class->nickname = "divide";
225 object_class->description = _( "divide two images" );
226
227 aclass->process_line = vips_divide_buffer;
228
229 vips_arithmetic_set_format_table( aclass, vips_divide_format_table );
230 }
231
232 static void
vips_divide_init(VipsDivide * divide)233 vips_divide_init( VipsDivide *divide )
234 {
235 }
236
237 /**
238 * vips_divide:
239 * @left: input image
240 * @right: input image
241 * @out: (out): output image
242 * @...: %NULL-terminated list of optional named arguments
243 *
244 * This operation calculates @in1 / @in2 and writes the result to @out. If any
245 * pixels in @in2 are zero, the corresponding pixel in @out is also zero.
246 *
247 * If the images differ in size, the smaller image is enlarged to match the
248 * larger by adding zero pixels along the bottom and right.
249 *
250 * If the number of bands differs, one of the images
251 * must have one band. In this case, an n-band image is formed from the
252 * one-band image by joining n copies of the one-band image together, and then
253 * the two n-band images are operated upon.
254 *
255 * The two input images are cast up to the smallest common format (see table
256 * Smallest common format in
257 * <link linkend="libvips-arithmetic">arithmetic</link>), then the
258 * following table is used to determine the output type:
259 *
260 * <table>
261 * <title>vips_divide() type promotion</title>
262 * <tgroup cols='2' align='left' colsep='1' rowsep='1'>
263 * <thead>
264 * <row>
265 * <entry>input type</entry>
266 * <entry>output type</entry>
267 * </row>
268 * </thead>
269 * <tbody>
270 * <row>
271 * <entry>uchar</entry>
272 * <entry>float</entry>
273 * </row>
274 * <row>
275 * <entry>char</entry>
276 * <entry>float</entry>
277 * </row>
278 * <row>
279 * <entry>ushort</entry>
280 * <entry>float</entry>
281 * </row>
282 * <row>
283 * <entry>short</entry>
284 * <entry>float</entry>
285 * </row>
286 * <row>
287 * <entry>uint</entry>
288 * <entry>float</entry>
289 * </row>
290 * <row>
291 * <entry>int</entry>
292 * <entry>float</entry>
293 * </row>
294 * <row>
295 * <entry>float</entry>
296 * <entry>float</entry>
297 * </row>
298 * <row>
299 * <entry>double</entry>
300 * <entry>double</entry>
301 * </row>
302 * <row>
303 * <entry>complex</entry>
304 * <entry>complex</entry>
305 * </row>
306 * <row>
307 * <entry>double complex</entry>
308 * <entry>double complex</entry>
309 * </row>
310 * </tbody>
311 * </tgroup>
312 * </table>
313 *
314 * In other words, the output type is just large enough to hold the whole
315 * range of possible values.
316 *
317 * See also: vips_multiply(), vips_linear(), vips_pow().
318 *
319 * Returns: 0 on success, -1 on error
320 */
321 int
vips_divide(VipsImage * left,VipsImage * right,VipsImage ** out,...)322 vips_divide( VipsImage *left, VipsImage *right, VipsImage **out, ... )
323 {
324 va_list ap;
325 int result;
326
327 va_start( ap, out );
328 result = vips_call_split( "divide", ap, left, right, out );
329 va_end( ap );
330
331 return( result );
332 }
333