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