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