Source code

001/*******************************************************************************
002 * This software is provided as a supplement to the authors' textbooks on digital
003 * image processing published by Springer-Verlag in various languages and editions.
004 * Permission to use and distribute this software is granted under the BSD 2-Clause
005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause).
006 * Copyright (c) 2006-2023 Wilhelm Burger, Mark J. Burge. All rights reserved.
007 * Visit https://imagingbook.com for additional details.
008 ******************************************************************************/
009
010package imagingbook.common.math;
011
012import static org.apache.commons.math3.util.FastMath.abs;
013import static org.apache.commons.math3.util.FastMath.atan2;
014import static org.apache.commons.math3.util.FastMath.cos;
015import static org.apache.commons.math3.util.FastMath.floor;
016import static org.apache.commons.math3.util.FastMath.floorMod;
017import static org.apache.commons.math3.util.FastMath.hypot;
018import static org.apache.commons.math3.util.FastMath.sin;
019
020/**
021 * This class defines static methods implementing arithmetic operations and predicates.
022 * 
023 * @author WB
024 *
025 */
026public abstract class Arithmetic {
027        
028        private Arithmetic() {}
029        
030        // machine accuracy for IEEE 754 float/double;
031        /** Default tolerance used for comparing {@code float} quantities. */
032        public static final float EPSILON_FLOAT = 1e-7f;        // 1.19 x 10^-7
033        /** Default tolerance used for comparing {@code double} quantities. */
034        public static final double EPSILON_DOUBLE = 2e-16;      // 2.22 x 10^-16
035
036        /**
037         * Returns the square of its argument.
038         * @param x argument
039         * @return square of argument
040         */
041        public static int sqr(int x) {
042                return x * x;
043        }
044        
045        /**
046         * Returns the square of its argument.
047         * @param x argument
048         * @return square of argument
049         */
050        public static long sqr(long x) {
051                return x * x;
052        }
053        
054        /**
055         * Returns the square of its argument.
056         * @param x argument
057         * @return square of argument
058         */
059        public static float sqr(float x) {
060                return x * x;
061        }
062        
063        /**
064         * Returns the square of its argument.
065         * @param x argument
066         * @return square of argument
067         */
068        public static double sqr(double x) {
069                return x * x;
070        }
071        
072        /**
073         * Returns the radius of the Cartesian point 
074         * @param x x-component
075         * @param y y-component
076         * @return the radius
077         */
078        public static double radius(double x, double y) {
079                return hypot(x, y);
080        }
081        
082        /**
083         * Returns the angle of the Cartesian point (x,y)
084         * @param x x-component
085         * @param y y-component
086         * @return the angle
087         */
088        public static double angle(double x, double y) {
089                return atan2(y, x);
090        }
091        
092        /**
093         * Returns the polar coordinates of the Cartesian point (x,y).
094         * @param x x-component
095         * @param y y-component
096         * @return a 2-element array holding the polar coordinates (radius, angle)
097         */
098        public static double[] toPolar(double x, double y) {
099                return new double[] {hypot(x, y), atan2(y, x)};
100        }
101        
102        /**
103         * Returns the polar coordinates of the Cartesian point xy.
104         * @param xy the Cartesian coordinates
105         * @return a 2-element array holding the polar coordinates (radius, angle)
106         */
107        public static double[] toPolar(double[] xy) {
108                return toPolar(xy[0], xy[1]);
109        }
110        
111        /**
112         * Converts polar point coordinates to Cartesian coordinates.
113         * @param radius the radius 
114         * @param angle the angle
115         * @return a 2-element array holding the Cartesian coordinates (x,y)
116         */
117        public static double[] toCartesian(double radius, double angle) {
118                return new double[] {radius * cos(angle), radius * sin(angle)};
119        }
120        
121        /**
122         * Converts polar point coordinates to Cartesian coordinates.
123         * @param polar the polar coordinates (radius, angle)
124         * @return a 2-element array holding the Cartesian coordinates (x,y)
125         */
126        public static double[] toCartesian(double[] polar) {
127                return toCartesian(polar[0], polar[1]);
128        }
129
130        /**
131         * Integer version of the modulus operator ({@code a mod b}). Also see <a
132         * href="http://en.wikipedia.org/wiki/Modulo_operation">here</a>. Calls {@code Math.floorMod(a,b)} (available in
133         * Java 8 and higher).
134         *
135         * @param a dividend
136         * @param b divisor (modulus), must be nonzero
137         * @return {@code a mod b}
138         */
139        public static int mod(int a, int b) {
140                return floorMod(a, b);
141        }
142
143        /**
144         * Non-integer version of modulus operator using floored division (see <a
145         * href="http://en.wikipedia.org/wiki/Modulo_operation">here</a>), with results identical to Mathematica. Calculates
146         * {@code a mod b} for floating-point arguments. An exception is thrown if  {@code b} is zero. Examples:
147         * <pre>
148         * mod( 3.5, 2.1) =  1.4
149         * mod(-3.5, 2.1) =  0.7
150         * mod( 3.5,-2.1) = -0.7
151         * mod(-3.5,-2.1) = -1.4</pre>
152         *
153         * @param a dividend
154         * @param b divisor (modulus), must be nonzero
155         * @return {@code a mod b}
156         */
157        public static double mod(double a, double b) {
158                if (isZero(b))
159                                throw new IllegalArgumentException("zero modulus in mod");
160                return a - b * floor(a / b);
161        }
162
163        /**
164         * Test for zero (float version) using a predefined tolerance. Returns true if the argument's absolute value is less
165         * than {@link #EPSILON_FLOAT}.
166         *
167         * @param x quantity to be tested
168         * @return true if argument is close to zero
169         */
170        public static boolean isZero(float x) {
171                return abs(x) < EPSILON_FLOAT;
172        }
173
174        /**
175         * Test for zero (float version) using a specified tolerance. Returns true if the argument's absolute value is less
176         * than the specified tolerance.
177         *
178         * @param x quantity to be tested
179         * @param tolerance the tolerance to be used
180         * @return true if argument is close to zero
181         */
182        public static boolean isZero(float x, float tolerance) {
183                return abs(x) < tolerance;
184        }
185
186        /**
187         * Test for zero (double version) using a predefined tolerance. Returns true if the argument's absolute value is
188         * less than {@link #EPSILON_DOUBLE}.
189         *
190         * @param x quantity to be tested
191         * @return true if argument is close to zero
192         */
193        public static boolean isZero(double x) {
194                return abs(x) < EPSILON_DOUBLE;
195        }
196
197        /**
198         * Test for zero (double version) using a specified tolerance. Returns true if the argument's absolute value is less
199         * than the specified tolerance.
200         *
201         * @param x quantity to be tested
202         * @param tolerance the tolerance to be used
203         * @return true if argument is close to zero
204         */
205        public static boolean isZero(double x, double tolerance) {
206                return abs(x) < tolerance;
207        }
208
209        /**
210         * Test for numerical equality (double version) using the default tolerance. Returns true if the absolute difference
211         * of the arguments is less than {@link #EPSILON_DOUBLE}.
212         *
213         * @param x first argument
214         * @param y second argument
215         * @return true if the absolute difference of the arguments is less than the tolerance
216         */
217        public static boolean equals(double x, double y) {
218                return isZero(x - y);
219        }
220
221        /**
222         * Test for numerical equality (double version) using a specific tolerance.
223         *
224         * @param x first argument
225         * @param y second argument
226         * @param tolerance the maximum (absolute) deviation
227         * @return true if the absolute difference of the arguments is less than the tolerance
228         */
229        public static boolean equals(double x, double y, double tolerance) {
230                return isZero(x - y, tolerance);
231        }
232
233        /**
234         * Test for numerical equality (float version) using the default tolerance. Returns true if the absolute difference
235         * of the arguments is less than {@link #EPSILON_FLOAT}.
236         *
237         * @param x first argument
238         * @param y second argument
239         * @return true if the absolute difference of the arguments is less than the tolerance
240         */
241        public static boolean equals(float x, float y) {
242                return isZero(x - y);
243        }
244
245        /**
246         * Test for numerical equality (float version) using a specific tolerance.
247         *
248         * @param x first argument
249         * @param y second argument
250         * @param tolerance the maximum (absolute) deviation
251         * @return true if the absolute difference of the arguments is less than the tolerance
252         */
253        public static boolean equals(float x, float y, float tolerance) {
254                return isZero(x - y, tolerance);
255        }
256        
257        // ---------------------------------------------------------------------------
258        
259        /**
260         * Returns the maximum of one or more integer values.
261         * @param a the first value
262         * @param vals more values
263         * @return the maximum value
264         */
265        public static int max(int a, int... vals) {
266                int maxVal = a;
267                for (int v : vals) {
268                        if (v > maxVal) {
269                                maxVal = v;
270                        }
271                }
272                return maxVal;
273        }
274        
275        /**
276         * Returns the minimum of one or more integer values.
277         * @param a the first value
278         * @param vals more values
279         * @return the minimum value
280         */
281        public static int min(int a, int... vals) {
282                int minVal = a;
283                for (int v : vals) {
284                        if (v < minVal) {
285                                minVal = v;
286                        }
287                }
288                return minVal;
289        }
290        
291        /**
292         * Returns the maximum of one or more float values.
293         * @param a the first value
294         * @param vals more values
295         * @return the maximum value
296         */
297        public static float max(float a, float... vals) {
298                float maxVal = a;
299                for (float v : vals) {
300                        if (v > maxVal) {
301                                maxVal = v;
302                        }
303                }
304                return maxVal;
305        }
306        
307        /**
308         * Returns the minimum of one or more float values.
309         * @param a the first value
310         * @param vals more values
311         * @return the minimum value
312         */
313        public static float min(float a, float... vals) {
314                float minVal = a;
315                for (float v : vals) {
316                        if (v < minVal) {
317                                minVal = v;
318                        }
319                }
320                return minVal;
321        }
322        
323        /**
324         * Returns the maximum of one or more double values.
325         * @param a the first value
326         * @param vals more values
327         * @return the maximum value
328         */
329        public static double max(double a, double... vals) {
330                double maxVal = a;
331                for (double v : vals) {
332                        if (v > maxVal) {
333                                maxVal = v;
334                        }
335                }
336                return maxVal;
337        }
338        
339        /**
340         * Returns the minimum of one or more double values.
341         * @param a the first value
342         * @param vals more values
343         * @return the minimum value
344         */
345        public static double min(double a, double... vals) {
346                double minVal = a;
347                for (double v : vals) {
348                        if (v < minVal) {
349                                minVal = v;
350                        }
351                }
352                return minVal;
353        }
354
355        /**
356         * Limits the first argument to the clipping interval specified by the two other arguments ({@code double} version).
357         * The clipped value is returned. Throws an exception if the clipping interval is empty.
358         *
359         * @param x the value to be clipped
360         * @param low the lower boundary of the clipping interval
361         * @param high the upper boundary of the clipping interval
362         * @return the clipped value
363         */
364        public static double clipTo(double x, double low, double high) {
365                if (low > high) {
366                        throw new IllegalArgumentException("clip interval low > high");
367                }
368                if (x < low) return low;
369                if (x > high) return high;
370                return x;
371        }
372
373        /**
374         * Limits the first argument to the clipping interval specified by the two other arguments ({@code float} version).
375         * The clipped value is returned. Throws an exception if the clipping interval is empty.
376         *
377         * @param x the value to be clipped
378         * @param low the lower boundary of the clipping interval
379         * @param high the upper boundary of the clipping interval
380         * @return the clipped value
381         */
382        public static float clipTo(float x, float low, float high) {
383                if (low > high) {
384                        throw new IllegalArgumentException("clip interval low > high");
385                }
386                if (x < low) return low;
387                if (x > high) return high;
388                return x;
389        }
390
391        /**
392         * Limits the first argument to the clipping interval specified by the two other arguments ({@code float} version).
393         * The clipped value is returned. Throws an exception if the clipping interval is empty.
394         *
395         * @param x the value to be clipped
396         * @param low the lower boundary of the clipping interval
397         * @param high the upper boundary of the clipping interval
398         * @return the clipped value
399         */
400        public static int clipTo(int x, int low, int high) {
401                if (low > high) {
402                        throw new IllegalArgumentException("clip interval low > high");
403                }
404                if (x < low) return low;
405                if (x > high) return high;
406                return x;
407        }
408
409        /**
410         * Returns the two real roots of the quadratic function f(x) = ax^2 + bx + c. Null is returned if roots are
411         * non-real.
412         *
413         * @param a function coefficient
414         * @param b function coefficient
415         * @param c function coefficient
416         * @return an array with the two roots x1, x2 (in no particular order) or null if the function has complex roots
417         */
418        public static double[] getRealRoots(double a, double b, double c) {
419                double d = sqr(b) - 4 * a * c;
420                if (d < 0) {
421                        return null;    // complex roots
422                }
423                double dr = Math.sqrt(d);       
424                double x1 = (-b - dr) / (2 * a);
425                double x2 = (-b + dr) / (2 * a);
426                return new double[] {x1, x2};
427        }
428
429}