1 /* subdiv1.c - subdivision algorithm #1:  cross / diamond
2  *
3  * Copyright (C) 2001 Patrice St-Gelais
4  *         patrstg@users.sourceforge.net
5  *
6  * This program is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19  */
20 
21 #include "hf.c"
22 #include "hf_calc.c"
23 #include "subdiv2.h"
24 
subdiv2_opt_new(gint * frq,gint nbrfrq)25 subdiv2_opt *subdiv2_opt_new(gint *frq,gint nbrfrq) {
26 	// frq: relative frequency level (from 0 to 100%)
27 	// nbrfrq: length of the frq vector, usually 1 or 12
28 	// 1: all the frequency level are the same
29 	// We use the values as a cycle, thanks to the modulo op.,
30 	// so it's not important if the value domain of nbrfrq is not enforced
31 	gint i;
32 	subdiv2_opt *opt;
33 	opt = (subdiv2_opt *) x_malloc(sizeof(subdiv2_opt), "subdiv2_opt");
34 	opt->seed = DEFAULT_SEED;
35 	opt->roughness = 0;
36 	opt->distribution = 1.0;
37 	opt->top_random = TRUE;
38 	opt->top_value = 0xF000;
39 	for (i=0; i<12; i++)
40 		opt->frq_percent[i] = *(frq+i%nbrfrq);
41 	return opt;
42 }
43 
44 
re_calc_2(hf_long val,hf_long * a,hf_long * b,hf_long * c,hf_long * d,int c_level)45 void re_calc_2( hf_long val, hf_long *a, hf_long *b,
46 		hf_long *c, hf_long *d, int c_level)
47 {
48 	static gfloat p1=0.2, p2 = 0.8;
49 
50 	val = (hf_long) (p1 * (gfloat) val);
51 	*(a) = val + (hf_long) (p2 * (gfloat) *(a));
52 	*(b) = val + (hf_long) (p2 * (gfloat) *(b));
53 	*(c) = val + (hf_long) (p2 * (gfloat) *(c));
54 	*(d) = val + (hf_long) (p2 * (gfloat) *(d));
55 }
56 
calc_r(gint level,gfloat froughness,gfloat distribution)57 static long int calc_r(gint level, gfloat froughness, gfloat distribution) {
58 	static long int r;
59 	static gfloat f;
60 	static gfloat maxf = (gfloat) 0xFFFF;
61 	r = 1 + (long int) ((hf_type) rand());
62 	if (distribution > 1.0) {
63 		f = (gfloat) r / maxf;
64 		r = (long int) (powf(f,distribution) * maxf);
65 	}
66 // printf("LEVEL: %d;  R1: %d;  ",level,r);
67 	r = r - (0xFFFF>>(level+2));
68 // printf("R2: %d;  ",r);
69  	r = (long int) (( (gfloat) r) * froughness);
70 // printf("R3: %d; level: %d; froughness: %f\n",r, level, froughness);
71 	return r;
72 }
73 
smooth_midpoints(hf_struct_type * hf,int x0,int y0,int x1,int y1,int c_level,hf_long * hfl)74 void smooth_midpoints (hf_struct_type *hf, int x0, int y0, int x1, int y1, int c_level, hf_long *hfl) {
75 	int mid_x, mid_y, x_1, y_1;
76 
77 	mid_x = WRAP((x0 + x1)>>1,hf->max_x);
78 	mid_y = WRAP((y0 + y1)>>1,hf->max_x);
79 
80 	x_1 = WRAP2(x0 - (mid_x - x0), hf->max_x);
81 	y_1 = WRAP2(mid_y - (y1 - y0),hf->max_y);
82 
83 	x0 = WRAP(x0,hf->max_x);
84 	y0 = WRAP(y0,hf->max_y);
85 	x1 = WRAP(x1,hf->max_y);
86 	y1 = WRAP(y1,hf->max_y);
87 
88 		re_calc_2(	*(hfl+VECTORIZE(mid_x,mid_y,hf->max_x))  ,
89 			hfl + VECTORIZE(x0,y0,hf->max_x),
90 			hfl + VECTORIZE(x0,y1,hf->max_x),
91 			hfl + VECTORIZE(x1,y0,hf->max_x),
92 			hfl + VECTORIZE(x1,y1,hf->max_x),
93 			c_level);
94 		re_calc_2(*(hfl+VECTORIZE(x0,mid_y,hf->max_x)) ,
95 			hfl +VECTORIZE(x0,y1,hf->max_x),
96 			hfl + VECTORIZE(mid_x,mid_y,hf->max_x),
97 			hfl + VECTORIZE(x0,y0,hf->max_x),
98 			hfl +VECTORIZE(x_1,mid_y,hf->max_x), c_level );
99 
100 		re_calc_2(*(hfl+VECTORIZE(mid_x,y0,hf->max_x)) ,
101 			hfl +VECTORIZE(mid_x,y_1,hf->max_x),
102 			hfl + VECTORIZE(x0,y0,hf->max_x),
103 			hfl + VECTORIZE(x1,y0,hf->max_x),
104 			hfl +VECTORIZE(mid_x,mid_y,hf->max_x), c_level );
105 }
106 
calc_midx_midy_2(hf_struct_type * hf,int x0,int y0,int x1,int y1,int c_level,subdiv2_opt * opt,gfloat froughness,hf_long * hfl)107 void calc_midx_midy_2(	hf_struct_type *hf, int x0, int y0, int x1, int y1, int c_level, subdiv2_opt *opt,
108 			gfloat froughness, hf_long *hfl) {
109 
110 	int mid_x, mid_y;
111 	hf_long avrg;
112 	long int val;
113 
114 	mid_x = WRAP((x0 + x1)>>1,hf->max_x);
115 	mid_y = WRAP((y0 + y1)>>1,hf->max_x);
116 	x0 = WRAP(x0,hf->max_x);
117 	y0 = WRAP(y0,hf->max_y);
118 	x1 = WRAP(x1,hf->max_y);
119 	y1 = WRAP(y1,hf->max_y);
120 
121 	avrg = calc_avrg (
122 		*(hfl + VECTORIZE(x0,y0,hf->max_x)),
123 		*(hfl +VECTORIZE(x0,y1,hf->max_x)),
124 		*(hfl +VECTORIZE(x1,y0,hf->max_x)),
125 		*(hfl +VECTORIZE(x1,y1,hf->max_x)) );
126 
127 	// mid_x, mid_y
128 /*
129 //	Disconnected, it doesn't work as intended, probably because of the smoothing process
130 	if ((!opt->top_random) && (c_level == 0)) {
131 // printf("TOP VALUE: %d; FRQ_PERCENT: %d\n", opt->top_value, opt->frq_percent[log2i(hf->max_x)-1]);
132 		*(hfl+VECTORIZE(mid_x,mid_y,hf->max_x)) =
133 				(hf_long) ( ( (gfloat) opt->frq_percent[log2i(hf->max_x)-1]/50.0)
134 				* (gfloat) opt->top_value) + (gfloat) ((~(hf_long) 0)>>25) ;
135 // printf("VALUE: %d\n",*(hfl+VECTORIZE(mid_x,mid_y,hf->max_x)) );
136 	  }
137 	else {
138 */
139 
140 	val = avrg OP_HF_2 calc_r (c_level, froughness, opt->distribution);
141 	*(hfl+VECTORIZE(mid_x,mid_y,hf->max_x))  = (hf_long) MAX(0,val);
142 
143 // if (c_level<2)
144 // printf("(mid_x,mid_y) == (%d,%d): AVRG: %u, VAL: %d; \n",mid_x,mid_y,avrg, val);
145 //	}
146 }
147 
calc_x_y_2(hf_struct_type * hf,int x0,int y0,int x1,int y1,int c_level,subdiv2_opt * opt,gfloat froughness,hf_long * hfl)148 void calc_x_y_2(	hf_struct_type *hf,
149 		int x0, int y0, int x1, int y1,
150 		int c_level,
151 		subdiv2_opt *opt,
152 		gfloat froughness,
153 		hf_long *hfl)
154 {
155 	int mid_x, mid_y, x_1, y_1;
156 	hf_long avrg;
157 	long int val;
158 
159 	mid_x = WRAP((x0 + x1)>>1,hf->max_x);
160 	mid_y = WRAP((y0 + y1)>>1,hf->max_x);
161 
162 	// Initialize height values - avrg taken from cross pos.
163 	x_1 = WRAP2(x0 - (mid_x - x0), hf->max_x);
164 	y_1 = WRAP2(mid_y - (y1 - y0),hf->max_y);
165 
166 	x0 = WRAP(x0,hf->max_x);
167 	y0 = WRAP(y0,hf->max_y);
168 	x1 = WRAP(x1,hf->max_x);
169 	y1 = WRAP(y1,hf->max_y);
170 
171 	// All points averaged from the mid points of the same level
172 	// All extremities are present, so we divide by 4 (>>2)
173 
174 	// For x, mid_y
175 	avrg = calc_avrg (*(hfl +( y1 * hf->max_x) + x0),
176 		*(hfl +( mid_y * hf->max_x) + mid_x),
177 		*(hfl + (y0 * hf->max_x) + x0),
178 		*(hfl +( mid_y * hf->max_x) + x_1 ));
179 
180 	// x0, mid_y
181 
182 	val = avrg OP_HF_2 calc_r (c_level, froughness,opt->distribution);
183 // if (c_level<2)
184 // printf("(x0,mid_y) == (%d,%d): AVRG: %u, VAL: %d; \n",x0,mid_y,avrg, val);
185 
186 	*(hfl+VECTORIZE(x0,mid_y,hf->max_x)) = (hf_long) MAX(val,0);
187 
188 	// For mid_x, y
189 	avrg = calc_avrg
190 		(*(hfl +( y_1 * hf->max_x) + mid_x),
191 		*(hfl +( y0 * hf->max_x) + x0),
192 		*(hfl + (y0 * hf->max_x) + x1),
193 		*(hfl +( mid_y * hf->max_x) + mid_x ));
194 
195 	// mid_x, y0
196 
197 	val = avrg OP_HF_2 calc_r (c_level, froughness, opt->distribution);
198 
199 // if (c_level<2)
200 // printf("(mid_x,y0) == (%d,%d): AVRG: %u, VAL: %d\n",mid_x,y0,avrg, val);
201 
202 	*(hfl+VECTORIZE(mid_x,y0,hf->max_x)) = (hf_long) MAX(0,val);
203 
204 }
205 
calc_level_2(hf_struct_type * hf,int level,subdiv2_opt * opt,hf_long * hfl)206 void calc_level_2(hf_struct_type *hf, int level, subdiv2_opt *opt, hf_long *hfl) {
207 	int chunk, i, j;
208 	gfloat froughness;
209 	chunk = hf->max_x >> level; // Works only for squares
210 
211 	froughness = ((gfloat) opt->frq_percent[log2i(hf->max_x)-level-1]/50.0) /
212 			powf(2.5,(gfloat) MAX(0.0,((gfloat) level) - 2.5 - ((gfloat) opt->roughness)/2.5)) ;
213 
214 // printf("LEVEL: %d; FROUGHNESS: %f\n",level,froughness);
215 	for (i = 0; i < hf->max_x; i= i+chunk) 			// diamond
216 		for (j = 0; j < hf->max_x; j = j + chunk)
217 			calc_midx_midy_2(hf, i, j, i+chunk, j+chunk,
218 				level, opt, froughness, hfl);
219 	for (i = 0 ; i < hf->max_x; i= i+chunk) 			// cross
220 		for (j = 0; j < hf->max_x; j = j + chunk)
221 			calc_x_y_2(hf, i, j, i+chunk, j+chunk,
222 				level, opt, froughness, hfl);
223 	for (i = 0; i <= hf->max_x; i= i+chunk) 			// smooth
224 		for (j = 0; j <= hf->max_x; j = j + chunk)
225 			smooth_midpoints(hf, i, j, i+chunk, j+chunk,
226 				level, hfl);
227 }
228 
subdiv2(hf_struct_type * hf,subdiv2_opt * opt)229 void subdiv2(hf_struct_type *hf, subdiv2_opt *opt) {
230 	// Plasma - style subdivision - TEST
231 	// Non recursive algorithm
232 	// For the given level:
233 	// 1. Calculate values for all mid points
234 	//	Average: square
235 	// 2. Calculate values for all mid_x, y0
236 	// 3. Calculate values for x0, mid_y
237 
238 	int i, levels;
239 	hf_long *hfl;
240 	hf_long min,max;
241 	gfloat ratio;
242 	srand(opt->seed);
243 
244 //	There could have been a size change, so we reallocate the memory
245 	hf->hf_buf = (hf_type *) x_realloc(hf->hf_buf, hf->max_x * hf->max_y * sizeof(hf_type), "hf_type (hf_buf in subdiv2)");
246 
247 //	long int buffer
248 	hfl = (hf_long *) x_malloc(hf->max_x * hf->max_y * sizeof(hf_long), "hf_long (hfl in subdiv2)");
249 
250 //	*(hf->hf_buf) = (hf_type) rand() ;
251 
252 //	Initialize the first pixel
253 	*(hfl) = (hf_long) ( ( (gfloat) opt->frq_percent[log2i(hf->max_x)-1]/50.0) *(hf_type) rand()) +
254 			 (gfloat) ((~(hf_long) 0)>>25) ;
255 
256 	levels = log2i(hf->max_x);
257 // printf("First pixel: %u; levels: %d\n",*(hfl), levels );
258 	for (i=0; i< levels; i++) {
259 		calc_level_2(hf, i, opt, hfl);
260 	}
261 //	Find max and min
262 	min = *(hfl);
263 	max = min;
264 	for (i=0; i<hf->max_x*hf->max_y; i++) {
265 		if (min>*(hfl+i) )
266 			min = *(hfl+i);
267 		else
268 			if (max<*(hfl+i))
269 				max = *(hfl+i);
270 	}
271 	if ((max-min)<=0xFFFF)
272 		ratio = 1.0;
273 	else
274 		ratio = ((gfloat) (max-min)) / (gfloat) 0xFFFF;
275 // printf("MIN: %u; MAX: %u; DIFF: %u;  RATIO: %f\n",min,max,max-min,ratio);
276 	for (i=0; i<hf->max_x*hf->max_y; i++) {
277 		*(hf->hf_buf+i) = (hf_type) (( (gfloat) (*(hfl+i)-min))/ratio);
278 	}
279 	x_free(hfl);
280 }
281 
282 
283 
284 
285