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