1 /* Updated by Michel Wurtz 12/99 */
2
3 /*****************************************************************************/
4
5 /*** ***/
6
7 /*** spec_syn() ***/
8
9 /*** Creates a fractal surface using spectral synthesis. ***/
10
11 /*** Algorithm adapted from Peitgen and Sauper (1988), p.108. ***/
12
13 /*** Jo Wood, V1.0, 19th October, 1994 ***/
14
15 /*** Modified to allow multiple realisations of same surface, ***/
16
17 /*** Jo Wood, V1.1 15th October, 1995. ***/
18
19 /*** ***/
20
21 /*****************************************************************************/
22
23 #include <unistd.h>
24 #include <math.h>
25 #include <grass/gis.h>
26 #include <grass/glocale.h>
27 #include "frac.h"
28 #include <grass/gmath.h>
29
30
specsyn(double * data[2],int nn)31 int specsyn(double *data[2], /* Array holding complex data to transform. */
32 int nn /* Size of array in one dimension. */
33 )
34 {
35 /* ---------------------------------------------------------------- */
36 /* --- Initialise --- */
37 /* ---------------------------------------------------------------- */
38
39 int row, col, /* Counts through half the data array. */
40 row0, col0, /* 'other half' of the array. */
41 coeff; /* No. of Fourier coeffficents to calc. */
42
43 double phase, rad, /* polar coordinates of Fourier coeff. */
44 *temp[2];
45
46 /* FIXME - allow seed to be specified for repeatability */
47 G_math_srand_auto(); /* Reset random number generator. */
48
49 temp[0] = (double *)G_malloc(nn * nn * sizeof(double));
50 temp[1] = (double *)G_malloc(nn * nn * sizeof(double));
51
52 /* ---------------------------------------------------------------- */
53 /* --- Create fractal surface --- */
54 /* ---------------------------------------------------------------- */
55
56
57 /* Calculate all the preliminary random coefficients. */
58 /* ================================================== */
59
60 G_message(_("Preliminary surface calculations..."));
61 data_reset(data, nn);
62
63 for (row = 0; row <= nn / 2; row++)
64 for (col = 0; col <= nn / 2; col++) {
65 /* Generate random Fourier coefficients. */
66
67 phase = TWOPI * G_math_rand();
68
69 if ((row != 0) || (col != 0))
70 rad =
71 pow(row * row + col * col,
72 -(H + 1) / 2.0) * G_math_rand_gauss(1.);
73 else
74 rad = 0.0;
75
76 /* Fill half the array with coefficients. */
77
78 *(data[0] + row * nn + col) = rad * cos(phase);
79 *(data[1] + row * nn + col) = rad * sin(phase);
80
81 /* Fill other half of array with coefficients. */
82
83 if (row == 0)
84 row0 = 0;
85 else
86 row0 = nn - row;
87
88 if (col == 0)
89 col0 = 0;
90 else
91 col0 = nn - col;
92
93 *(data[0] + row0 * nn + col0) = rad * cos(phase);
94 *(data[1] + row0 * nn + col0) = -rad * sin(phase);
95 }
96
97 *(temp[1] + nn / 2) = 0.0;
98 *(temp[1] + nn * nn / 2) = 0.0;
99 *(temp[1] + nn * nn / 2 + nn / 2) = 0.0;
100
101 for (row = 1; row < nn / 2; row++)
102 for (col = 1; col < nn / 2; col++) {
103 phase = TWOPI * G_math_rand();
104 rad =
105 pow(row * row + col * col,
106 -(H + 1) / 2.0) * G_math_rand_gauss(1.);
107
108 *(data[0] + row * nn + nn - col) = rad * cos(phase);
109 *(data[1] + row * nn + nn - col) = rad * sin(phase);
110
111 *(data[0] + (nn - row) * nn + col) = rad * cos(phase);
112 *(data[1] + (nn - row) * nn + col) = -rad * sin(phase);
113 }
114
115 /* Transfer random coeffients to array before ifft transform */
116 /* ========================================================= */
117
118 for (coeff = 0; coeff < Steps; coeff++) {
119 G_message(_("Calculating surface %d (of %d)..."), coeff + 1, Steps);
120 data_reset(temp, nn);
121
122 for (row = 0; row <= (coeff + 1) * nn / (Steps * 2); row++)
123 for (col = 0; col <= (coeff + 1) * nn / (Steps * 2); col++) {
124 if (row == 0)
125 row0 = 0;
126 else
127 row0 = nn - row;
128
129 if (col == 0)
130 col0 = 0;
131 else
132 col0 = nn - col;
133
134 *(temp[0] + row * nn + col) = *(data[0] + row * nn + col);
135 *(temp[1] + row * nn + col) = *(data[1] + row * nn + col);
136
137 *(temp[0] + row0 * nn + col0) = *(data[0] + row0 * nn + col0);
138 *(temp[1] + row0 * nn + col0) = *(data[1] + row0 * nn + col0);
139 }
140
141 for (row = 1; row < (coeff + 1) * nn / (Steps * 2); row++)
142 for (col = 1; col < (coeff + 1) * nn / (Steps * 2); col++) {
143 *(temp[0] + row * nn + nn - col) =
144 *(data[0] + row * nn + nn - col);
145 *(temp[1] + row * nn + nn - col) =
146 *(data[1] + row * nn + nn - col);
147
148 *(temp[0] + (nn - row) * nn + col) =
149 *(data[0] + (nn - row) * nn + col);
150 *(temp[1] + (nn - row) * nn + col) =
151 *(data[1] + (nn - row) * nn + col);
152 }
153
154 fft(1, temp, nn * nn, nn, nn); /* Perform inverse FFT and */
155 write_rast(temp, nn, coeff + 1); /* write out raster. */
156 }
157
158 /* Free memory. */
159 /* ============ */
160
161 G_free(temp[0]);
162 G_free(temp[1]);
163
164 return 0;
165 }
166