00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
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
00046
00047
00048
00049
00050
00051 #define SMOOTHSTEP(x) ((x) * (x) * (3 - 2*(x)))
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
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
00078 if (src[i] < 0) dst[i] = -src[i];
00079 else dst[i] = src[i];
00080
00081
00082
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
00108
00109
00110
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
00122
00123
00124
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
00156
00157
00158
00159
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
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
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
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233 double
00234 bn_noise_perlin(fastf_t *point)
00235 {
00236 register int jx, jy, jz;
00237 int ix, iy, iz;
00238 double x, y, z;
00239 double fx, fy, fz;
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
00250
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;
00265 jy = iy + 1;
00266 jz = iz + 1;
00267
00268 sx = SMOOTHSTEP(fx);
00269 sy = SMOOTHSTEP(fy);
00270 sz = SMOOTHSTEP(fz);
00271
00272
00273 tx = 1.0 - sx;
00274 ty = 1.0 - sy;
00275 tz = 1.0 - sz;
00276
00277
00278
00279
00280
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
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;
00317 double x, y, z;
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
00329
00330
00331
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
00349 tx = 1.0 - sx;
00350 ty = 1.0 - sy;
00351 tz = 1.0 - sz;
00352
00353
00354
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
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
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
00456
00457
00458
00459
00460
00461
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
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
00483 for (frequency = 1.0, i=0; i < octaves; i++) {
00484
00485 spec_wgts[i] = pow(frequency, -h_val);
00486 frequency *= lacunarity;
00487 }
00488
00489 etbl_next++;
00490 return ep;
00491 }
00492
00493
00494
00495
00496
00497
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
00518
00519
00520
00521 bu_semaphore_acquire(BU_SEM_BN_NOISE);
00522
00523
00524
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
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
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
00576
00577
00578
00579
00580
00581 ep = find_spec_wgt(h_val, lacunarity, octaves);
00582
00583
00584
00585 value = 0.0;
00586
00587 PCOPY(pt, point);
00588
00589 spec_wgts = ep->spec_wgts;
00590
00591
00592 oct=(int)octaves;
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
00601
00602
00603 value += noise_remainder * bn_noise_perlin(pt) * spec_wgts[i];
00604 }
00605
00606 return value;
00607
00608 }
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
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
00647
00648
00649
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
00658
00659 value = 0.0;
00660
00661
00662
00663
00664 PCOPY(pt, point);
00665 spec_wgts = ep->spec_wgts;
00666
00667
00668 oct=(int)octaves;
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
00677
00678
00679 value += noise_remainder * bn_noise_perlin(pt) * spec_wgts[i];
00680 }
00681 #else
00682 PCOPY(pt, point);
00683
00684 value = 0.0;
00685 frequency = 1.0;
00686
00687 oct=(int)octaves;
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
00697
00698
00699 value += noise_remainder * bn_noise_perlin(pt) * pow(frequency, -h_val);
00700 }
00701 #endif
00702 return value;
00703
00704 }
00705
00706
00707
00708
00709
00710
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
00721
00722
00723
00724
00725
00726 ep = find_spec_wgt(h_val, lacunarity, octaves);
00727
00728
00729
00730
00731 PCOPY(pt, point);
00732 spec_wgts = ep->spec_wgts;
00733
00734
00735
00736 noise_signal = bn_noise_perlin(pt);
00737
00738
00739 if (noise_signal < 0.0) noise_signal = -noise_signal;
00740
00741
00742 noise_signal = offset - noise_signal;
00743
00744
00745 noise_signal *= noise_signal;
00746
00747
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
00760 noise_signal *= weight;
00761 result += noise_signal * spec_wgts[i];
00762 }
00763 return result;
00764 }
00765
00766
00767
00768
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
00780
00781
00782
00783
00784
00785 ep = find_spec_wgt(h_val, lacunarity, octaves);
00786
00787
00788
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
00814
00815
00816
00817
00818
00819
00820