noise.c

Go to the documentation of this file.
00001 /*                         N O I S E . C
00002  * BRL-CAD
00003  *
00004  * Copyright (c) 2004-2012 United States Government as represented by
00005  * the U.S. Army Research Laboratory.
00006  *
00007  * This library is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU Lesser General Public License
00009  * version 2.1 as published by the Free Software Foundation.
00010  *
00011  * This library is distributed in the hope that it will be useful, but
00012  * WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014  * Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public
00017  * License along with this file; see the file named COPYING for more
00018  * information.
00019  */
00020 /** @addtogroup noise */
00021 /** @{ */
00022 /** @file libbn/noise.c
00023  *
00024  * These noise functions provide mostly random noise at the integer
00025  * lattice points.  The functions should be evaluated at non-integer
00026  * locations for their nature to be realized.
00027  *
00028  * Conatins contributed code from:
00029  * F. Kenton Musgrave
00030  * Robert Skinner
00031  *
00032  */
00033 
00034 #include "common.h"
00035 
00036 #include <stdio.h>
00037 #include <stdlib.h>
00038 #include <math.h>
00039 #include "bu.h"
00040 #include "vmath.h"
00041 #include "bn.h"
00042 
00043 
00044 /**
00045  * @brief interpolate smoothly from 0 .. 1
00046  *
00047  * SMOOTHSTEP() takes a value in the range [0:1] and provides a number
00048  * in the same range indicating the amount of (a) present in a smooth
00049  * interpolation transition between (a) and (b)
00050  */
00051 #define SMOOTHSTEP(x) ((x) * (x) * (3 - 2*(x)))
00052 
00053 
00054 /**
00055  * @brief
00056  * fold space to extend the domain over which we can take noise
00057  * values.
00058  *
00059  *@n x, y, z are set to the noise space location for the source point.
00060  *@n ix, iy, iz are the integer lattice point (integer portion of x, y, z)
00061  *@n fx, fy, fz are the fractional lattice distance above ix, iy, iz
00062  *
00063  * The noise function has a finite domain, which can be exceeded when
00064  * using fractal textures with very high frequencies.  This routine is
00065  * designed to extend the effective domain of the function, albeit by
00066  * introducing periodicity. -FKM 4/93
00067  */
00068 static void
00069 filter_args(fastf_t *src, fastf_t *p, fastf_t *f, int *ip)
00070 {
00071     register int i;
00072     point_t dst;
00073     static unsigned long max2x = ~((unsigned long)0);
00074     static unsigned long max = (~((unsigned long)0)) >> 1;
00075 
00076     for (i=0; i < 3; i++) {
00077         /* assure values are positive */
00078         if (src[i] < 0) dst[i] = -src[i];
00079         else dst[i] = src[i];
00080 
00081 
00082         /* fold space */
00083         while (dst[i] > max || dst[i]<0) {
00084             if (dst[i] > max) {
00085                 dst[i] = max2x - dst[i];
00086             } else {
00087                 dst[i] = -dst[i];
00088             }
00089         }
00090 
00091     }
00092 
00093     p[X] = dst[0];
00094     ip[X] = floor(p[X]);
00095     f[X] = p[X] - ip[X];
00096 
00097     p[Y] = dst[1];
00098     ip[Y] = floor(p[Y]);
00099     f[Y] = p[Y] - ip[Y];
00100 
00101     p[Z] = dst[2];
00102     ip[Z] = floor(p[Z]);
00103     f[Z] = p[Z] - ip[Z];
00104 }
00105 
00106 
00107 #define MAXSIZE 267     /* 255 + 3 * (4 values) */
00108 
00109 /**
00110  * The RTable maps numbers into the range [-1..1]
00111  */
00112 static double RTable[MAXSIZE];
00113 
00114 #define INCRSUM(m, s, x, y, z)  ((s)*(RTable[m]*0.5             \
00115                                       + RTable[m+1]*(x)         \
00116                                       + RTable[m+2]*(y)         \
00117                                       + RTable[m+3]*(z)))
00118 
00119 
00120 /**
00121  * A heavily magic-number protected version of the hashtable.
00122  *
00123  * This table is used to convert integers into repeatable random
00124  * results for indicies into RTable.
00125  */
00126 struct str_ht {
00127     uint32_t magic;
00128     char hashTableValid;
00129     uint32_t *hashTableMagic1;
00130     uint32_t *hashTable;
00131     uint32_t *hashTableMagic2;
00132     uint32_t magic_end;
00133 };
00134 
00135 static struct str_ht ht;
00136 
00137 #define TABLE_SIZE 4096
00138 
00139 #define MAGIC_STRHT1 1771561
00140 #define MAGIC_STRHT2 1651771
00141 #define MAGIC_TAB1 9823
00142 #define MAGIC_TAB2 784642
00143 #define CK_HT() {                                                       \
00144         BU_CKMAG(&ht.magic, MAGIC_STRHT1, "struct str_ht ht 1");        \
00145         BU_CKMAG(&ht.magic_end, MAGIC_STRHT2, "struct str_ht ht 2");    \
00146         BU_CKMAG(ht.hashTableMagic1, MAGIC_TAB1, "hashTable Magic 1");  \
00147         BU_CKMAG(ht.hashTableMagic2, MAGIC_TAB2, "hashTable Magic 2");  \
00148         if (ht.hashTable != &ht.hashTableMagic1[1])                     \
00149             bu_bomb("ht.hashTable changed rel ht.hashTableMagic1");     \
00150         if (ht.hashTableMagic2 != &ht.hashTable[TABLE_SIZE])            \
00151             bu_bomb("ht.hashTable changed rel ht.hashTableMagic2");     \
00152     }
00153 
00154 /**
00155  * Map integer point into repeatable random number [0..TABLE_SIZE-1]
00156  * (i.e., [0-4095]).  We actually only use the first 8 bits of the
00157  * final value extracted from this table.  It's not quite clear that
00158  * we really need this big a table.  The extra size does provide some
00159  * extra randomness for intermediate results.
00160  */
00161 #define Hash3d(a, b, c)                                 \
00162     ht.hashTable[                                       \
00163         ht.hashTable[                                   \
00164             ht.hashTable[(a) & 0xfff] ^ ((b) & 0xfff)   \
00165             ]  ^ ((c) & 0xfff)                          \
00166         ]
00167 
00168 
00169 void
00170 bn_noise_init(void)
00171 {
00172     uint32_t i, j, k, temp;
00173     uint32_t rndtabi = BN_RAND_TABSIZE - 1;
00174 
00175     bu_semaphore_acquire(BU_SEM_BN_NOISE);
00176 
00177     if (ht.hashTableValid) {
00178         bu_semaphore_release(BU_SEM_BN_NOISE);
00179         return;
00180     }
00181 
00182     BN_RANDSEED(rndtabi, (BN_RAND_TABSIZE-1));
00183 
00184     /* alloc table size plus two magic numbers */
00185     ht.hashTableMagic1 = (uint32_t *) bu_calloc(1, 2*sizeof(uint32_t) + TABLE_SIZE*sizeof(uint32_t), "noise hashTable");
00186 
00187     ht.hashTable = &ht.hashTableMagic1[1];
00188     ht.hashTableMagic2 = &ht.hashTable[TABLE_SIZE];
00189 
00190     *ht.hashTableMagic1 = (uint32_t)MAGIC_TAB1;
00191     *ht.hashTableMagic2 = (uint32_t)MAGIC_TAB2;
00192 
00193     ht.magic_end = (uint32_t)MAGIC_STRHT2;
00194     ht.magic = (uint32_t)MAGIC_STRHT1;
00195 
00196     for (i = 0; i < TABLE_SIZE; i++)
00197         ht.hashTable[i] = i;
00198 
00199     /* scramble the hash table */
00200     for (i = TABLE_SIZE-1; i > 0; i--) {
00201         j = (uint32_t)(BN_RANDOM(rndtabi) * (fastf_t)TABLE_SIZE);
00202 
00203         temp = ht.hashTable[i];
00204         ht.hashTable[i] = ht.hashTable[j];
00205         ht.hashTable[j] = temp;
00206     }
00207 
00208     BN_RANDSEED(k, 13);
00209 
00210     for (i = 0; i < MAXSIZE; i++)
00211         RTable[i] = BN_RANDOM(k) * 2.0 - 1.0;
00212 
00213     ht.hashTableValid = 1;
00214 
00215     bu_semaphore_release(BU_SEM_BN_NOISE);
00216 
00217     CK_HT();
00218 }
00219 
00220 
00221 /**
00222  *@brief
00223  * Robert Skinner's Perlin-style "Noise" function
00224  *
00225  * Results are in the range [-0.5 .. 0.5].  Unlike many
00226  * implementations, this function provides random noise at the integer
00227  * lattice values.  However this produces much poorer quality and
00228  * should be avoided if possible.
00229  *
00230  * The power distribution of the result has no particular shape,
00231  * though it isn't as flat as the literature would have one believe.
00232  */
00233 double
00234 bn_noise_perlin(fastf_t *point)
00235 {
00236     register int jx, jy, jz;
00237     int ix, iy, iz;     /* lower integer lattice point */
00238     double x, y, z;     /* corrected point */
00239     double fx, fy, fz;  /* distance above integer lattice point */
00240     double sx, sy, sz, tx, ty, tz;
00241     double sum;
00242     short m;
00243     point_t p, f;
00244     int ip[3];
00245 
00246     if (!ht.hashTableValid)
00247         bn_noise_init();
00248 
00249     /* IS: const fastf_t *, point_t, point_t, int[3] */
00250     /* NE: fastf_t *, fastf_t *, fastf_t *, int * */
00251     filter_args(point, p, f, ip);
00252     ix = ip[X];
00253     iy = ip[Y];
00254     iz = ip[Z];
00255 
00256     fx = f[X];
00257     fy = f[Y];
00258     fz = f[Z];
00259 
00260     x = p[X];
00261     y = p[Y];
00262     z = p[Z];
00263 
00264     jx = ix + 1; /* (jx, jy, jz) = integer lattice point above (ix, iy, iz) */
00265     jy = iy + 1;
00266     jz = iz + 1;
00267 
00268     sx = SMOOTHSTEP(fx);
00269     sy = SMOOTHSTEP(fy);
00270     sz = SMOOTHSTEP(fz);
00271 
00272     /* the complement values of sx, sy, sz */
00273     tx = 1.0 - sx;
00274     ty = 1.0 - sy;
00275     tz = 1.0 - sz;
00276 
00277     /*
00278      * interpolate!
00279      */
00280     /* get a repeatable random # 0..TABLE_SIZE & 0xFF*/
00281     m = Hash3d(ix, iy, iz) & 0xFF;
00282     sum = INCRSUM(m, (tx*ty*tz), (x-ix), (y-iy), (z-iz));
00283 
00284     m = Hash3d(jx, iy, iz) & 0xFF;
00285     sum += INCRSUM(m, (sx*ty*tz), (x-jx), (y-iy), (z-iz));
00286 
00287     m = Hash3d(ix, jy, iz) & 0xFF;
00288     sum += INCRSUM(m, (tx*sy*tz), (x-ix), (y-jy), (z-iz));
00289 
00290     m = Hash3d(jx, jy, iz) & 0xFF;
00291     sum += INCRSUM(m, (sx*sy*tz), (x-jx), (y-jy), (z-iz));
00292 
00293     m = Hash3d(ix, iy, jz) & 0xFF;
00294     sum += INCRSUM(m, (tx*ty*sz), (x-ix), (y-iy), (z-jz));
00295 
00296     m = Hash3d(jx, iy, jz) & 0xFF;
00297     sum += INCRSUM(m, (sx*ty*sz), (x-jx), (y-iy), (z-jz));
00298 
00299     m = Hash3d(ix, jy, jz) & 0xFF;
00300     sum += INCRSUM(m, (tx*sy*sz), (x-ix), (y-jy), (z-jz));
00301 
00302     m = Hash3d(jx, jy, jz) & 0xFF;
00303     sum += INCRSUM(m, (sx*sy*sz), (x-jx), (y-jy), (z-jz));
00304 
00305     return sum;
00306 
00307 }
00308 
00309 /**
00310  * Vector-valued "Noise"
00311  */
00312 void
00313 bn_noise_vec(fastf_t *point, fastf_t *result)
00314 {
00315     register int jx, jy, jz;
00316     int ix, iy, iz;             /* lower integer lattice point */
00317     double x, y, z;             /* corrected point */
00318     double px, py, pz, s;
00319     double sx, sy, sz, tx, ty, tz;
00320     short m;
00321     point_t p, f;
00322     int ip[3];
00323 
00324 
00325     if (! ht.hashTableValid) bn_noise_init();
00326 
00327 
00328     /* sets:
00329      * x, y, z to range [0..maxval],
00330      * ix, iy, iz to integer portion,
00331      * fx, fy, fz to fractional portion
00332      */
00333     filter_args(point, p, f, ip);
00334     ix = ip[X];
00335     iy = ip[Y];
00336     iz = ip[Z];
00337 
00338     x = p[X];
00339     y = p[Y];
00340     z = p[Z];
00341 
00342     jx = ix+1;   jy = iy + 1;   jz = iz + 1;
00343 
00344     sx = SMOOTHSTEP(x - ix);
00345     sy = SMOOTHSTEP(y - iy);
00346     sz = SMOOTHSTEP(z - iz);
00347 
00348     /* the complement values of sx, sy, sz */
00349     tx = 1.0 - sx;
00350     ty = 1.0 - sy;
00351     tz = 1.0 - sz;
00352 
00353     /*
00354      * interpolate!
00355      */
00356     m = Hash3d(ix, iy, iz) & 0xFF;
00357     px = x-ix;
00358     py = y-iy;
00359     pz = z-iz;
00360     s = tx*ty*tz;
00361     result[0] = INCRSUM(m, s, px, py, pz);
00362     result[1] = INCRSUM(m+4, s, px, py, pz);
00363     result[2] = INCRSUM(m+8, s, px, py, pz);
00364 
00365     m = Hash3d(jx, iy, iz) & 0xFF;
00366     px = x-jx;
00367     s = sx*ty*tz;
00368     result[0] += INCRSUM(m, s, px, py, pz);
00369     result[1] += INCRSUM(m+4, s, px, py, pz);
00370     result[2] += INCRSUM(m+8, s, px, py, pz);
00371 
00372     m = Hash3d(jx, jy, iz) & 0xFF;
00373     py = y-jy;
00374     s = sx*sy*tz;
00375     result[0] += INCRSUM(m, s, px, py, pz);
00376     result[1] += INCRSUM(m+4, s, px, py, pz);
00377     result[2] += INCRSUM(m+8, s, px, py, pz);
00378 
00379     m = Hash3d(ix, jy, iz) & 0xFF;
00380     px = x-ix;
00381     s = tx*sy*tz;
00382     result[0] += INCRSUM(m, s, px, py, pz);
00383     result[1] += INCRSUM(m+4, s, px, py, pz);
00384     result[2] += INCRSUM(m+8, s, px, py, pz);
00385 
00386     m = Hash3d(ix, jy, jz) & 0xFF;
00387     pz = z-jz;
00388     s = tx*sy*sz;
00389     result[0] += INCRSUM(m, s, px, py, pz);
00390     result[1] += INCRSUM(m+4, s, px, py, pz);
00391     result[2] += INCRSUM(m+8, s, px, py, pz);
00392 
00393     m = Hash3d(jx, jy, jz) & 0xFF;
00394     px = x-jx;
00395     s = sx*sy*sz;
00396     result[0] += INCRSUM(m, s, px, py, pz);
00397     result[1] += INCRSUM(m+4, s, px, py, pz);
00398     result[2] += INCRSUM(m+8, s, px, py, pz);
00399 
00400     m = Hash3d(jx, iy, jz) & 0xFF;
00401     py = y-iy;
00402     s = sx*ty*sz;
00403     result[0] += INCRSUM(m, s, px, py, pz);
00404     result[1] += INCRSUM(m+4, s, px, py, pz);
00405     result[2] += INCRSUM(m+8, s, px, py, pz);
00406 
00407     m = Hash3d(ix, iy, jz) & 0xFF;
00408     px = x-ix;
00409     s = tx*ty*sz;
00410     result[0] += INCRSUM(m, s, px, py, pz);
00411     result[1] += INCRSUM(m+4, s, px, py, pz);
00412     result[2] += INCRSUM(m+8, s, px, py, pz);
00413 }
00414 /*************************************************************
00415  *@brief
00416  * Spectral Noise functions
00417  *
00418  *************************************************************
00419  *
00420  * The Spectral Noise functions cache the values of the term:
00421  *
00422  *@code
00423                     (-h_val)
00424                 freq
00425                 *@endcode
00426                 * Which on some systems is rather expensive to compute.
00427                 *
00428                 *************************************************************/
00429 struct fbm_spec {
00430     uint32_t magic;
00431     double octaves;
00432     double lacunarity;
00433     double h_val;
00434     double remainder;
00435     double *spec_wgts;
00436 };
00437 #define MAGIC_fbm_spec_wgt 0x837592
00438 
00439 static struct fbm_spec *etbl = (struct fbm_spec *)NULL;
00440 static int etbl_next = 0;
00441 static int etbl_size = 0;
00442 
00443 #define PSCALE(_p, _s) _p[0] *= _s; _p[1] *= _s; _p[2] *= _s
00444 #define PCOPY(_d, _s) _d[0] = _s[0]; _d[1] = _s[1]; _d[2] = _s[2]
00445 
00446 
00447 static struct fbm_spec *
00448 build_spec_tbl(double h_val, double lacunarity, double octaves)
00449 {
00450     struct fbm_spec *ep;
00451     double *spec_wgts;
00452     double frequency;
00453     int i;
00454 
00455     /* The right spectral weights table for these parameters has not
00456      * been pre-computed.  As a result, we compute the table now and
00457      * save it with the knowledge that we'll likely want it again
00458      * later.
00459      */
00460 
00461     /* allocate storage for new tables if needed */
00462     if (etbl_next >= etbl_size) {
00463         if (etbl_size) {
00464             etbl_size *= 2;
00465             etbl = (struct fbm_spec *)bu_realloc((char *)etbl,
00466                                                  etbl_size*sizeof(struct fbm_spec),
00467                                                  "spectral weights table");
00468         } else
00469             etbl = (struct fbm_spec *)bu_calloc(etbl_size = 10,
00470                                                 sizeof(struct fbm_spec),
00471                                                 "spectral weights table");
00472     }
00473 
00474     /* set up the next available table */
00475     ep = &etbl[etbl_next];
00476     ep->magic = MAGIC_fbm_spec_wgt;     ep->octaves = octaves;
00477     ep->h_val = h_val;          ep->lacunarity = lacunarity;
00478     spec_wgts = ep->spec_wgts =
00479         (double *)bu_malloc(((int)(octaves+1)) * sizeof(double),
00480                             "spectral weights");
00481 
00482     /* precompute and store spectral weights table */
00483     for (frequency = 1.0, i=0; i < octaves; i++) {
00484         /* compute weight for each frequency */
00485         spec_wgts[i] = pow(frequency, -h_val);
00486         frequency *= lacunarity;
00487     }
00488 
00489     etbl_next++; /* saved for last in case we're running multi-threaded */
00490     return ep;
00491 }
00492 
00493 /**
00494  * The first order of business is to see if we have pre-computed the
00495  * spectral weights table for these parameters in a previous
00496  * invocation.  If not, the we compute them and save them for possible
00497  * future use
00498  */
00499 struct fbm_spec *
00500 find_spec_wgt(double h, double l, double o)
00501 {
00502     struct fbm_spec *ep;
00503     int i;
00504 
00505 
00506     for (ep = etbl, i=0; i < etbl_next; i++, ep++) {
00507         if (ep->magic != MAGIC_fbm_spec_wgt)
00508             bu_bomb("find_spec_wgt");
00509         if (ZERO(ep->lacunarity - l)
00510             && ZERO(ep->h_val - h)
00511             && ep->octaves > (o - SMALL_FASTF))
00512         {
00513             return ep;
00514         }
00515     }
00516 
00517     /* we didn't find the table we wanted so we've got to semaphore on
00518      * the list to wait our turn to add what we want to the table.
00519      */
00520 
00521     bu_semaphore_acquire(BU_SEM_BN_NOISE);
00522 
00523     /* We search the list one more time in case the last process to
00524      * hold the semaphore just created the table we were about to add
00525      */
00526     for (ep = etbl, i=0; i < etbl_next; i++, ep++) {
00527         if (ep->magic != MAGIC_fbm_spec_wgt)
00528             bu_bomb("find_spec_wgt");
00529         if (ZERO(ep->lacunarity - l)
00530             && ZERO(ep->h_val - h)
00531             && ep->octaves > (o - SMALL_FASTF))
00532             break;
00533     }
00534 
00535     if (i >= etbl_next) ep = build_spec_tbl(h, l, o);
00536 
00537     bu_semaphore_release(BU_SEM_BN_NOISE);
00538 
00539     return ep;
00540 }
00541 /**
00542  * @brief
00543  * Procedural fBm evaluated at "point"; returns value stored in
00544  * "value".
00545  *
00546  * @param point          location to sample noise
00547  * @param ``h_val''      fractal increment parameter
00548  * @param ``lacunarity'' gap between successive frequencies
00549  * @param ``octaves''    number of frequencies in the fBm
00550  *
00551  * The spectral properties of the result are in the APPROXIMATE range
00552  * [-1..1] Depending upon the number of octaves computed, this range
00553  * may be exceeded.  Applications should clamp or scale the result to
00554  * their needs.  The results have a more-or-less gaussian
00555  * distribution.  Typical results for 1M samples include:
00556  *
00557  * @li Min -1.15246
00558  * @li Max  1.23146
00559  * @li Mean -0.0138744
00560  * @li s.d.  0.306642
00561  * @li Var 0.0940295
00562  *
00563  * The function call pow() is relatively expensive.  Therfore, this
00564  * function pre-computes and saves the spectral weights in a table for
00565  * re-use in successive invocations.
00566  */
00567 double
00568 bn_noise_fbm(fastf_t *point, double h_val, double lacunarity, double octaves)
00569 {
00570     struct fbm_spec *ep;
00571     double value, noise_remainder, *spec_wgts;
00572     point_t pt;
00573     int i, oct;
00574 
00575     /* The first order of business is to see if we have pre-computed
00576      * the spectral weights table for these parameters in a previous
00577      * invocation.  If not, the we compute them and save them for
00578      * possible future use
00579      */
00580 
00581     ep = find_spec_wgt(h_val, lacunarity, octaves);
00582 
00583     /* now we're ready to compute the fBm value */
00584 
00585     value = 0.0;            /* initialize vars to proper values */
00586     /* copy the point so we don't corrupt the caller's version */
00587     PCOPY(pt, point);
00588 
00589     spec_wgts = ep->spec_wgts;
00590 
00591     /* inner loop of spectral construction */
00592     oct=(int)octaves; /* save repeating double->int cast */
00593     for (i=0; i < oct; i++) {
00594         value += bn_noise_perlin(pt) * spec_wgts[i];
00595         PSCALE(pt, lacunarity);
00596     }
00597 
00598     noise_remainder = octaves - (int)octaves;
00599     if (!ZERO(noise_remainder)) {
00600         /* add in ``octaves'' noise_remainder ``i'' and spatial freq. are
00601          * preset in loop above
00602          */
00603         value += noise_remainder * bn_noise_perlin(pt) * spec_wgts[i];
00604     }
00605 
00606     return value;
00607 
00608 } /* bn_noise_fbm() */
00609 
00610 
00611 /**
00612  * @brief
00613  * Procedural turbulence evaluated at "point";
00614  *
00615  * @return turbulence value for point
00616  *
00617  * @param point          location to sample noise at
00618  * @param ``h_val''      fractal increment parameter
00619  * @param ``lacunarity'' gap between successive frequencies
00620  * @param ``octaves''    number of frequencies in the fBm
00621  *
00622  * The result is characterized by sharp, narrow trenches in low values
00623  * and a more fbm-like quality in the mid-high values.  Values are in
00624  * the APPROXIMATE range [0 .. 1] depending upon the number of octaves
00625  * evaluated.  Typical results:
00626  @code
00627  * Min  0.00857137
00628  * Max  1.26712
00629  * Mean 0.395122
00630  * s.d. 0.174796
00631  * Var  0.0305536
00632  @endcode
00633  * The function call pow() is relatively expensive.  Therfore, this
00634  * function pre-computes and saves the spectral weights in a table for
00635  * re-use in successive invocations.
00636  */
00637 double
00638 bn_noise_turb(fastf_t *point, double h_val, double lacunarity, double octaves)
00639 {
00640     struct fbm_spec *ep;
00641     double value, noise_remainder, *spec_wgts;
00642     point_t pt;
00643     int i, oct;
00644 
00645 
00646     /* The first order of business is to see if we have pre-computed
00647      * the spectral weights table for these parameters in a previous
00648      * invocation.  If not, the we compute them and save them for
00649      * possible future use
00650      */
00651 
00652 #define CACHE_SPECTRAL_WGTS 1
00653 #ifdef CACHE_SPECTRAL_WGTS
00654 
00655     ep = find_spec_wgt(h_val, lacunarity, octaves);
00656 
00657     /* now we're ready to compute the fBm value */
00658 
00659     value = 0.0;            /* initialize vars to proper values */
00660 
00661     /* copy the point so we don't corrupt the caller's copy of the
00662      * variable
00663      */
00664     PCOPY(pt, point);
00665     spec_wgts = ep->spec_wgts;
00666 
00667     /* inner loop of spectral construction */
00668     oct=(int)octaves; /* save repeating double->int cast */
00669     for (i=0; i < oct; i++) {
00670         value += fabs(bn_noise_perlin(pt)) * spec_wgts[i];
00671         PSCALE(pt, lacunarity);
00672     }
00673 
00674     noise_remainder = octaves - (int)octaves;
00675     if (!ZERO(noise_remainder)) {
00676         /* add in ``octaves'' noise_remainder ``i'' and spatial freq. are
00677          * preset in loop above
00678          */
00679         value += noise_remainder * bn_noise_perlin(pt) * spec_wgts[i];
00680     }
00681 #else
00682     PCOPY(pt, point);
00683 
00684     value = 0.0;      /* initialize vars to proper values */
00685     frequency = 1.0;
00686 
00687     oct=(int)octaves; /* save repeating double->int cast */
00688     for (i=0; i < oct; i++) {
00689         value += fabs(bn_noise_perlin(pt)) * pow(frequency, -h_val);
00690         frequency *= lacunarity;
00691         PSCALE(pt, lacunarity);
00692     }
00693 
00694     noise_remainder = octaves - (int)octaves;
00695     if (noise_remainder) {
00696         /* add in ``octaves'' noise_remainder ``i'' and spatial freq. are
00697          * preset in loop above
00698          */
00699         value += noise_remainder * bn_noise_perlin(pt) * pow(frequency, -h_val);
00700     }
00701 #endif
00702     return value;
00703 
00704 } /* bn_noise_turb() */
00705 
00706 /**
00707  *@brief
00708  * A ridged noise pattern
00709  *
00710  * From "Texturing and Modeling, A Procedural Approach" 2nd ed p338
00711  */
00712 double
00713 bn_noise_ridged(fastf_t *point, double h_val, double lacunarity, double octaves, double offset)
00714 {
00715     struct fbm_spec *ep;
00716     double result, weight, noise_signal, *spec_wgts;
00717     point_t pt;
00718     int i;
00719 
00720     /* The first order of business is to see if we have pre-computed
00721      * the spectral weights table for these parameters in a previous
00722      * invocation.  If not, the we compute them and save them for
00723      * possible future use
00724      */
00725 
00726     ep = find_spec_wgt(h_val, lacunarity, octaves);
00727 
00728     /* copy the point so we don't corrupt the caller's copy of the
00729      * variable
00730      */
00731     PCOPY(pt, point);
00732     spec_wgts = ep->spec_wgts;
00733 
00734 
00735     /* get first octave */
00736     noise_signal = bn_noise_perlin(pt);
00737 
00738     /* get absolute value of noise_signal (this creates the ridges) */
00739     if (noise_signal < 0.0) noise_signal = -noise_signal;
00740 
00741     /* invert and translate (note that "offset shoudl be ~= 1.0 */
00742     noise_signal = offset - noise_signal;
00743 
00744     /* square the noise_signal, to increase "sharpness" of ridges */
00745     noise_signal *= noise_signal;
00746 
00747     /* assign initial value */
00748     result = noise_signal;
00749     weight = 1.0;
00750 
00751     for (i=1; i < octaves; i++) {
00752         PSCALE(pt, lacunarity);
00753 
00754         noise_signal = bn_noise_perlin(pt);
00755 
00756         if (noise_signal < 0.0) noise_signal = - noise_signal;
00757         noise_signal = offset - noise_signal;
00758 
00759         /* weight the contribution */
00760         noise_signal *= weight;
00761         result += noise_signal * spec_wgts[i];
00762     }
00763     return result;
00764 }
00765 
00766 /**
00767  *
00768  * From "Texturing and Modeling, A Procedural Approach" 2nd ed
00769  */
00770 double
00771 bn_noise_mf(fastf_t *point, double h_val, double lacunarity, double octaves, double offset)
00772 {
00773     double frequency = 1.0;
00774     struct fbm_spec *ep;
00775     double result, weight, noise_signal, *spec_wgts;
00776     point_t pt;
00777     int i;
00778 
00779     /* The first order of business is to see if we have pre-computed
00780      * the spectral weights table for these parameters in a previous
00781      * invocation.  If not, the we compute them and save them for
00782      * possible future use
00783      */
00784 
00785     ep = find_spec_wgt(h_val, lacunarity, octaves);
00786 
00787     /* copy the point so we don't corrupt the caller's copy of the
00788      * variable.
00789      */
00790     PCOPY(pt, point);
00791     spec_wgts = ep->spec_wgts;
00792     offset = 1.0;
00793 
00794     result = (bn_noise_perlin(pt) + offset) * spec_wgts[0];
00795     weight = result;
00796 
00797     for (i=1; i < octaves; i++) {
00798         PSCALE(pt, lacunarity);
00799 
00800         if (weight > 1.0) weight = 1.0;
00801 
00802         noise_signal = (bn_noise_perlin(pt) + offset) * spec_wgts[i];
00803 
00804         noise_signal += fabs(bn_noise_perlin(pt)) * pow(frequency, -h_val);
00805         frequency *= lacunarity;
00806         PSCALE(pt, lacunarity);
00807     }
00808     return result;
00809 }
00810 
00811 /** @} */
00812 /*
00813  * Local Variables:
00814  * mode: C
00815  * tab-width: 8
00816  * indent-tabs-mode: t
00817  * c-file-style: "stroustrup"
00818  * End:
00819  * ex: shiftwidth=4 tabstop=8
00820  */
Generated on Tue Dec 11 13:14:27 2012 for LIBBN by  doxygen 1.6.3