Changeset 3726


Ignore:
Timestamp:
Dec 25, 2014, 4:59:34 PM (7 years ago)
Author:
sam
Message:

math: shitloads of tweaks, optimisations, fixes and comments to the
simplex noise code.

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/demos/test/simplex.cpp

    r3724 r3726  
    1818using namespace lol;
    1919
     20/* Constants to tweak */
     21float const zoom = 0.03f;
     22int const octaves = 1;
     23
    2024int main(int argc, char **argv)
    2125{
    2226    UNUSED(argc, argv);
    23 
    24     float const zoom = 0.03f;
    25     int const octaves = 10;
    2627
    2728    /* Create an image */
     
    3031    array2d<vec4> &data = img.Lock2D<PixelFormat::RGBA_F32>();
    3132
    32     /* Fill image with simplex noise */
     33    /* Declare plenty of allocators */
    3334    simplex_interpolator<2> s2;
    3435    simplex_interpolator<3> s3;
     36    simplex_interpolator<4> s4;
     37    simplex_interpolator<5> s5;
     38    simplex_interpolator<6> s6;
    3539    simplex_interpolator<7> s7;
     40    simplex_interpolator<8> s8;
     41    simplex_interpolator<9> s9;
     42    simplex_interpolator<10> s10;
     43    simplex_interpolator<11> s11;
     44    simplex_interpolator<12> s12;
    3645
     46    /* Fill image with simplex noise */
    3747    for (int j = 0; j < size.y; ++j)
    3848    for (int i = 0; i < size.x; ++i)
    3949    {
     50        int cell = j / (size.y / 2) * 3 + i / (size.x / 3);
     51
    4052        float x = (float)i, y = (float)j;
    41         float sum = 0.f;
    42         int maxoctave = (j < size.y / 2) ? 1 : (1 << octaves);
     53        float sum = 0.f, coeff = 0.f;
     54        int const max_k = 1 << octaves;
    4355
    44         for (int k = 1; k <= maxoctave; k *= 2)
     56        for (int k = 1; k < max_k; k *= 2)
    4557        {
    46             if (i < size.x / 3)
    47                 sum += 0.5f / k * s2.Interp(zoom * k * vec2(x, y));
    48             else if (i < size.x * 2 / 3)
    49                 sum += 0.5f / k * s3.Interp(zoom * k * vec3(x, y, 0.f));
    50             else
    51                 sum += 0.5f / k * s7.Interp(zoom * k * vec7(x, 0.f, 0.f, 0.f, y, 0.f, 0.f));
     58            float t = 0.f;
     59
     60            switch (cell)
     61            {
     62            case 0:
     63                t = s2.Interp(zoom * k * vec2(x, y));
     64                break;
     65            case 1:
     66                t = s3.Interp(zoom * k * vec3(x, 1.f, y));
     67                break;
     68            case 2:
     69                t = s4.Interp(zoom * k * vec4(x, 1.f, y, 2.f));
     70                break;
     71            case 3:
     72                t = s6.Interp(zoom * k * vec6(x, 1.f, 2.f, y, 3.f, 4.f));
     73                break;
     74            case 4:
     75                t = s8.Interp(zoom * k * vec8(x, 1.f, 2.f, 3.f, y, 4.f,
     76                                              5.f, 6.f));
     77                break;
     78            case 5:
     79                t = s12.Interp(zoom * k * vec12(x, 1.f, 2.f, 3.f, 4.f, 5.f,
     80                                                y, 6.f, 7.f, 8.f, 9.f, 10.f));
     81                break;
     82            default:
     83                break;
     84            }
     85
     86            sum += 1.f / k * t;
     87            coeff += 1.f / k;
    5288        }
    5389
    54         float c = saturate(0.5f + sum);
     90        float c = saturate(0.5f + 0.5f * sum / coeff);
    5591        data[i][j] = vec4(c, c, c, 1.f);
     92        //data[i][j] = Color::HSVToRGB(vec4(c, 1.0f, 0.5f, 1.f));
    5693    }
     94
     95#if 0
     96    /* Mark simplex vertices */
     97    vec2 diagonal = normalize(vec2(1.f));
     98    vec2 dx = mat2::rotate(60.f) * diagonal;
     99    vec2 dy = mat2::rotate(-60.f) * diagonal;
     100    for (int j = -100; j < 100; ++j)
     101    for (int i = -100; i < 100; ++i)
     102    {
     103        auto putpixel = [&](ivec2 p, vec4 c = vec4(1.f, 0.f, 1.f, 1.f))
     104        {
     105            if (p.x >= 0 && p.x < size.x - 1 && p.y >= 0 && p.y < size.y - 1)
     106                data[p.x][p.y] = c;
     107        };
     108
     109        ivec2 coord = ivec2(i / zoom * dx + j / zoom * dy);
     110
     111        vec2 g = s2.GetGradient((i + 0.1f) * dx + (j + 0.1f) * dy);
     112        for (int t = 0; t < 40; ++t)
     113            putpixel(coord + (ivec2)(g * (t / 2.f)), vec4(0.f, 1.f, 0.f, 1.f));
     114
     115        putpixel(coord);
     116        putpixel(coord + ivec2(1, 0));
     117        putpixel(coord + ivec2(0, 1));
     118        putpixel(coord + ivec2(1, 1));
     119    }
     120#endif
    57121
    58122    /* Save image */
  • trunk/src/lol/math/simplex_interpolator.h

    r3722 r3726  
    1818{
    1919
     20/*
     21 * Simplex noise in dimension N
     22 * ----------------------------
     23 *
     24 *  The N-dimensional regular hypercube can be split into N! simplices that
     25 * all have the main diagonal as a shared edge.
     26 *  - number of simplices: N!
     27 *  - number of vertices per simplex: N+1
     28 *  - number of edges: N(N+1)/2
     29 *  - minimum edge length: 1           (hypercube edges, e.g. [1,0,0,…,0])
     30 *  - maximum edge length: sqrt(N)     (hypercube diagonal, i.e. [1,1,1,…,1])
     31 *
     32 *  We skew the simplicial grid along the main diagonal by a factor f =
     33 * sqrt(N+1), which means the diagonal of the initial parallelotope has
     34 * length sqrt(N/(N+1)). The edges of that parallelotope have length
     35 * sqrt(N/(N+1)), too. A formula for the maximum edge length was found
     36 * empirically:
     37 *  - minimum edge length: sqrt(N/(N+1))   (parallelotope edges and diagonal)
     38 *  - maximum edge length: sqrt(floor((N+1)²/4)/(N+1))
     39 */
     40
    2041template<int N>
    2142class simplex_interpolator
     
    2546      : m_seed(seed)
    2647    {
    27         /* Matrix coordinate transformation to skew simplex referential is done
    28          * by inverting the base matrix M which is written as follows:
    29          *
    30          *  M = | a b b b … |        M^(-1) = | c d d d … |
    31          *      | b a b b … |                 | d c d d … |
    32          *      | b b a b … |                 | d d c d … |
    33          *      | …         |                 | …         |
    34          *
    35          *  where a, b, c, d are computed below ↴
    36          */
    37 
    38         float b = (1.f - lol::sqrt(N + 1.f)) / lol::sqrt((float)N * N * N);
    39         float a = b + lol::sqrt((N + 1.f) / N);
    40 
    41         // determinant of matrix M
    42         float det = a * (a + (N - 2) * b) - (b * b) * (N - 1);
    43 
    44         float c = (a + (N - 2) * b) / det;
    45         float d = -b / det;
    46 
    47         for (int i = 0; i < N; ++i)
    48         {
    49             for (int j = 0; j < N; ++j)
    50             {
    51                 m_base[i][j] = (i == j ? a : b);
    52                 m_base_inverse[i][j] = (i == j ? c : d);
    53             }
     48#if 0
     49        // Print some debug information
     50        printf("Simplex Noise of Dimension %d\n", N);
     51
     52        long long int n = 1; for (int i = 1; i <= N; ++i) n *= i;
     53        printf(" - each hypercube cell has %lld simplices "
     54               "with %d vertices and %d edges\n", n, N + 1, N * (N + 1) / 2);
     55
     56        vec_t<float, N> diagonal(1.f);
     57        printf(" - length of hypercube edges is 1, diagonal is %f\n",
     58               length(diagonal));
     59        printf(" - length of parallelotope edges and diagonal is %f\n",
     60               length(unskew(diagonal)));
     61
     62        vec_t<float, N> vertices[N + 1];
     63        vertices[0] = vec_t<float, N>(0.f);
     64        for (int i = 0; i < N; ++i)
     65        {
     66            vertices[i + 1] = vertices[i];
     67            vertices[i + 1][i] += 1.f;
    5468        }
    55     }
    56 
    57     // Single interpolation
     69        float minlength = 10.0f;
     70        float maxlength = 0.0f;
     71        for (int i = 0; i < N + 1; ++i)
     72            for (int j = i + 1; j < N + 1; ++j)
     73            {
     74                float l = length(unskew(vertices[i] - vertices[j]));
     75                minlength = min(minlength, l);
     76                maxlength = max(maxlength, l);
     77            }
     78        printf(" - edge lengths between %f and %f\n", minlength, maxlength);
     79        printf(" - predicted: %f and %f\n", sqrt((float)N/(N+1)),
     80                                            sqrt((N+1)*(N+1)/4/(float)(N+1)));
     81        printf("\n");
     82#endif
     83    }
     84
     85    /* Single interpolation */
    5886    inline float Interp(vec_t<float, N> position) const
    5987    {
    60         // Computing position in simplex referential
    61         vec_t<float, N> simplex_position = m_base_inverse * position;
    62 
    63         // Retrieving the closest floor point and decimals
     88        // Retrieve the containing hypercube origin and associated decimals
    6489        vec_t<int, N> origin;
    6590        vec_t<float, N> pos;
    66 
    67         this->ExtractFloorDecimal(simplex_position, origin, pos);
    68 
    69         vec_t<int, N> index_order = this->GetIndexOrder(pos);
    70 
    71         // Resetting decimal point in regular orthonormal coordinates
    72         pos = m_base * pos;
    73 
    74         return this->LastInterp(origin, pos, index_order);
     91        get_origin(skew(position), origin, pos);
     92
     93        return get_noise(origin, pos);
     94    }
     95
     96    /* Only for debug purposes: return the gradient vector */
     97    inline vec_t<float, N> GetGradient(vec_t<float, N> position) const
     98    {
     99        vec_t<int, N> origin;
     100        vec_t<float, N> pos;
     101        get_origin(skew(position), origin, pos);
     102        return get_gradient(origin);
    75103    }
    76104
    77105protected:
    78     inline float LastInterp(vec_t<int, N> origin,
    79                             vec_t<float, N> const & pos,
    80                             vec_t<int, N> const & index_order) const
    81     {
    82         /* “corner” will traverse the simplex along its edges in
    83          * orthonormal coordinates. */
    84         vec_t<float, N> corner(0);
    85         float result = 0;
     106    inline float get_noise(vec_t<int, N> origin,
     107                           vec_t<float, N> const & pos) const
     108    {
     109        /* For a given position [0…1]^N inside a regular N-hypercube, find
     110         * the N-simplex which contains that position, and return a path
     111         * along the hypercube edges from (0,0,…,0) to (1,1,…,1) which
     112         * uniquely describes that simplex. */
     113        vec_t<int, N> traversal_order;
     114        for (int i = 0; i < N; ++i)
     115            traversal_order[i] = i;
     116
     117        /* Naïve bubble sort — enough for now */
     118        for (int i = 0; i < N; ++i)
     119            for (int j = i + 1; j < N; ++j)
     120                if (pos[traversal_order[i]] < pos[traversal_order[j]])
     121                    std::swap(traversal_order[i], traversal_order[j]);
     122
     123
     124        /* Get the position in world coordinates, too */
     125        vec_t<float, N> world_pos = unskew(pos);
     126
     127        /* “corner” will traverse the simplex along its edges in world
     128         * coordinates. */
     129        vec_t<float, N> world_corner(0.f);
     130        float result = 0.f, coeff = 0.f;
    86131
    87132        for (int i = 0; i < N + 1; ++i)
    88133        {
    89             vec_t<float, N> const delta = pos - corner;
    90             vec_t<float, N> const &gradient = GetGradient(origin);
    91 
    92             /* We use 4/3 because the maximum radius of influence for
    93              * a given corner is sqrt(3/4). FIXME: check whether this
    94              * holds for higher dimensions. */
    95             float dist = 1.0f - 4.f / 3.f * sqlength(delta);
    96             if (dist > 0)
    97             {
    98                 dist *= dist;
    99                 result += dist * dist * dot(gradient, delta);
     134#if 1
     135            // FIXME: In “Noise Hardware” (2-17) Perlin uses 0.6 but Gustavson
     136            // makes an exception for dimension 2 and uses 0.5 instead.
     137            // I think m should actually increase with the dimension count.
     138            float const m = N == 2 ? 0.5f : 0.6f;
     139            float d = m - sqlength(world_pos - world_corner);
     140#else
     141            // DEBUG: this is the linear contribution of each vertex
     142            // in the skewed simplex. Unfortunately it creates artifacts.
     143            float d = ((i == 0) ? 1.f : pos[traversal_order[i - 1]])
     144                    - ((i == N) ? 0.f : pos[traversal_order[i]]);
     145#endif
     146
     147            if (d > 0)
     148            {
     149                vec_t<float, N> const &gradient = get_gradient(origin);
     150
     151                // Perlin uses 8d⁴ whereas Gustavson uses d⁴ and a final
     152                // multiplication factor at the end.
     153                //d = 8.f * d * d * d * d;
     154                d = d * d * d * d;
     155
     156                //d = (3.f - 2.f * d) * d * d;
     157                //d = ((6 * d - 15) * d + 10) * d * d * d;
     158
     159                result += d * dot(gradient, world_pos - world_corner);
     160                coeff += d;
    100161            }
    101162
    102163            if (i < N)
    103164            {
    104                 corner += m_base[index_order[i]];
    105                 origin[index_order[i]] += 1;
     165                vec_t<float, N> v(0.f);
     166                v[traversal_order[i]] = 1.f;
     167                world_corner += unskew(v);
     168                origin[traversal_order[i]] += 1;
    106169            }
    107170        }
    108171
    109         // FIXME: another paper uses the value 70 for dimension 2, 32 for
    110         // dimension 3, and 27 for dimension 4; find where this comes from.
    111         return result * 70.f / 16.f;
    112     }
    113 
    114     /* For a given position [0…1]^n inside a square/cube/hypercube etc.,
    115      * find the simplex which contains that position, and return the path
    116      * from (0,0,…,0) to (1,1,…,1) that describes that simplex. */
    117     inline vec_t<int, N> GetIndexOrder(vec_t<float, N> const & pos) const
    118     {
    119         vec_t<int, N> result;
    120         for (int i = 0; i < N; ++i)
    121             result[i] = i;
    122 
    123         /* Naïve bubble sort — enough for now */
    124         for (int i = 0; i < N; ++i)
    125             for (int j = i + 1; j < N; ++j)
    126                 if (pos[result[i]] < pos[result[j]])
    127                     std::swap(result[i], result[j]);
    128 
    129         return result;
    130     }
    131 
    132     inline vec_t<float, N> GetGradient(vec_t<int, N> origin) const
     172        // FIXME: Gustavson uses the value 70 for dimension 2, 32 for
     173        // dimension 3, and 27 for dimension 4; find where this comes from
     174        // and maybe find a reasonable formula.
     175        //return 70.f * result / coeff;
     176        float const k = N == 2 ? 70 : N == 3 ? 32 : N == 4 ? 27 : 20;
     177        return k * result;
     178    }
     179
     180    inline vec_t<float, N> get_gradient(vec_t<int, N> origin) const
    133181    {
    134182        /* Quick shuffle table:
     
    137185        static int const shuffle[256] =
    138186        {
    139             111, 14, 180, 186, 221, 114, 17, 79, 66, 46, 11, 81, 246, 200, 141,
    140             172, 85, 244, 112, 92, 34, 106, 218, 205, 236, 7, 121, 115, 109,
    141             131, 10, 96, 188, 148, 219, 107, 94, 182, 235, 163, 143, 213, 248,
    142             202, 52, 154, 37, 241, 53, 129, 25, 60, 242, 38, 171, 63, 203, 255,
    143             193, 6, 42, 209, 28, 176, 210, 159, 54, 144, 3, 71, 89, 116, 12,
    144             237, 67, 216, 252, 178, 174, 164, 98, 234, 32, 26, 175, 24, 130,
    145             128, 113, 99, 212, 62, 152, 75, 185, 73, 93, 31, 30, 151, 122, 173,
    146             139, 91, 136, 162, 194, 208, 56, 101, 68, 69, 211, 44, 97, 55, 83,
    147             33, 50, 119, 156, 149, 41, 157, 253, 247, 161, 47, 230, 166, 225,
    148             204, 224, 13, 110, 123, 142, 64, 65, 155, 215, 120, 197, 140, 58,
    149             77, 214, 126, 195, 179, 220, 232, 125, 147, 8, 39, 187, 27, 217,
    150             100, 134, 199, 88, 206, 231, 250, 74, 2, 135, 9, 245, 118, 21, 243,
    151             82, 183, 238, 87, 158, 61, 4, 177, 146, 153, 117, 249, 254, 233,
    152             90, 222, 207, 48, 15, 18, 20, 239, 133, 0, 165, 138, 127, 169, 72,
    153             1, 201, 145, 191, 192, 16, 49, 19, 95, 226, 228, 84, 181, 251, 36,
    154             150, 22, 43, 70, 45, 105, 5, 189, 160, 196, 40, 59, 57, 190, 80,
    155             104, 167, 78, 124, 103, 240, 184, 170, 137, 29, 23, 223, 108, 102,
    156             86, 198, 227, 35, 229, 76, 168, 132, 51,
     187            111, 14, 180, 186, 221, 114, 219, 79, 66, 46, 152, 81, 246, 200,
     188            141, 172, 85, 244, 112, 92, 34, 106, 218, 205, 236, 7, 121, 115,
     189            109, 131, 10, 96, 188, 148, 17, 107, 94, 182, 235, 163, 143, 63,
     190            248, 202, 52, 154, 37, 241, 53, 129, 25, 159, 242, 38, 171, 213,
     191            6, 203, 255, 193, 42, 209, 28, 176, 210, 60, 54, 144, 3, 71, 89,
     192            116, 12, 237, 67, 216, 252, 178, 174, 164, 98, 234, 32, 26, 175,
     193            24, 130, 128, 113, 99, 212, 62, 11, 75, 185, 73, 93, 31, 30, 44,
     194            122, 173, 139, 91, 136, 162, 194, 41, 56, 101, 68, 69, 211, 151,
     195            97, 55, 83, 33, 50, 119, 156, 149, 208, 157, 253, 247, 161, 133,
     196            230, 166, 225, 204, 224, 13, 110, 123, 142, 64, 65, 155, 215, 9,
     197            197, 140, 58, 77, 214, 126, 195, 179, 220, 232, 125, 147, 8, 39,
     198            187, 27, 217, 100, 134, 199, 88, 206, 231, 250, 74, 2, 135, 120,
     199            21, 245, 118, 243, 82, 183, 238, 150, 158, 61, 4, 177, 146, 153,
     200            117, 249, 254, 233, 90, 222, 207, 48, 15, 18, 20, 16, 47, 0, 51,
     201            165, 138, 127, 169, 72, 1, 201, 145, 191, 192, 239, 49, 19, 160,
     202            226, 228, 84, 181, 251, 36, 87, 22, 43, 70, 45, 105, 5, 189, 95,
     203            40, 196, 59, 57, 190, 80, 104, 167, 78, 124, 103, 240, 184, 170,
     204            137, 29, 23, 223, 108, 102, 86, 198, 227, 35, 229, 76, 168, 132,
    157205        };
    158206
    159         /* 16 random vectors; this should be enough for small dimensions */
    160         auto v = [&]()
    161         {
    162             vec_t<float, N> ret;
    163             for (int i = 0; i < N; ++i)
    164                 ret[i] = rand(-1.f, 1.f);
    165             return normalize(ret);
     207        /* Generate 2^(N+2) random vectors, but at least 2^5 (32) and not
     208         * more than 2^20 (~ 1 million). */
     209        int const gradient_count = 1 << min(max(N + 2, 5), 20);
     210
     211        static auto build_gradients = [&]()
     212        {
     213            array<vec_t<float, N>> ret;
     214            for (int k = 0; k < gradient_count; ++k)
     215            {
     216                vec_t<float, N> v;
     217                for (int i = 0; i < N; ++i)
     218                    v[i] = rand(-1.f, 1.f);
     219                ret << normalize(v);
     220            }
     221            return ret;
    166222        };
    167223
    168         static vec_t<float, N> const gradients[16] =
    169         {
    170             v(), v(), v(), v(), v(), v(), v(), v(),
    171             v(), v(), v(), v(), v(), v(), v(), v(),
    172         };
     224        static array<vec_t<float, N>> const gradients = build_gradients();
    173225
    174226        int idx = m_seed;
    175227        for (int i = 0; i < N; ++i)
    176             idx = shuffle[(idx + origin[i]) & 255];
    177 
    178         return gradients[(idx ^ (idx >> 4)) & 15];
     228            idx ^= shuffle[(idx + origin[i]) & 255];
     229
     230        idx &= (gradient_count - 1);
     231#if 0
     232        // DEBUG: only output a few gradients
     233        if (idx > 2)
     234            return vec_t<float, N>(0);
     235#endif
     236        return gradients[idx];
     237    }
     238
     239    static inline vec_t<float, N> skew(vec_t<float, N> const &v)
     240    {
     241        /* Quoting Perlin in “Hardware Noise” (2-18):
     242         *   The “skew factor” f should be set to f = sqrt(N+1), so that
     243         *   the point (1,1,...1) is transformed to the point (f,f,...f). */
     244        float const sum = dot(v, vec_t<float, N>(1));
     245        float const f = sqrt(1.f + N);
     246        return v + vec_t<float, N>(sum * (f - 1) / N);
     247    }
     248
     249    static inline vec_t<float, N> unskew(vec_t<float, N> const &v)
     250    {
     251        float const sum = dot(v, vec_t<float, N>(1));
     252        float const f = sqrt(1.f + N);
     253        return v + vec_t<float, N>(sum * (1 / f - 1) / N);
    179254    }
    180255
    181256    /* For a given world position, extract grid coordinates (origin) and
    182257     * the corresponding delta position (pos). */
    183     inline void ExtractFloorDecimal(vec_t<float, N> const & simplex_position,
    184                                     vec_t<int, N> & origin,
    185                                     vec_t<float, N> & pos) const
     258    inline void get_origin(vec_t<float, N> const & simplex_position,
     259                           vec_t<int, N> & origin, vec_t<float, N> & pos) const
    186260    {
    187261        // Finding floor point index
     
    193267    }
    194268
    195     mat_t<float, N, N> m_base, m_base_inverse;
     269    /* A user-provided random seed. Defaults to zero. */
    196270    int m_seed;
    197271};
Note: See TracChangeset for help on using the changeset viewer.