1 /* fault.c: Drawing a fault in a height field
2 *
3 * Copyright (C) 2003-2004 Patrice St-Gelais
4 * patrstg@users.sourceforge.net
5 * www.oricom.ca/patrice.st-gelais
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20 */
21
22 #include "line.h"
23 #include "fault.h"
24 #include "hf.h"
25 #include "hf_calc.h"
26 #include "img_process.h"
27 #include "fill.h"
28
init_fault_buffer(fault_struct * f,gint hf_size)29 void init_fault_buffer (fault_struct *f, gint hf_size) {
30
31 // Make sure the fault buffer is large enough for the current HF
32 if (hf_size == f->buf_size)
33 return;
34 f->buf_size = hf_size;
35 f->buffer = (hf_type *)
36 x_realloc(f->buffer, sizeof(hf_type) * f->buf_size * f->buf_size, "hf_type (f->buffer in init_fault_buffer)");
37 }
38
prepare_fault_result_buf(fault_struct * f,hf_struct_type * hf)39 void prepare_fault_result_buf (fault_struct *f, hf_struct_type *hf) {
40 // Public function
41 // Used when beginning a fault, changing the altitude difference
42 // or separating planes (cracks)
43
44 hf_type diff;
45
46 // printf("PREPARE_RESULT_BUF\n");
47
48 memcpy(hf->result_buf, hf->tmp_buf, hf->max_x * hf->max_y * sizeof(hf_type));
49
50 // "Clamp" the HF values, removing the absolute altitude difference,
51 // so that we can merge - add the HF with the line buffer later
52
53 diff = (hf_type) ABS(( ( (gfloat) MAX_HF_VALUE) * (gfloat) f->altitude_difference) / 100.0);
54 if (diff < (MAX_HF_VALUE - hf->max)) // We have enough room
55 return;
56 if (diff < (MAX_HF_VALUE - hf->max + hf->min)) // We'll have enough room if we clamp the min to 0
57 hf_clamp_buffer (hf->result_buf, hf->max_x * hf->max_y, 0, hf->max - hf->min);
58 else
59 hf_clamp_buffer (hf->result_buf, hf->max_x * hf->max_y, 0, MAX_HF_VALUE - diff);
60
61 f->buffer_reset_required = FALSE;
62 }
63
begin_fault(fault_struct * f,hf_struct_type * hf)64 void begin_fault (fault_struct *f, hf_struct_type *hf) {
65
66 // "begin_fault" is called when the user clicks on the HF to draw the control line
67
68 // printf("BEGIN_FAULT\n");
69 hf_min_max (hf);
70
71 init_fault_buffer (f, hf->max_x);
72
73 // We need 2 result buffers
74 // The original HF (hf_buf) is copied in tmp_buf
75 // result_buf is allocated, if not already available
76 // It will contain the "clamped" tmp_buf, before the merge with the line buffer
77
78 if (!hf->tmp_buf)
79 hf_backup(hf);
80 if (!hf->result_buf) {
81 hf->result_buf = (hf_type *) x_malloc(hf->max_x*hf->max_y*sizeof(hf_type), "hf_type (hf->result_buf in begin_fault)");
82 }
83
84 prepare_fault_result_buf (f, hf);
85
86 }
87
compute_fault(fault_struct * f,hf_struct_type * hf,gint x0,gint y0,gint x1,gint y1,dist_matrix_struct * dm,gdouble ** gauss_list)88 void compute_fault (fault_struct *f, hf_struct_type *hf,
89 gint x0, gint y0, gint x1, gint y1,
90 dist_matrix_struct *dm, gdouble **gauss_list) {
91
92 // Computing the fault
93 // 1. Draw the polyline fault in the fault buffer
94 // (it was only drawn in a display buffer before)
95 // 2. Raise & lower the planes each side of the line, in the buffer
96 // This gives an "altitude difference map"
97 // 3. Smooth the buffer with the given radius
98 // 4. Merge the buffer with hf->result_buf, put the result in hf->hf_buf
99 // The calling function must unactivate the control line
100 // and activate the "accept" button
101
102 // f->mode = FAULT_PREVIEW (without smoothing) or FAULT_FINAL (with)
103
104 // static gchar *CORNERS[] = {"TOP_LEFT","TOP_RIGHT","BOTTOM_RIGHT",
105 // "BOTTOM_LEFT","LEFT_RIGHT","TOP_BOTTOM"}; // For testing purposes
106
107 gint diff, dx, dy, x, y, i, corner, minx, miny, maxx,maxy, ixs, iys;
108
109 glong lvalue;
110
111 gdouble ddist, startx=0.0, starty=0.0, endx=0.0, endy=0.0, ratio, dxs, dys;
112 // Ratio = ratio by which we multiply dx,dy to get the perpendicular separation
113
114 hf_type abs_diff, *sep_buf=NULL, *buf1=NULL, *buf2=NULL;
115
116 static hf_type value = 0x0001;
117
118 dx = x1-x0;
119 dy = y1-y0;
120
121 if ( (dx==dy) && !dx)
122 return;
123 // printf("COMPUTE_FAULT - MODE: %d\n",f->mode);
124
125 // printf("BUF_SIZE: %d; HF_SIZE: %d\n",f->buf_size, hf->max_x);
126
127 if (f->buffer_reset_required)
128 prepare_fault_result_buf (f,hf); // Resets the flag to FALSE
129 if (f->separation)
130 f->buffer_reset_required = TRUE; // For the next "compute_fault" call
131
132
133 // 1. Initialize the buffer to its default value
134 for (i=0; i<f->buf_size*f->buf_size; i++) {
135 *(f->buffer+i) = 0;
136 }
137
138 // 2. Extend the control line (mouse controlled)
139 // to find where the polyline intersects with the HF edges
140 if (! intersect_rect (0,0,hf->max_x-1, hf->max_y-1, x0, y0, x1, y1,
141 &startx, &starty, &endx, &endy) )
142 return;
143
144 ddist = DIST2(startx, starty, endx,endy);
145 set_line_scale (f->fractal_line->polyline, ddist, ddist);
146
147 // 3. Compute the absolute difference between the planes
148 diff = (gint) ( ( (gfloat) MAX_HF_VALUE) * (gfloat) f->altitude_difference) / 100.0;
149 abs_diff = (hf_type) ABS(diff);
150
151 // printf("\n***********LINE (%d,%d) - (%d,%d) from (%5.2f, %5.2f) to (%5.2f, %5.2f)\n",x0, y0, x1, y1, startx, starty, endx, endy);
152 // 4. Draw the polyline
153 draw_n_transform_all_segments (f->fractal_line->polyline,
154 endx-startx, endy-starty,
155 startx, starty,
156 (gpointer) f->buffer,
157 f->buf_size,
158 f->buf_size,
159 (gpointer) &value,
160 HF_TYPE_ID,
161 OVERFLOW_ZERO);
162
163 if (f->separation) {
164 sep_buf = (hf_type *) x_malloc(sizeof(hf_type) * f->buf_size * f->buf_size, "hf_type (sep_buf in compute_fault)");
165 memcpy(sep_buf, f->buffer,sizeof(hf_type) * f->buf_size * f->buf_size);
166 }
167
168 // 5. Apply the altitude difference
169 // Lower plane is always 0
170 // Difference ranges from -50% to +50%
171 // When difference is positive, the top or left plane is higher
172 // When difference is negative, the bottom or right plane is higher
173
174 miny = (gint) MIN(starty,endy);
175 minx = (gint) MIN(startx,endx);
176 maxy = (gint) MAX(starty+0.5,endy+0.5);
177 maxx = (gint) MAX(startx+0.5,endx+0.5);
178 if ( (!miny) && (!minx) )
179 corner = TOP_LEFT;
180 else
181 if ( (!miny) && (maxx==f->buf_size-1) )
182 corner = TOP_RIGHT;
183 else
184 if ( (!minx) && (maxy==f->buf_size-1) )
185 corner = BOTTOM_LEFT;
186 else
187 if ( (maxx==f->buf_size-1) &&
188 (maxy==f->buf_size-1) )
189 corner = BOTTOM_RIGHT;
190 else
191 if ( (!miny) && (maxy==f->buf_size-1) )
192 corner = LEFT_RIGHT;
193 else
194 corner = TOP_BOTTOM;
195
196 if ( (corner==TOP_LEFT) || (corner==LEFT_RIGHT) ||
197 (corner==TOP_BOTTOM) || (corner==BOTTOM_RIGHT) )
198 if (diff>0)
199 fill_area (f->buffer,f->buffer, f->buf_size, f->buf_size, abs_diff, 0, 0, 0);
200 else
201 fill_area (f->buffer, f->buffer, f->buf_size, f->buf_size, abs_diff, 0, f->buf_size-1, f->buf_size-1);
202 else
203 if (corner==TOP_RIGHT)
204 if (diff>0)
205 fill_area (f->buffer,f->buffer, f->buf_size, f->buf_size, abs_diff, 0, f->buf_size-1, 0);
206 else
207 fill_area (f->buffer, f->buffer, f->buf_size, f->buf_size, abs_diff, 0, 0, 0);
208 else // BOTTOM_LEFT
209 if (diff>0)
210 fill_area (f->buffer, f->buffer,f->buf_size, f->buf_size, abs_diff, 0, 0, f->buf_size-1);
211 else
212 fill_area (f->buffer,f->buffer, f->buf_size, f->buf_size, abs_diff, 0, f->buf_size-1,f->buf_size-1);
213
214 // printf( "(minx,miny): (%d,%d); (maxx,maxy): (%d,%d); corner: %s\n",minx,miny,maxx,maxy,CORNERS[corner]);
215
216 // Updated 2004-09-09 for processing separation (cracks)
217
218 // If f->mode == FAULT_PREVIEW, we merge the fault buffer and the HF, then separate
219 // If f->mode == FAULT_FINAL, we separate both the fault buffer and the HF,
220 // then smooth the fault buffer, then merge
221
222 if (f->separation) {
223 // 5.5 Calculate delta x and delta y resulting from the separation
224 // Initialize the separation buffer sep_buf with the polyline
225
226 // 1. Find the length of the vector perpendicular to the control line
227 // 2. Scale down this vector to f->separation, giving dxs,dys, in pixel units
228 ratio = ((gdouble) f->separation) / sqrt( (gdouble) (dx*dx + dy*dy) ) ;
229 dxs = ratio * (gdouble) dy;
230 dys = ratio * (gdouble) dx;
231
232 // printf("(dx,dy): (%d,%d); (dxs,dys): (%5.2f, %5.2f); ratio: %5.2f\n",dx,dy,dxs,dys,ratio);
233
234 // printf("(startx,starty): (%5.2f,%5.2f); (endx,endy): (%5.2f,%5.2f); CORNER: %s\n",startx,starty,endx,endy,CORNERS[corner]);
235
236 // For final computation, we try to get values in both planes, so that smoothing will work
237 if (f->mode==FAULT_FINAL)
238 for (i=0; i<f->buf_size*f->buf_size; i++)
239 *(f->buffer+i) += MAX_HF_VALUE>>2;
240
241 // We choose the corner to move so that the move doesn't create "blank" areas
242 // The delta sign depends on where the control line cuts the edges
243 switch (corner) {
244 case TOP_LEFT:
245 dxs = - ABS(dxs);
246 dys = - ABS(dys);
247 fill_area (sep_buf,sep_buf, f->buf_size, f->buf_size, value, 0, 0, 0);
248 break;
249 case TOP_RIGHT:
250 dxs = ABS(dxs);
251 dys = - ABS(dys);
252 fill_area (sep_buf,sep_buf, f->buf_size, f->buf_size, value, 0, f->buf_size-1, 0);
253 break;
254 case BOTTOM_LEFT:
255 dxs = - ABS(dxs);
256 dys = ABS(dys);
257 fill_area (sep_buf,sep_buf, f->buf_size, f->buf_size, value, 0, 0, f->buf_size-1);
258 break;
259 case BOTTOM_RIGHT:
260 dxs = ABS(dxs);
261 dys = ABS(dys);
262 fill_area (sep_buf,sep_buf, f->buf_size, f->buf_size, value, 0, f->buf_size-1, f->buf_size-1);
263 break;
264 case LEFT_RIGHT:
265 dxs = (gdouble) f->separation; // Temporary solution - we'll "stretch" one plane in the future
266 dys = 0.0;
267 fill_area (sep_buf,sep_buf, f->buf_size, f->buf_size, value, 0, f->buf_size-1, f->buf_size>>1);
268 break;
269 default: // TOP_BOTTOM
270 dys = (gdouble) f->separation;
271 dxs = 0.0;
272 fill_area (sep_buf,sep_buf, f->buf_size, f->buf_size, value, 0, f->buf_size>>1, f->buf_size-1);
273
274 }
275
276 // Integer values for preview (separation without interpolation)
277 ixs = (gint) floor(dxs+0.5);
278 iys = (gint) floor(dys+0.5);
279
280 // Separate the source and the control buffers, given sep_buf
281
282 // printf("(DX,DY): (%d, %d); (DXS, DYS): (%5.2f, %5.2f);\n", dx, dy, dxs,dys);
283 if (f->mode == FAULT_FINAL) {
284 // if (FALSE) {
285 buf1 = (hf_type *) x_calloc(f->buf_size*f->buf_size, sizeof(hf_type), "hf_type (buf1 in compute_fault)");
286 buf2 = (hf_type *) x_calloc(f->buf_size*f->buf_size, sizeof(hf_type), "hf_type (buf2 in compute_fault)");
287 // We separate then we merge
288 for (x=0; x < f->buf_size; x++) {
289 for (y=0; y< f->buf_size; y++) {
290 i =VECTORIZE(x,y,f->buf_size);
291 if (*(sep_buf+i) ) {
292 translate_real_forward_mapping (
293 f->buffer,
294 buf1,
295 HF_TYPE_ID,
296 f->buf_size, f->buf_size,
297 x,y,
298 dxs,dys);
299
300 translate_real_forward_mapping (
301 hf->result_buf,
302 buf2,
303 HF_TYPE_ID,
304 f->buf_size, f->buf_size,
305 x,y,
306 dxs,dys);
307
308 } // if sep_buf
309 else { // No translation!
310 *(buf1+i) = *(f->buffer+i);
311 *(buf2+i) = *(hf->result_buf+i);
312 } // if !sep_buf
313 } // for y
314 } // for x
315
316 swap_buffers ((gpointer *) &f->buffer, (gpointer *) &buf1);
317 swap_buffers ((gpointer *) &hf->result_buf, (gpointer *) &buf2);
318
319 } // End FAULT_FINAL
320
321 else {
322 // FAULT_PREVIEW
323 // We merge then we separate
324
325 // Merge
326 for (i=0; i<(f->buf_size*f->buf_size); i++)
327 *(hf->result_buf+i) = *(hf->result_buf+i) + *(f->buffer+i);
328 // We need an empty buffer
329 for (i=0; i<f->buf_size*f->buf_size; i++) {
330 *(hf->hf_buf+i) = 0;
331 }
332 // Separate (no interpolation)
333 for (x=0; x < f->buf_size; x++) {
334 for (y=0; y< f->buf_size; y++) {
335 i =VECTORIZE(x,y,f->buf_size);
336 if (*(sep_buf+i) ) {
337
338 translate_forward_mapping (
339 hf->result_buf,
340 hf->hf_buf,
341 HF_TYPE_ID,
342 f->buf_size, f->buf_size,
343 x,y,
344 ixs,iys);
345
346 } // if sep_buf
347 else { // No translation!
348 *(hf->hf_buf+i) = *(hf->result_buf+i);
349 } // if !sep_buf
350 } // for y
351 } // for x
352 } // End FAULT_PREVIEW
353
354 } // End with separation
355
356 else {
357
358 }
359
360
361 // 6. Smooth the line buffer
362 if (f->smoothing && (f->mode==FAULT_FINAL) ) {
363 swap_buffers ((gpointer*) &hf->hf_buf, (gpointer*) &f->buffer);
364 hf_smooth (hf, f->smoothing,dm, FALSE, gauss_list);
365 swap_buffers ((gpointer*) &hf->hf_buf, (gpointer*) &f->buffer);
366 }
367
368 // 7. Merge the line buffer with the HF
369
370 if (f->separation) {
371 if (f->mode==FAULT_FINAL) {
372 for (i=0; i<(f->buf_size*f->buf_size); i++) {
373 lvalue = ((glong) *(hf->result_buf+i)) + ((glong) *(f->buffer+i))
374 - (glong) (MAX_HF_VALUE>>2) ;
375 *(hf->hf_buf+i) = (hf_type) MAX(0,lvalue);
376 }
377 }
378 }
379 else
380 for (i=0; i<(f->buf_size*f->buf_size); i++)
381 *(hf->hf_buf+i) =*(hf->result_buf+i) + *(f->buffer+i);
382
383 // memcpy(hf->hf_buf, f->buffer, f->buf_size * f->buf_size * sizeof(hf_type)); // Test
384
385 // The next function draws the line in the buffer, for testing purposes
386 /*
387 draw_n_transform_all_segments (f->fractal_line,
388 endx-startx, endy-starty,
389 startx, starty,
390 (gpointer) hf->hf_buf,
391 f->buf_size,
392 f->buf_size,
393 (gpointer) &value,
394 HF_TYPE_ID,
395 OVERFLOW_ZERO);
396 */
397 if (sep_buf) {
398 x_free(sep_buf);
399 if (buf1)
400 x_free(buf1);
401 if (buf2)
402 x_free(buf2);
403 }
404 }
405
fault_pen_free(fault_struct * f)406 void fault_pen_free (fault_struct *f) {
407 if (f) {
408 if (f->fractal_line)
409 fractal_polyline_free (f->fractal_line);
410 x_free (f);
411 }
412 }
413
fault_pen_new()414 fault_struct *fault_pen_new() {
415 fault_struct *f;
416 gint i;
417 f = (fault_struct *) x_malloc(sizeof(fault_struct), "fault_struct");
418 f->fractal_line = fractal_polyline_new (7, // steps
419 0,50, // random_x, random_y
420 1,0, // width, random_w
421 123456, // seed
422 2.0, // distribution
423 FALSE, // if_cracks
424 50, 2, // cracks_depth, cracks_width
425 1); // cracks_branching_steps (for future enhancement)
426
427 f->smoothing = 4;
428 f->altitude_difference = 25;
429 f->separation = 0;
430 f->buffer_reset_required = TRUE;
431 f->mode = FAULT_FINAL;
432 f->buffer =(hf_type *) x_malloc(FAULT_PREVIEW_SIZE*FAULT_PREVIEW_SIZE * sizeof(hf_type), "hf_type (f->buffer in fault_pen_new)") ;
433 for (i=0; i<FAULT_PREVIEW_SIZE*FAULT_PREVIEW_SIZE; i++)
434 *(f->buffer+i) = MAX_HF_VALUE;
435 f->buf_size = FAULT_PREVIEW_SIZE;
436 return f;
437 }
438