fournier@Navajo.ARPA (07/31/85)
/* computes approximations to 3d fractional Brownian motion */
/* try as input: 4 0.65 1.0 4157633211 */
/* to watch results, easiest is to try 64*64*64, then load the resulting
*/
/* array into a frame buffer as 64 64x64 frames. Then use the hardware */
/* pan/scroll/zoom feature (if your frame buffer does not have it, junk
*/
/* it) to "slice" through the cube. You will get a nice animation */
/* of clouds or fire, depending on the colour map you use. You might */
/* have to rescale the numbers you write out to fit your frame buffer */
#include <stdio.h>
#define TYPES 7
#define SIZE 65
/* At the beginning there are 8 CORNERS */
/* Then are computed 8 types of points. The first type is CENTER,
computed from */
/* the four CORNERS. The next three are FACE CENTERS, computed from two
CENTERS */
/* and four CORNERS. There are three kinds of FACE CENTERS, in order */
/* in the XY plane, in the XZ plane and in the YZ plane. */
/* Last are EDGE CENTERS, computed from two CORNERS and four FACE
CENTERS. */
/* There are three kinds of EDGE CENTERS, in order Z//, X// and Y// */
float fbm[SIZE][SIZE][SIZE];
int start[TYPES][3] = {
1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0
};
int neigh[TYPES] = {
8, 6, 6, 6, 6, 6, 6
};
int offset[TYPES][8][3] = {
1, 1, 1, 1, 1, -1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, -1, -1, -1, 1, -1, -1, -1,
0, 0, 1, 0, 0, -1, 1, 1, 0, 1, -1, 0, -1, 1, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, -1, 0, 1, 0, 1, 1, 0, -1, -1, 0, 1, -1, 0, -1, 0, 0, 0, 0, 0, 0,
1, 0, 0, -1, 0, 0, 0, 1, 1, 0, 1, -1, 0, -1, 1, 0, -1, -1, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, -1, 0, -1, 0, 0, 1, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0,
1, 0, 0, -1, 0, 0, 0, 0, -1, 0, 0, 1, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, -1, 0, 0, 0, 1, 0, 0, -1, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0
};
/* the last point in the loops is U-start for each coordinate */
/* the steps are 2*dist in every case */
/* the orders of the points is important, since they use some of the
precedings */
float weight[TYPES][8] = {
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0,
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0,
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0,
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0,
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0,
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0
};
/* for the time being, the weights are all one. They should be related to
the distance */
int fudge[3] = {
2514521, 891653, 8765142
};
/* fudge(s) are coefficients to relate the pseudo-random seeds to the indices
*/
main () {
float h,
scale,
std,
ratio;
int level,
seed,
max,
clevel,
type,
i,
j,
k,
dist;
int inseed;
float gauss (), expect ();
int scrseed ();
int scaledvalue;
double pow ();
int fd;
printf (" Give maximum level, h, scale and seed: ");
scanf ("%d %f %f %d", &level, &h, &scale, &inseed);
if (level < 0 || level > 6) {
fprintf (stderr, " Level (%d) out of range.\n", level);
exit ();
};
if (h > 1.0 || h < 0.0) {
fprintf (stderr, " h (%f) out of range.\n", h);
exit ();
};
if (scale < 0.0) {
fprintf (stderr, " Scale (%f) out of range.\n", scale);
exit ();
};
/* determine max */
max = 1 << level;
std = scale;
dist = max;
ratio = pow (2.0, -h);
/* compute corners (all eight of them) */
seed = scrseed (0, 0, 0, inseed, max);
fbm[0][0][0] = gauss (0.0, std, &seed);
seed = scrseed (max, 0, 0, inseed, max);
fbm[max][0][0] = gauss (0.0, std, &seed);
seed = scrseed (0, max, 0, inseed, max);
fbm[0][max][0] = gauss (0.0, std, &seed);
seed = scrseed (0, 0, max, inseed, max);
fbm[0][0][max] = gauss (0.0, std, &seed);
seed = scrseed (max, max, 0, inseed, max);
fbm[max][max][0] = gauss (0.0, std, &seed);
seed = scrseed (0, max, max, inseed, max);
fbm[0][max][max] = gauss (0.0, std, &seed);
seed = scrseed (max, 0, max, inseed, max);
fbm[max][0][max] = gauss (0.0, std, &seed);
seed = scrseed (max, max, max, inseed, max);
fbm[max][max][max] = gauss (0.0, std, &seed);
/* compute all others */
for (clevel = 1; clevel <= level; clevel++) {
std = std * ratio;
dist = dist / 2;
for (type = 0; type < TYPES; type++) {
for (i = start[type][0] * dist; i <= max - start[type][0] * dist; i += 2 * dist) {
for (j = start[type][1] * dist; j <= max - start[type][1] * dist; j += 2 * dist) {
for (k = start[type][2] * dist; k <= max - start[type][2] * dist; k += 2 * dist) {
/* compute value */
seed = scrseed (i, j, k, inseed, max);
fbm[i][j][k] = gauss (expect (type, i, j, k, dist, max), std, &seed);
};
};
};
};
};
/* save to file */
/* replace name "fbm" by your favourite string */
fd = creat ("fbm", 2);
if (fd < 1) {
fprintf (stderr, "Can't open file. \n");
exit ();
};
write (fd, &max, sizeof (int));
for (k = 0; k < max; k++) {
for (i = 0; i < max; i++) {
for (j = 0; j < max; j++) {
/* somewhat arbitrary scaling, based on -2:+2 range */
scaledvalue = 2048 + 1024 * fbm[i][j][k];
write (fd, &scaledvalue, sizeof (int));
};
};
};
close (fd);
}
float expect (type, i, j, k, dist, max)
/* computes the expected values, given the
type of point and the distances */
int type,
i,
j,
k,
dist,
max;
{
float exp;
int in,
jn,
kn,
neighbour;
float w,
sumw;
exp = 0.0;
sumw = 0.0;
for (neighbour = 0; neighbour < neigh[type]; neighbour++) {
/* if an index is out of range, take the mirror image with negative
weight */
w = weight[type][neighbour];
in = i + dist * offset[type][neighbour][0];
if (in > max || in < 0) {
w = -w;
in = i - dist * offset[type][neighbour][0];
};
jn = j + dist * offset[type][neighbour][1];
if (jn > max || jn < 0) {
w = -w;
jn = j - dist * offset[type][neighbour][1];
};
kn = k + dist * offset[type][neighbour][2];
if (kn > max || kn < 0) {
w = -w;
kn = k - dist * offset[type][neighbour][2];
};
/* above, only one of the three indices can be out of range */
exp = exp + fbm[in][jn][kn] * w;
sumw = sumw + w;
};
return exp / sumw;
}
float gauss (mean, std, seed) /* note that std is standard dev., not
variance */
/* also seed is pointer, since through ran it returns a new value */
/* you should substitute a better normal distribution generator */
float mean,
std;
int *seed;
{
float sum;
float ran ();
int i;
sum = 0.0;
for (i = 1; i <= 12; i++)
sum = sum + ran (seed);
return (mean + (sum - 6.0) * std);
}
float ran (seed)
int *seed;
/* use your own function here (should return uniformly distributed */
/* float between 0 and 1) */
{
int max = 017777777777;
srand (*seed);
*seed = rand ();
return (*seed / (float) max);
}
int scrseed (i, j, k, inseed, max)
int i,
j,
k,
inseed,
max;
/* compute seed as function of indices */
{
return (i * fudge[0] + j * fudge[1] + (k % max) * fudge[2] + inseed);
}