001/*
002Copyright 2006 Jerry Huxtable
003
004Licensed under the Apache License, Version 2.0 (the "License");
005you may not use this file except in compliance with the License.
006You may obtain a copy of the License at
007
008   http://www.apache.org/licenses/LICENSE-2.0
009
010Unless required by applicable law or agreed to in writing, software
011distributed under the License is distributed on an "AS IS" BASIS,
012WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
013See the License for the specific language governing permissions and
014limitations under the License.
015*/
016
017package com.jhlabs.image;
018
019/**
020 * A class containing static math methods useful for image processing.
021 */
022public class ImageMath {
023
024    /**
025     * The value of pi as a float.
026     */
027        public final static float PI = (float)Math.PI;
028
029    /**
030     * The value of half pi as a float.
031     */
032        public final static float HALF_PI = (float)Math.PI/2.0f;
033
034    /**
035     * The value of quarter pi as a float.
036     */
037        public final static float QUARTER_PI = (float)Math.PI/4.0f;
038
039    /**
040     * The value of two pi as a float.
041     */
042        public final static float TWO_PI = (float)Math.PI*2.0f;
043
044        /**
045         * Apply a bias to a number in the unit interval, moving numbers towards 0 or 1
046         * according to the bias parameter.
047         * @param a the number to bias
048         * @param b the bias parameter. 0.5 means no change, smaller values bias towards 0, larger towards 1.
049         * @return the output value
050         */
051        public static float bias(float a, float b) {
052//              return (float)Math.pow(a, Math.log(b) / Math.log(0.5));
053                return a/((1.0f/b-2)*(1.0f-a)+1);
054        }
055
056        /**
057         * A variant of the gamma function.
058         * @param a the number to apply gain to
059         * @param b the gain parameter. 0.5 means no change, smaller values reduce gain, larger values increase gain.
060         * @return the output value
061         */
062        public static float gain(float a, float b) {
063/*
064                float p = (float)Math.log(1.0 - b) / (float)Math.log(0.5);
065
066                if (a < .001)
067                        return 0.0f;
068                else if (a > .999)
069                        return 1.0f;
070                if (a < 0.5)
071                        return (float)Math.pow(2 * a, p) / 2;
072                else
073                        return 1.0f - (float)Math.pow(2 * (1. - a), p) / 2;
074*/
075                float c = (1.0f/b-2.0f) * (1.0f-2.0f*a);
076                if (a < 0.5)
077                        return a/(c+1.0f);
078                else
079                        return (c-a)/(c-1.0f);
080        }
081
082        /**
083         * The step function. Returns 0 below a threshold, 1 above.
084         * @param a the threshold position
085         * @param x the input parameter
086         * @return the output value - 0 or 1
087         */
088        public static float step(float a, float x) {
089                return (x < a) ? 0.0f : 1.0f;
090        }
091
092        /**
093         * The pulse function. Returns 1 between two thresholds, 0 outside.
094         * @param a the lower threshold position
095         * @param b the upper threshold position
096         * @param x the input parameter
097         * @return the output value - 0 or 1
098         */
099        public static float pulse(float a, float b, float x) {
100                return (x < a || x >= b) ? 0.0f : 1.0f;
101        }
102
103        /**
104         * A smoothed pulse function. A cubic function is used to smooth the step between two thresholds.
105         * @param a1 the lower threshold position for the start of the pulse
106         * @param a2 the upper threshold position for the start of the pulse
107         * @param b1 the lower threshold position for the end of the pulse
108         * @param b2 the upper threshold position for the end of the pulse
109         * @param x the input parameter
110         * @return the output value
111         */
112        public static float smoothPulse(float a1, float a2, float b1, float b2, float x) {
113                if (x < a1 || x >= b2)
114                        return 0;
115                if (x >= a2) {
116                        if (x < b1)
117                                return 1.0f;
118                        x = (x - b1) / (b2 - b1);
119                        return 1.0f - (x*x * (3.0f - 2.0f*x));
120                }
121                x = (x - a1) / (a2 - a1);
122                return x*x * (3.0f - 2.0f*x);
123        }
124
125        /**
126         * A smoothed step function. A cubic function is used to smooth the step between two thresholds.
127         * @param a the lower threshold position
128         * @param b the upper threshold position
129         * @param x the input parameter
130         * @return the output value
131         */
132        public static float smoothStep(float a, float b, float x) {
133                if (x < a)
134                        return 0;
135                if (x >= b)
136                        return 1;
137                x = (x - a) / (b - a);
138                return x*x * (3 - 2*x);
139        }
140
141        /**
142         * A "circle up" function. Returns y on a unit circle given 1-x. Useful for forming bevels.
143         * @param x the input parameter in the range 0..1
144         * @return the output value
145         */
146        public static float circleUp(float x) {
147                x = 1-x;
148                return (float)Math.sqrt(1-x*x);
149        }
150
151        /**
152         * A "circle down" function. Returns 1-y on a unit circle given x. Useful for forming bevels.
153         * @param x the input parameter in the range 0..1
154         * @return the output value
155         */
156        public static float circleDown(float x) {
157                return 1.0f-(float)Math.sqrt(1-x*x);
158        }
159
160        /**
161         * Clamp a value to an interval.
162         * @param a the lower clamp threshold
163         * @param b the upper clamp threshold
164         * @param x the input parameter
165         * @return the clamped value
166         */
167        public static float clamp(float x, float a, float b) {
168                return (x < a) ? a : (x > b) ? b : x;
169        }
170
171        /**
172         * Clamp a value to an interval.
173         * @param a the lower clamp threshold
174         * @param b the upper clamp threshold
175         * @param x the input parameter
176         * @return the clamped value
177         */
178        public static int clamp(int x, int a, int b) {
179                return (x < a) ? a : (x > b) ? b : x;
180        }
181
182        /**
183         * Return a mod b. This differs from the % operator with respect to negative numbers.
184         * @param a the dividend
185         * @param b the divisor
186         * @return a mod b
187         */
188        public static double mod(double a, double b) {
189                int n = (int)(a/b);
190                
191                a -= n*b;
192                if (a < 0)
193                        return a + b;
194                return a;
195        }
196
197        /**
198         * Return a mod b. This differs from the % operator with respect to negative numbers.
199         * @param a the dividend
200         * @param b the divisor
201         * @return a mod b
202         */
203        public static float mod(float a, float b) {
204                int n = (int)(a/b);
205                
206                a -= n*b;
207                if (a < 0)
208                        return a + b;
209                return a;
210        }
211
212        /**
213         * Return a mod b. This differs from the % operator with respect to negative numbers.
214         * @param a the dividend
215         * @param b the divisor
216         * @return a mod b
217         */
218        public static int mod(int a, int b) {
219                int n = a/b;
220                
221                a -= n*b;
222                if (a < 0)
223                        return a + b;
224                return a;
225        }
226
227        /**
228         * The triangle function. Returns a repeating triangle shape in the range 0..1 with wavelength 1.0
229         * @param x the input parameter
230         * @return the output value
231         */
232        public static float triangle(float x) {
233                float r = mod(x, 1.0f);
234                return 2.0f*(r < 0.5 ? r : 1-r);
235        }
236
237        /**
238         * Linear interpolation.
239         * @param t the interpolation parameter
240         * @param a the lower interpolation range
241         * @param b the upper interpolation range
242         * @return the interpolated value
243         */
244        public static float lerp(float t, float a, float b) {
245                return a + t * (b - a);
246        }
247        
248        /**
249         * Linear interpolation.
250         * @param t the interpolation parameter
251         * @param a the lower interpolation range
252         * @param b the upper interpolation range
253         * @return the interpolated value
254         */
255        public static int lerp(float t, int a, int b) {
256                return (int)(a + t * (b - a));
257        }
258
259        /**
260         * Linear interpolation of ARGB values.
261         * @param t the interpolation parameter
262         * @param rgb1 the lower interpolation range
263         * @param rgb2 the upper interpolation range
264         * @return the interpolated value
265         */
266        public static int mixColors(float t, int rgb1, int rgb2) {
267                int a1 = (rgb1 >> 24) & 0xff;
268                int r1 = (rgb1 >> 16) & 0xff;
269                int g1 = (rgb1 >> 8) & 0xff;
270                int b1 = rgb1 & 0xff;
271                int a2 = (rgb2 >> 24) & 0xff;
272                int r2 = (rgb2 >> 16) & 0xff;
273                int g2 = (rgb2 >> 8) & 0xff;
274                int b2 = rgb2 & 0xff;
275                a1 = lerp(t, a1, a2);
276                r1 = lerp(t, r1, r2);
277                g1 = lerp(t, g1, g2);
278                b1 = lerp(t, b1, b2);
279                return (a1 << 24) | (r1 << 16) | (g1 << 8) | b1;
280        }
281
282        /**
283         * Bilinear interpolation of ARGB values.
284         * @param x the X interpolation parameter 0..1
285         * @param y the y interpolation parameter 0..1
286         * @param rgb array of four ARGB values in the order NW, NE, SW, SE
287         * @return the interpolated value
288         */
289        public static int bilinearInterpolate(float x, float y, int nw, int ne, int sw, int se) {
290                float m0, m1;
291                int a0 = (nw >> 24) & 0xff;
292                int r0 = (nw >> 16) & 0xff;
293                int g0 = (nw >> 8) & 0xff;
294                int b0 = nw & 0xff;
295                int a1 = (ne >> 24) & 0xff;
296                int r1 = (ne >> 16) & 0xff;
297                int g1 = (ne >> 8) & 0xff;
298                int b1 = ne & 0xff;
299                int a2 = (sw >> 24) & 0xff;
300                int r2 = (sw >> 16) & 0xff;
301                int g2 = (sw >> 8) & 0xff;
302                int b2 = sw & 0xff;
303                int a3 = (se >> 24) & 0xff;
304                int r3 = (se >> 16) & 0xff;
305                int g3 = (se >> 8) & 0xff;
306                int b3 = se & 0xff;
307
308                float cx = 1.0f-x;
309                float cy = 1.0f-y;
310
311                m0 = cx * a0 + x * a1;
312                m1 = cx * a2 + x * a3;
313                int a = (int)(cy * m0 + y * m1);
314
315                m0 = cx * r0 + x * r1;
316                m1 = cx * r2 + x * r3;
317                int r = (int)(cy * m0 + y * m1);
318
319                m0 = cx * g0 + x * g1;
320                m1 = cx * g2 + x * g3;
321                int g = (int)(cy * m0 + y * m1);
322
323                m0 = cx * b0 + x * b1;
324                m1 = cx * b2 + x * b3;
325                int b = (int)(cy * m0 + y * m1);
326
327                return (a << 24) | (r << 16) | (g << 8) | b;
328        }
329
330        /**
331         * Return the NTSC gray level of an RGB value.
332         * @param rgb1 the input pixel
333         * @return the gray level (0-255)
334         */
335        public static int brightnessNTSC(int rgb) {
336                int r = (rgb >> 16) & 0xff;
337                int g = (rgb >> 8) & 0xff;
338                int b = rgb & 0xff;
339                return (int)(r*0.299f + g*0.587f + b*0.114f);
340        }
341        
342        // Catmull-Rom splines
343        private final static float m00 = -0.5f;
344        private final static float m01 =  1.5f;
345        private final static float m02 = -1.5f;
346        private final static float m03 =  0.5f;
347        private final static float m10 =  1.0f;
348        private final static float m11 = -2.5f;
349        private final static float m12 =  2.0f;
350        private final static float m13 = -0.5f;
351        private final static float m20 = -0.5f;
352        private final static float m21 =  0.0f;
353        private final static float m22 =  0.5f;
354        private final static float m23 =  0.0f;
355        private final static float m30 =  0.0f;
356        private final static float m31 =  1.0f;
357        private final static float m32 =  0.0f;
358        private final static float m33 =  0.0f;
359
360        /**
361         * Compute a Catmull-Rom spline.
362         * @param x the input parameter
363         * @param numKnots the number of knots in the spline
364         * @param knots the array of knots
365         * @return the spline value
366         */
367        public static float spline(float x, int numKnots, float[] knots) {
368                int span;
369                int numSpans = numKnots - 3;
370                float k0, k1, k2, k3;
371                float c0, c1, c2, c3;
372                
373                if (numSpans < 1)
374                        throw new IllegalArgumentException("Too few knots in spline");
375                
376                x = clamp(x, 0, 1) * numSpans;
377                span = (int)x;
378                if (span > numKnots-4)
379                        span = numKnots-4;
380                x -= span;
381
382                k0 = knots[span];
383                k1 = knots[span+1];
384                k2 = knots[span+2];
385                k3 = knots[span+3];
386                
387                c3 = m00*k0 + m01*k1 + m02*k2 + m03*k3;
388                c2 = m10*k0 + m11*k1 + m12*k2 + m13*k3;
389                c1 = m20*k0 + m21*k1 + m22*k2 + m23*k3;
390                c0 = m30*k0 + m31*k1 + m32*k2 + m33*k3;
391                
392                return ((c3*x + c2)*x + c1)*x + c0;
393        }
394        
395        /**
396         * Compute a Catmull-Rom spline, but with variable knot spacing.
397         * @param x the input parameter
398         * @param numKnots the number of knots in the spline
399         * @param xknots the array of knot x values
400         * @param yknots the array of knot y values
401         * @return the spline value
402         */
403        public static float spline(float x, int numKnots, int[] xknots, int[] yknots) {
404                int span;
405                int numSpans = numKnots - 3;
406                float k0, k1, k2, k3;
407                float c0, c1, c2, c3;
408                
409                if (numSpans < 1)
410                        throw new IllegalArgumentException("Too few knots in spline");
411                
412                for (span = 0; span < numSpans; span++)
413                        if (xknots[span+1] > x)
414                                break;
415                if (span > numKnots-3)
416                        span = numKnots-3;
417                float t = (float)(x-xknots[span]) / (xknots[span+1]-xknots[span]);
418                span--;
419                if (span < 0) {
420                        span = 0;
421                        t = 0;
422                }
423
424                k0 = yknots[span];
425                k1 = yknots[span+1];
426                k2 = yknots[span+2];
427                k3 = yknots[span+3];
428                
429                c3 = m00*k0 + m01*k1 + m02*k2 + m03*k3;
430                c2 = m10*k0 + m11*k1 + m12*k2 + m13*k3;
431                c1 = m20*k0 + m21*k1 + m22*k2 + m23*k3;
432                c0 = m30*k0 + m31*k1 + m32*k2 + m33*k3;
433                
434                return ((c3*t + c2)*t + c1)*t + c0;
435        }
436
437        /**
438         * Compute a Catmull-Rom spline for RGB values.
439         * @param x the input parameter
440         * @param numKnots the number of knots in the spline
441         * @param knots the array of knots
442         * @return the spline value
443         */
444        public static int colorSpline(float x, int numKnots, int[] knots) {
445                int span;
446                int numSpans = numKnots - 3;
447                float k0, k1, k2, k3;
448                float c0, c1, c2, c3;
449                
450                if (numSpans < 1)
451                        throw new IllegalArgumentException("Too few knots in spline");
452                
453                x = clamp(x, 0, 1) * numSpans;
454                span = (int)x;
455                if (span > numKnots-4)
456                        span = numKnots-4;
457                x -= span;
458
459                int v = 0;
460                for (int i = 0; i < 4; i++) {
461                        int shift = i * 8;
462                        
463                        k0 = (knots[span] >> shift) & 0xff;
464                        k1 = (knots[span+1] >> shift) & 0xff;
465                        k2 = (knots[span+2] >> shift) & 0xff;
466                        k3 = (knots[span+3] >> shift) & 0xff;
467                        
468                        c3 = m00*k0 + m01*k1 + m02*k2 + m03*k3;
469                        c2 = m10*k0 + m11*k1 + m12*k2 + m13*k3;
470                        c1 = m20*k0 + m21*k1 + m22*k2 + m23*k3;
471                        c0 = m30*k0 + m31*k1 + m32*k2 + m33*k3;
472                        int n = (int)(((c3*x + c2)*x + c1)*x + c0);
473                        if (n < 0)
474                                n = 0;
475                        else if (n > 255)
476                                n = 255;
477                        v |= n << shift;
478                }
479                
480                return v;
481        }
482
483        /**
484         * Compute a Catmull-Rom spline for RGB values, but with variable knot spacing.
485         * @param x the input parameter
486         * @param numKnots the number of knots in the spline
487         * @param xknots the array of knot x values
488         * @param yknots the array of knot y values
489         * @return the spline value
490         */
491        public static int colorSpline(int x, int numKnots, int[] xknots, int[] yknots) {
492                int span;
493                int numSpans = numKnots - 3;
494                float k0, k1, k2, k3;
495                float c0, c1, c2, c3;
496                
497                if (numSpans < 1)
498                        throw new IllegalArgumentException("Too few knots in spline");
499                
500                for (span = 0; span < numSpans; span++)
501                        if (xknots[span+1] > x)
502                                break;
503                if (span > numKnots-3)
504                        span = numKnots-3;
505                float t = (float)(x-xknots[span]) / (xknots[span+1]-xknots[span]);
506                span--;
507                if (span < 0) {
508                        span = 0;
509                        t = 0;
510                }
511
512                int v = 0;
513                for (int i = 0; i < 4; i++) {
514                        int shift = i * 8;
515                        
516                        k0 = (yknots[span] >> shift) & 0xff;
517                        k1 = (yknots[span+1] >> shift) & 0xff;
518                        k2 = (yknots[span+2] >> shift) & 0xff;
519                        k3 = (yknots[span+3] >> shift) & 0xff;
520                        
521                        c3 = m00*k0 + m01*k1 + m02*k2 + m03*k3;
522                        c2 = m10*k0 + m11*k1 + m12*k2 + m13*k3;
523                        c1 = m20*k0 + m21*k1 + m22*k2 + m23*k3;
524                        c0 = m30*k0 + m31*k1 + m32*k2 + m33*k3;
525                        int n = (int)(((c3*t + c2)*t + c1)*t + c0);
526                        if (n < 0)
527                                n = 0;
528                        else if (n > 255)
529                                n = 255;
530                        v |= n << shift;
531                }
532                
533                return v;
534        }
535
536        /**
537         * An implementation of Fant's resampling algorithm.
538         * @param source the source pixels
539         * @param dest the destination pixels
540         * @param length the length of the scanline to resample
541         * @param offset the start offset into the arrays
542         * @param stride the offset between pixels in consecutive rows
543         * @param out an array of output positions for each pixel
544         */
545        public static void resample(int[] source, int[] dest, int length, int offset, int stride, float[] out) {
546                int i, j;
547                float sizfac;
548                float inSegment;
549                float outSegment;
550                int a, r, g, b, nextA, nextR, nextG, nextB;
551                float aSum, rSum, gSum, bSum;
552                float[] in;
553                int srcIndex = offset;
554                int destIndex = offset;
555                int lastIndex = source.length;
556                int rgb;
557
558                in = new float[length+2];
559                i = 0;
560                for (j = 0; j < length; j++) {
561                        while (out[i+1] < j)
562                                i++;
563                        in[j] = i + (float) (j - out[i]) / (out[i + 1] - out[i]);
564//                      in[j] = ImageMath.clamp( in[j], 0, length-1 );
565                }
566                in[length] = length;
567                in[length+1] = length;
568
569                inSegment  = 1.0f;
570                outSegment = in[1];
571                sizfac = outSegment;
572                aSum = rSum = gSum = bSum = 0.0f;
573                rgb = source[srcIndex];
574                a = (rgb >> 24) & 0xff;
575                r = (rgb >> 16) & 0xff;
576                g = (rgb >> 8) & 0xff;
577                b = rgb & 0xff;
578                srcIndex += stride;
579                rgb = source[srcIndex];
580                nextA = (rgb >> 24) & 0xff;
581                nextR = (rgb >> 16) & 0xff;
582                nextG = (rgb >> 8) & 0xff;
583                nextB = rgb & 0xff;
584                srcIndex += stride;
585                i = 1;
586
587                while (i <= length) {
588                        float aIntensity = inSegment * a + (1.0f - inSegment) * nextA;
589                        float rIntensity = inSegment * r + (1.0f - inSegment) * nextR;
590                        float gIntensity = inSegment * g + (1.0f - inSegment) * nextG;
591                        float bIntensity = inSegment * b + (1.0f - inSegment) * nextB;
592                        if (inSegment < outSegment) {
593                                aSum += (aIntensity * inSegment);
594                                rSum += (rIntensity * inSegment);
595                                gSum += (gIntensity * inSegment);
596                                bSum += (bIntensity * inSegment);
597                                outSegment -= inSegment;
598                                inSegment = 1.0f;
599                                a = nextA;
600                                r = nextR;
601                                g = nextG;
602                                b = nextB;
603                                if (srcIndex < lastIndex)
604                                        rgb = source[srcIndex];
605                                nextA = (rgb >> 24) & 0xff;
606                                nextR = (rgb >> 16) & 0xff;
607                                nextG = (rgb >> 8) & 0xff;
608                                nextB = rgb & 0xff;
609                                srcIndex += stride;
610                        } else {
611                                aSum += (aIntensity * outSegment);
612                                rSum += (rIntensity * outSegment);
613                                gSum += (gIntensity * outSegment);
614                                bSum += (bIntensity * outSegment);
615                                dest[destIndex] = 
616                                        ((int)Math.min(aSum/sizfac, 255) << 24) |
617                                        ((int)Math.min(rSum/sizfac, 255) << 16) |
618                                        ((int)Math.min(gSum/sizfac, 255) << 8) |
619                                        (int)Math.min(bSum/sizfac, 255);
620                                destIndex += stride;
621                                aSum = rSum = gSum = bSum = 0.0f;
622                                inSegment -= outSegment;
623                                outSegment = in[i+1] - in[i];
624                                sizfac = outSegment;
625                                i++;
626                        }
627                }
628        }
629
630    /**
631         * Premultiply a block of pixels
632         */
633        public static void premultiply( int[] p, int offset, int length ) {
634        length += offset;
635                for ( int i = offset; i < length; i ++ ) {
636            int rgb = p[i];
637            int a = (rgb >> 24) & 0xff;
638            int r = (rgb >> 16) & 0xff;
639            int g = (rgb >> 8) & 0xff;
640            int b = rgb & 0xff;
641            float f = a * (1.0f / 255.0f);
642            r *= f;
643            g *= f;
644            b *= f;
645            p[i] = (a << 24) | (r << 16) | (g << 8) | b;
646        }
647    }
648
649    /**
650         * Premultiply a block of pixels
651         */
652        public static void unpremultiply( int[] p, int offset, int length ) {
653        length += offset;
654                for ( int i = offset; i < length; i ++ ) {
655            int rgb = p[i];
656            int a = (rgb >> 24) & 0xff;
657            int r = (rgb >> 16) & 0xff;
658            int g = (rgb >> 8) & 0xff;
659            int b = rgb & 0xff;
660            if ( a != 0 && a != 255 ) {
661                float f = 255.0f / a;
662                r *= f;
663                g *= f;
664                b *= f;
665                if ( r > 255 )
666                    r = 255;
667                if ( g > 255 )
668                    g = 255;
669                if ( b > 255 )
670                    b = 255;
671                p[i] = (a << 24) | (r << 16) | (g << 8) | b;
672            }
673        }
674    }
675}