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 **/
019package lucee.runtime.img.math;
020
021import java.util.Random;
022
023/**
024 * Sparse Convolution Noise. This is computationally very expensive, but worth it.
025 */
026public class SCNoise implements Function1D, Function2D, Function3D {
027
028        private static Random randomGenerator = new Random();
029        
030        public float evaluate(float x) {
031                return evaluate(x, .1f);
032        }
033        
034        public float evaluate(float x, float y) {
035                int i, j,  h, n;
036                int ix, iy;
037                float sum = 0;
038                float fx, fy, dx, dy, distsq;
039
040                if (impulseTab == null)
041                        impulseTab = impulseTabInit(665);
042
043                ix = floor(x); fx = x - ix;
044                iy = floor(y); fy = y - iy;
045                
046                /* Perform the sparse convolution. */
047                int m = 2;
048                for (i = -m; i <= m; i++) {
049                  for (j = -m; j <= m; j++) {
050                        /* Compute voxel hash code. */
051                        h = perm[(ix+i + perm[(iy+j)&TABMASK])&TABMASK];
052                        
053                        for (n = NIMPULSES; n > 0; n--, h = (h+1) & TABMASK) {
054                                /* Convolve filter and impulse. */
055                                int h4 = h*4;
056                                dx = fx - (i + impulseTab[h4++]);
057                                dy = fy - (j + impulseTab[h4++]);
058                                distsq = dx*dx + dy*dy;
059                                sum += catrom2(distsq) * impulseTab[h4];
060                        }
061                  }
062                }
063
064                return sum / NIMPULSES;
065        }
066        
067        public float evaluate(float x, float y, float z) {
068                int i, j, k, h, n;
069                int ix, iy, iz;
070                float sum = 0;
071                float fx, fy, fz, dx, dy, dz, distsq;
072
073                if (impulseTab == null)
074                        impulseTab = impulseTabInit(665);
075
076                ix = floor(x); fx = x - ix;
077                iy = floor(y); fy = y - iy;
078                iz = floor(z); fz = z - iz;
079                
080                /* Perform the sparse convolution. */
081                int m = 2;
082                for (i = -m; i <= m; i++) {
083                  for (j = -m; j <= m; j++) {
084                        for (k = -m; k <= m; k++) {
085                                /* Compute voxel hash code. */
086                                h = perm[(ix+i + perm[(iy+j + perm[(iz+k)&TABMASK])&TABMASK])&TABMASK];
087                                
088                                for (n = NIMPULSES; n > 0; n--, h = (h+1) & TABMASK) {
089                                        /* Convolve filter and impulse. */
090                                        int h4 = h*4;
091                                        dx = fx - (i + impulseTab[h4++]);
092                                        dy = fy - (j + impulseTab[h4++]);
093                                        dz = fz - (k + impulseTab[h4++]);
094                                        distsq = dx*dx + dy*dy + dz*dz;
095                                        sum += catrom2(distsq) * impulseTab[h4];
096                                }
097                        }
098                  }
099                }
100
101                return sum / NIMPULSES;
102        }
103        
104        public short[] perm = {
105                        225,155,210,108,175,199,221,144,203,116, 70,213, 69,158, 33,252,
106                          5, 82,173,133,222,139,174, 27,  9, 71, 90,246, 75,130, 91,191,
107                        169,138,  2,151,194,235, 81,  7, 25,113,228,159,205,253,134,142,
108                        248, 65,224,217, 22,121,229, 63, 89,103, 96,104,156, 17,201,129,
109                         36,  8,165,110,237,117,231, 56,132,211,152, 20,181,111,239,218,
110                        170,163, 51,172,157, 47, 80,212,176,250, 87, 49, 99,242,136,189,
111                        162,115, 44, 43,124, 94,150, 16,141,247, 32, 10,198,223,255, 72,
112                         53,131, 84, 57,220,197, 58, 50,208, 11,241, 28,  3,192, 62,202,
113                         18,215,153, 24, 76, 41, 15,179, 39, 46, 55,  6,128,167, 23,188,
114                        106, 34,187,140,164, 73,112,182,244,195,227, 13, 35, 77,196,185,
115                         26,200,226,119, 31,123,168,125,249, 68,183,230,177,135,160,180,
116                         12,  1,243,148,102,166, 38,238,251, 37,240,126, 64, 74,161, 40,
117                        184,149,171,178,101, 66, 29, 59,146, 61,254,107, 42, 86,154,  4,
118                        236,232,120, 21,233,209, 45, 98,193,114, 78, 19,206, 14,118,127,
119                         48, 79,147, 85, 30,207,219, 54, 88,234,190,122, 95, 67,143,109,
120                        137,214,145, 93, 92,100,245,  0,216,186, 60, 83,105, 97,204, 52
121        };
122
123        private final static int TABSIZE = 256;
124        private final static int TABMASK = (TABSIZE-1);
125        private final static int NIMPULSES = 3;
126
127        private static float[] impulseTab;
128
129        public static int floor(float x) {
130                int ix = (int)x;
131                if (x < 0 && x != ix)
132                        return ix-1;
133                return ix;
134        }
135
136        private final static int SAMPRATE = 100;  /* table entries per unit distance */
137        private final static int NENTRIES = (4*SAMPRATE+1);
138        private static float[] table;
139
140        public float catrom2(float d) {
141                float x;
142                int i;
143
144                if (d >= 4)
145                        return 0;
146
147                if (table == null) {
148                        table = new float[NENTRIES];
149                        for (i = 0; i < NENTRIES; i++) {
150                                x = i/(float)SAMPRATE;
151                                x = (float)Math.sqrt(x);
152                                if (x < 1)
153                                        table[i] = 0.5f * (2+x*x*(-5+x*3));
154                                else
155                                        table[i] = 0.5f * (4+x*(-8+x*(5-x)));
156                        }
157                }
158
159                d = d*SAMPRATE + 0.5f;
160                i = floor(d);
161                if (i >= NENTRIES)
162                        return 0;
163                return table[i];
164        }
165
166        static float[] impulseTabInit(int seed) {
167                float[] impulseTab = new float[TABSIZE*4];
168
169                randomGenerator = new Random(seed); /* Set random number generator seed. */
170                for (int i = 0; i < TABSIZE; i++) {
171                        impulseTab[i++] = randomGenerator.nextFloat();
172                        impulseTab[i++] = randomGenerator.nextFloat();
173                        impulseTab[i++] = randomGenerator.nextFloat();
174                        impulseTab[i++] = 1.0f - 2.0f*randomGenerator.nextFloat();
175                }
176                
177                return impulseTab;
178        }
179        
180}