1 /* img_process.c - functions for filtering and transforming images
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 <math.h>
22 #include "img_process.h"
23 #include "hf.h"
24 #include "hf_calc.h"
25 #include "hf_filters.h"
26 #include "../fourier/fourier.h"
27
hf_auto_contrast(hf_struct_type * hf)28 void hf_auto_contrast (hf_struct_type *hf) {
29 // Stretch contrast so that the full value range is used (0 to 0xFFFF)
30 hf_clamp (hf, 0, 0xFFFF);
31 }
32
hf_brightness_contrast(hf_struct_type * hf,gint contrast_level,gint br_level,gboolean keep_luminosity,gboolean bound_overflow)33 void hf_brightness_contrast(hf_struct_type *hf, gint contrast_level, gint br_level,
34 gboolean keep_luminosity, gboolean bound_overflow) {
35 //
36 // Increases / decreases the contrast & brightness of a square size X size 16 greyscale image
37 // Levels are a % from -100 to 100
38 // The brightness process is additive
39 // At -100%, we subtract 0xFFFF>>1; at 100%, we add 0xFFFF>>1.
40 // The contrast is relative to the medium value between the min and max
41 // The input is hf->tmp_buf (backup of the displayed buffer before getting into
42 // the contrast dialog) and the output is the displayed buffer hf->hf_buf
43 // Min and Max should be updated BEFORE, when entering into the dialog
44
45 // The contrast process adjusts the delta between min and max value, so that
46 // the new delta is 0 at -100% level and 2x at 100% level
47
48 gfloat flevel;
49 gint i, max = (gint) 0xFFFF, max2 , value, offset, delta_brightness ;
50 max2 = 2 * max - 1;
51 if (!hf->tmp_buf)
52 hf_backup(hf);
53 delta_brightness = (gint) ( (0xFFFF>>1) * ((gfloat) br_level / 100));
54 // flevel ranges from 0.0 to 2.0
55 flevel = ((gfloat) (100 + contrast_level)) / 100 ;
56 if (keep_luminosity)
57 offset = hf->avrg * (1.0-flevel);
58 else
59 offset = 0;
60 for (i=0; i<hf->max_x*hf->max_y; i++) {
61 // If contrast_level < 0, we apply brightness 2nd
62 if (contrast_level < 0) {
63 if (bound_overflow) {
64 // Contrast
65 *(hf->hf_buf + i) = MAX(0,MIN(0xFFFF, offset
66 + ((gfloat) *(hf->tmp_buf + i)) * flevel));
67 // Brightness
68 if (br_level)
69 *(hf->hf_buf + i) =
70 (hf_type) MIN(0xFFFF,MAX(0,
71 (gint) *(hf->hf_buf + i) + delta_brightness));
72 }
73 else { // Reverted overflow
74 // Contrast
75 value = offset + ((gfloat) *(hf->tmp_buf+i)) * flevel;
76 // Brightness
77 if (br_level)
78 value = value + delta_brightness;
79 value = ABS(value);
80 value = (value>=max)?(max2-value):value;
81 *(hf->hf_buf+i) = (hf_type) value;
82 }
83 }
84 else {
85 if (bound_overflow) {
86 // Brightness
87 *(hf->hf_buf + i) =
88 (hf_type) MIN(0xFFFF,MAX(0,
89 (gint) *(hf->tmp_buf + i) + delta_brightness));
90 // Contrast
91 if (contrast_level)
92 *(hf->hf_buf + i) = MAX(0,MIN(0xFFFF, offset
93 + ((gfloat) *(hf->hf_buf + i)) * flevel));
94 }
95 else { // Reverted overflow
96 // Brightness
97 value = *(hf->tmp_buf+i) + delta_brightness;
98 // Contrast
99 if (contrast_level)
100 value = offset + ((gfloat) value) * flevel;
101 value = ABS(value);
102 value = (value>=max)?(max2-value):value;
103 *(hf->hf_buf+i) = (hf_type) value;
104 }
105 }
106 }
107 }
108
hf_remove_artifacts(hf_struct_type * hf,gint radius,dist_matrix_struct * dm)109 void hf_remove_artifacts (hf_struct_type *hf, gint radius, dist_matrix_struct *dm) {
110 // Removing artifacts in a quantized image (terrace)
111 // Move a window of a given radius over the image
112 // and fill it with the border color if all the border pixels are the same
113 // In the simplest form, the window is square
114 // The function "wraps"
115 gint i, j, k, l;
116 hf_type value;
117 gboolean test;
118 for (i=0; i<hf->max_x; i++) {
119 for (j=0; j<hf->max_y; j++) {
120 value = *(hf->hf_buf + i + j*hf->max_x);
121 test = TRUE;
122 for (k=-radius; k<radius+1; k++) {
123 if (value != *(hf->hf_buf + WRAP2(i - radius,hf->max_x)
124 + WRAP2((j+k),hf->max_y) *hf->max_x)) {
125 test = FALSE;
126 break;
127 }
128 if (value != *(hf->hf_buf + WRAP2(i + radius,hf->max_x)
129 + WRAP2((j+k),hf->max_y) *hf->max_x)) {
130 test = FALSE;
131 break;
132 }
133 } // End for k (1st)
134 if (!test)
135 continue;
136 for (k = -radius+1; k<radius; k++) {
137 if (value != *(hf->hf_buf + WRAP2(i + k,hf->max_x)
138 + WRAP2((j-radius),hf->max_y)*hf->max_x)) {
139 test = FALSE;
140 break;
141 }
142 if (value != *(hf->hf_buf + WRAP2(i + k,hf->max_x)
143 + WRAP2((j+radius),hf->max_y)*hf->max_x)) {
144 test = FALSE;
145 break;
146 }
147 } // End for k (2nd)
148 if (!test)
149 continue;
150 // We fill up the window with "value"
151 for (k=-radius+1; k<radius; k++)
152 for (l=-radius+1; l<radius; l++)
153 *(hf->hf_buf+VECTORIZE(
154 WRAP2(i+l,hf->max_x),
155 WRAP2(j+k,hf->max_y),
156 hf->max_x) )= value;
157
158 } // Enf for j
159 } // End for i
160 }
161
hf_terrace(hf_struct_type * hf,gint levels,gint seed,gint percent_random,gint smooth_radius,gboolean if_wrap,dist_matrix_struct * dist_matrix,gint artifact_radius,gboolean apply_postprocess,gdouble ** gauss_list)162 void hf_terrace (hf_struct_type *hf, gint levels, gint seed, gint percent_random,
163 gint smooth_radius, gboolean if_wrap,
164 dist_matrix_struct *dist_matrix,
165 gint artifact_radius, gboolean apply_postprocess,
166 gdouble **gauss_list) {
167
168 gint i,*disturbance_vector, index, slice, index_offset;
169
170 if (!levels)
171 return;
172 slice = (hf->max-hf->min) / levels;
173 // No process if slice is 0 - HF is uniform
174 if (!slice)
175 return;
176 // Modified 2005-06-23: s_offset and r_offset are now %,not absolute values
177 /*
178 s_offset = (slice* source_offset) / max_f; // source_offset and result_offset are % of slice
179 r_offset = (slice* result_offset) / max_f ;
180 printf("slice: %d; s_offset: %f; r_offset: %f\n",slice, s_offset, r_offset);
181 */
182 if (!hf->tmp_buf)
183 hf_backup(hf);
184 // We need a second temporary buffer for storing
185 // the result of the quantization process before the postprocessing
186 if (!hf->result_buf) {
187 hf->result_buf = (hf_type *)x_malloc(hf->max_x*hf->max_y*sizeof(hf_type), "hf_type (result_buf in hf_terrace)");
188 }
189 // Variation implementing random disturbance of each level
190 // Disturbance is bounded by 0 and slice_size
191 // (NOTE: the terrace effect is calculating by "flooring" a value
192 // to the nearest lower multiple of slice_size, so the total H
193 // should never be more than 0xFFFF)
194 /* disturbance_vector = (gint *) x_malloc(sizeof(gint)*(levels+1), "gint (disturbance_vector in hf_terrace)");
195 srand(seed);
196 *(disturbance_vector) = 0; // We want the "floor" to stay black
197 for (i=1; i<levels+1; i++) {
198 *(disturbance_vector+i) = (percent_random*(rand()%slice)) / 100;
199 }
200 for (i=0; i<hf->max_x*hf->max_y; i++) {
201 index = *(hf->tmp_buf + i) / slice;
202 *(hf->hf_buf + i) = ((gfloat) *(hf->tmp_buf +i)* (100-percent) / 100) +
203 ((gfloat) ( (*(disturbance_vector+index) +
204 (slice * index) )* percent) / 100);
205 }
206 */
207 // 2nd variation
208 // The disturbance vector is applied when quantizing, not after
209 disturbance_vector = (gint *) x_malloc(sizeof(gint)*(levels+1), "gint (disturbance_vector in hf_terrace)");
210 srand(seed);
211 *(disturbance_vector) = hf->min; // We want the "floor" to stay at minimum
212
213 for (i=1; i<levels+1; i++) {
214 *(disturbance_vector+i) = (percent_random*(rand()%slice)) / 100;
215 // For the quantization algorithm to work,
216 // each slice displacement must be strictly < to slice_size (never =)
217 *(disturbance_vector+i) = MIN(-1+((i+1)*slice), *(disturbance_vector+i));
218 // Contrary to the simpler preceding algorithm,
219 // the quantization requires cumulative values, so that we can
220 // rapidly check upper and lower bounds
221 *(disturbance_vector+i) = MIN(0xFFFF,hf->min + (i*slice)+*(disturbance_vector+i));
222 // printf("DIST-VECT[%d] = %d\n",i,*(disturbance_vector+i));
223 }
224
225 // When the minimum HF value is > slice size, we need to adjust the index
226 index_offset = hf->min / slice;
227 // If a smooth radius or an artifact removal radius is given,
228 // we must decompose the process in 2 passes:
229 if ( apply_postprocess && (smooth_radius || artifact_radius) ) {
230 for (i=0; i<hf->max_x*hf->max_y; i++) {
231 // First, we approximate quantization
232 index = (*(hf->tmp_buf + i) / slice) - index_offset;
233 // "index*slice_size" gives the 1st value to test
234 // Since displacement is strictly less than the slice_size,
235 // the looked for value can only be index or index-1
236 if (*(disturbance_vector+index)>*(hf->tmp_buf+i))
237 index = MAX(0,index-1); // This MAX is not supposed to be required!!
238 *(hf->hf_buf + i) = *(disturbance_vector+index);
239 }
240 // The 2 next processes operate on hf->hf_buf
241 if (artifact_radius)
242 hf_remove_artifacts(hf, artifact_radius, dist_matrix);
243 if (smooth_radius)
244 hf_smooth(hf, smooth_radius, dist_matrix, if_wrap, gauss_list);
245 // Swap result_buf and hf_buf, so that result_buf becomes the input
246 // and hf_buf is prepared for a new output
247 swap_buffers ((gpointer*) &hf->result_buf, (gpointer*) &hf->hf_buf);
248 } // end of "if (apply_postprocess)"
249 else { // if ( ! apply_postprocess)
250
251 for (i=0; i<hf->max_x*hf->max_y; i++) {
252 index = (*(hf->tmp_buf + i) / slice) - index_offset;
253 if (*(disturbance_vector+index)>*(hf->tmp_buf+i))
254 index = MAX(0,index-1);
255 *(hf->result_buf + i) = *(disturbance_vector+index);
256 }
257 }
258 if (disturbance_vector)
259 x_free(disturbance_vector);
260 // write_png ("tmp_buf.png", 16, (unsigned char *) hf->tmp_buf, hf->max_x, hf->max_y);
261 // write_png ("result_buf.png", 16, (unsigned char *) hf->result_buf, hf->max_x, hf->max_y);
262 // write_png ("hf_buf.png", 16, (unsigned char *) hf->hf_buf, hf->max_x, hf->max_y);
263
264 }
265
hf_threshold(hf_struct_type * hf,hf_type min,hf_type max,gint percent)266 void hf_threshold (hf_struct_type *hf, hf_type min, hf_type max, gint percent) {
267 gint i;
268 if (!hf->tmp_buf)
269 hf_backup(hf);
270 for (i=0; i<hf->max_x*hf->max_y; i++) {
271 *(hf->hf_buf + i) = (hf_type) ((long int) *(hf->tmp_buf +i)* (100-percent) / 100) +
272 ( ( (long int) ( MAX(min,MIN(max,*(hf->tmp_buf + i))) ) * percent) / 100);
273 }
274 }
275
create_gaussian_kernel(gint radius,gdouble base,gdouble exponent,gdouble relative_value_offset,gdouble scale,gboolean not_inverted)276 gdouble *create_gaussian_kernel (gint radius, gdouble base, gdouble exponent, gdouble relative_value_offset, gdouble scale, gboolean not_inverted) {
277 // Create a gaussian kernel, suitable for hf_convolve
278 // To get a "sharpen" or "noisify" kernel, use a negative
279 // relative_value_offset and set not_inverted to FALSE
280 gdouble *kernel, min, max, value_offset, distance, somme;
281 gint x, y, size;
282 size = 1 + 2*radius;
283 if (scale<=0.0)
284 scale=1.0;
285 kernel = (gdouble *) x_malloc(sizeof(gdouble)*size*size, "gdouble (kernel in create_gaussian_kernel)");
286 max = BELLD(base, exponent, (gdouble) 0.0, (gdouble) size*scale);
287 min = BELLD(base, exponent, (gdouble) DIST2(0, 0, radius, radius), (gdouble) size*scale);
288 // printf("Distance: %f\n",(gdouble) DIST2(0, 0, radius, radius));
289 value_offset = relative_value_offset * (max-min);
290 // printf("Radius: %d; Min: %f; Max: %f; Offset: %f\n",radius, min, max, value_offset);
291 somme = 0.0;
292 for (y=0; y<size; y++) {
293 for (x=0; x<size; x++) {
294 distance = DIST2(x,y,radius, radius);
295 if (not_inverted)
296 *(kernel + y*size + x) = value_offset + BELLD(base, exponent, distance, scale*size) - min;
297 else
298 *(kernel + y*size + x) = value_offset + max - BELLD(base, exponent, distance, scale*size);
299 // printf(" %7.4f ",*(kernel + y*size + x) );
300 somme += *(kernel + y*size + x);
301 // if (x==(size-1))
302 // printf("\n");
303 }
304 }
305 // printf("Somme: %7.4f\n", somme);
306 return kernel;
307 }
308
create_hat_kernel(gint radius,gdouble relative_value_offset,gdouble scale,gboolean not_inverted)309 gdouble *create_hat_kernel (gint radius, gdouble relative_value_offset, gdouble scale, gboolean not_inverted) {
310 gdouble *kernel, value_offset, somme, edges_value, sum_edges, center_value;
311 gint x, y, size;
312 size = 1 + 2*radius;
313 if (scale<=0.0)
314 scale=1.0;
315 kernel = (gdouble *) x_malloc(sizeof(gdouble)*size*size, "gdouble (kernel in create_hat_kernel)");
316 edges_value = 1.0;
317 sum_edges = (4.0*(size-2.0) + 4.0) * edges_value;
318 if (not_inverted)
319 center_value = scale;
320 else
321 center_value = -(sum_edges*scale - 1.0) /(gdouble)((size-2) * (size-2));
322 value_offset = relative_value_offset;
323 printf("Radius: %d; Sum_edgest: %f; Center_value: %f; Offset: %f\n",radius, sum_edges, center_value, value_offset);
324 somme = 0.0;
325 for (y=0; y<size; y++) {
326 for (x=0; x<size; x++) {
327 if ((y==0) || (x==0) || (x==(size-1)) || (y==(size-1)))
328 *(kernel + y*size + x) = edges_value;
329 else
330 *(kernel + y*size + x) = center_value;
331 printf(" %7.4f ",*(kernel + y*size + x) );
332 // printf("(%d,%d) = %7.4f ",x,y,*(kernel + y*size + x) );
333 somme += *(kernel + y*size + x);
334 if (x==(size-1))
335 printf("\n");
336 }
337 }
338 printf("Somme: %7.4f\n", somme);
339 return kernel;
340 }
341
create_pixel_kernel(gint radius,gdouble relative_value_offset,gdouble scale,gboolean not_inverted)342 gdouble *create_pixel_kernel (gint radius, gdouble relative_value_offset, gdouble scale, gboolean not_inverted) {
343 gdouble *kernel, value_offset, somme, edges_value, sum_edges, center_value;
344 gint x, y, size;
345 size = 1 + 2*radius;
346 if (scale<=0.0)
347 scale=1.0;
348 kernel = (gdouble *) x_malloc(sizeof(gdouble)*size*size, "gdouble (kernel in create_pixel_kernel)");
349 edges_value = scale * 1.0;
350 sum_edges = ((gdouble) (size*size)-1.0) * edges_value;
351 if (not_inverted)
352 center_value = scale;
353 else
354 center_value = -sum_edges + 1.0;
355 value_offset = relative_value_offset;
356 printf("Radius: %d; Sum_edgest: %f; Center_value: %f; Offset: %f\n",radius, sum_edges, center_value, value_offset);
357 somme = 0.0;
358 for (y=0; y<size; y++) {
359 for (x=0; x<size; x++) {
360 if ((y==x) && (x==radius))
361 *(kernel + y*size + x) = center_value;
362 else
363 *(kernel + y*size + x) = edges_value;
364 printf(" %7.4f ",*(kernel + y*size + x) );
365 // printf("(%d,%d) = %7.4f ",x,y,*(kernel + y*size + x) );
366 somme += *(kernel + y*size + x);
367 if (x==(size-1))
368 printf("\n");
369 }
370 }
371 printf("Somme: %7.4f\n", somme);
372 return kernel;
373 }
374
create_tent_kernel(gint radius,gdouble value_offset,gdouble scale,gboolean not_inverted)375 gdouble *create_tent_kernel (gint radius, gdouble value_offset, gdouble scale, gboolean not_inverted) {
376 gdouble *kernel, somme, dradius, di;
377 gint i, j, x, y, size;
378 size = 1 + 2*radius;
379 if (scale<=0.0)
380 scale=1.0;
381 kernel = (gdouble *) x_malloc(sizeof(gdouble)*size*size, "gdouble (kernel in create_tent_kernel)");
382 printf("Radius: %d; Offset: %f\n",radius, value_offset);
383 somme = 0.0;
384 dradius = (gdouble) radius;
385 for (i=0; i<radius; i++) {
386 for (j=i; j<(size-i); j++) {
387 di = (gdouble) i;
388 if (not_inverted) {
389 *(kernel + i*size + j) = scale * (di+value_offset);
390 *(kernel + (size-i-1)*size + j) = scale * (di+value_offset);
391 }
392 else {
393 *(kernel + i*size + j) = scale * (value_offset-di);
394 *(kernel + (size-i-1)*size + j) = scale * (value_offset-di);
395 }
396 }
397 for (j=(i+1); j<(size-i-1); j++) {
398 if (not_inverted) {
399 *(kernel + j*size + i) = scale * (di+value_offset);
400 *(kernel + j*size + size - i - 1) = scale * (di+value_offset);
401 }
402 else {
403 *(kernel + j*size + i) = scale * (value_offset-di);
404 *(kernel + j*size + size - i - 1) = scale * (value_offset-di);
405 }
406 }
407 }
408 if (not_inverted)
409 *(kernel + radius*size + radius) = scale * (dradius + value_offset);
410 else
411 *(kernel + radius*size + radius) = scale * (value_offset - dradius);
412
413 for (y=0; y<size; y++) {
414 for (x=0; x<size; x++) {
415 printf(" %7.4f ",*(kernel + y*size + x) );
416 // printf("(%d,%d) = %7.4f ",x,y,*(kernel + y*size + x) );
417 somme += *(kernel + y*size + x);
418 if (x==(size-1))
419 printf("\n");
420 }
421 }
422 printf("Somme: %7.4f\n", somme);
423 return kernel;
424 }
425
hf_convolve(hf_struct_type * hf,gint radius,gdouble * matrix,gboolean if_tiles)426 void hf_convolve (hf_struct_type *hf, gint radius, gdouble *matrix, gboolean if_tiles) {
427 // Apply a convolution filter, normalizing the sum of the weights to 1
428 // No normalization takes place when the weights sum is 0 or 1
429 // matrix size = (2*radius+1)*(2*radius+1)
430 gint i,j,k,l, ii, x0,y0,width;
431 gdouble value, weight=0.0, maxd = (gdouble) MAX_HF_VALUE;
432 hf_type *hf_buf;
433 gboolean if_weight = TRUE;
434 if ( (!radius) || (!matrix))
435 return;
436 width = 2*radius + 1;
437 hf_buf = (hf_type *) x_malloc(sizeof(hf_type)*hf->max_x*hf->max_y, "hf_type (hf_buf in hf_convolve)");
438 for (i=0;i<width*width;i++) {
439 weight += *(matrix+i);
440 }
441 if ((weight==0.0) || (weight==1.0))
442 if_weight = FALSE;
443 // printf("W: %s... %f\n",if_weight?"TRUE":"FALSE",weight);
444 for (i=0; i<hf->max_y; i++) {
445 ii = i*hf->max_x;
446 for (j=0; j<hf->max_x; j++) {
447 value = 0.0;
448 y0 = i - radius;
449 x0 = j - radius;
450 // Add the values contained in the (2*radius)+1 *... square centered at current dot
451 for (k = 0; k < width; k++)
452 for (l = 0; l < width; l++) {
453 if ((*matrix+k*width+l)>0.0) {
454 if (if_tiles)
455 value = value + *(matrix+k*width+l) *
456 (gdouble) *(hf->hf_buf +
457 WRAP2(y0+k,hf->max_y)*(hf->max_y) +
458 WRAP2(x0+l,hf->max_x)) ;
459 else
460 value = value + *(matrix+k*width+l) *
461 (gdouble) *(hf->hf_buf +
462 MAX(0,MIN(y0+k,hf->max_y-1))*(hf->max_y) +
463 MAX(0,MIN(x0+l,hf->max_x-1)) ) ;
464 }
465 }
466 if (if_weight)
467 value = value / weight;
468 *(hf_buf + ii+j) = (hf_type) MIN(maxd,MAX(0.0,value)) ;
469 }
470 }
471 memcpy(hf->hf_buf,hf_buf,sizeof(hf_type)*hf->max_x*hf->max_y);
472 x_free(hf_buf);
473 }
474
convolve(gpointer buf,gint max_x,gint max_y,gint data_type,gint radius,gdouble * matrix,gboolean if_tiles)475 void convolve (gpointer buf, gint max_x, gint max_y, gint data_type, gint radius, gdouble *matrix, gboolean if_tiles) {
476 // Apply a convolution filter, normalizing the sum of the weights to 1
477 // No normalization takes place when the weights sum is 0 or 1
478 // matrix size = (2*radius+1)*(2*radius+1)
479 gint i,j,k,l, ii, x0,y0,width;
480 gdouble value, weight=0.0, maxd = (gdouble) MAX_HF_VALUE;
481 gdouble *dbuf, *dbufin;
482 gboolean if_weight = TRUE;
483 if ( (!radius) || (!matrix))
484 return;
485 width = 2*radius + 1;
486 dbufin = (gdouble *) x_malloc(sizeof(gdouble)*max_x*max_y, "gdouble (dbufin in convolve)");
487 dbuf = (gdouble *) x_malloc(sizeof(gdouble)*max_x*max_y, "gdouble (dbuf in convolve)");
488 switch (data_type) {
489 case GINT_ID:
490 for (i=0; i<(max_x*max_y); i++) {
491 *(dbufin+i) = (gdouble) *(((gint*) buf)+i);
492 }
493 break;
494 case HF_TYPE_ID:
495 for (i=0; i<(max_x*max_y); i++) {
496 *(dbufin+i) = (gdouble) *(((hf_type*) buf)+i);
497 }
498 break;
499 case GDOUBLE_ID:
500 for (i=0; i<(max_x*max_y); i++) {
501 *(dbufin+i) = *(((gdouble*) buf)+i);
502 }
503 break;
504 case UNSIGNED_CHAR_ID:
505 for (i=0; i<(max_x*max_y); i++) {
506 *(dbufin+i) = (gdouble) *(((unsigned char*) buf)+i);
507 }
508 break;
509 default:
510 printf("Unexpected data type in convolve\n");
511 exit(0);
512 }
513 for (i=0;i<width*width;i++) {
514 weight += *(matrix+i);
515 }
516 if ((weight==0.0) || (weight==1.0))
517 if_weight = FALSE;
518 // printf("W: %s... %f\n",if_weight?"TRUE":"FALSE",weight);
519 for (i=0; i<max_y; i++) {
520 ii = i*max_x;
521 for (j=0; j<max_x; j++) {
522 value = 0.0;
523 y0 = i - radius;
524 x0 = j - radius;
525 // Add the values contained in the (2*radius)+1 *... square centered at current dot
526 for (k = 0; k < width; k++)
527 for (l = 0; l < width; l++) {
528 if ((*matrix+k*width+l)>0.0) {
529 if (if_tiles)
530 value = value + *(matrix+k*width+l) *
531 (gdouble) *(dbufin +
532 WRAP2(y0+k,max_y)*max_y +
533 WRAP2(x0+l,max_x)) ;
534 else
535 value = value + *(matrix+k*width+l) *
536 (gdouble) *(dbufin +
537 MAX(0,MIN(y0+k,max_y-1))*max_y +
538 MAX(0,MIN(x0+l,max_x-1)) ) ;
539 }
540 }
541 if (if_weight)
542 value = value / weight;
543 *(dbuf + ii+j) = MIN(maxd,MAX(0.0,value)) ;
544 }
545 }
546 switch (data_type) {
547 case GINT_ID:
548 for (i=0; i<(max_x*max_y); i++) {
549 *(((gint*) buf)+i) = (gint) *(dbuf+i);
550 }
551 break;
552 case HF_TYPE_ID:
553 for (i=0; i<(max_x*max_y); i++) {
554 *(((hf_type*) buf)+i) = (hf_type) *(dbuf+i);
555 }
556 break;
557 case GDOUBLE_ID:
558 for (i=0; i<(max_x*max_y); i++) {
559 *(((gdouble*) buf)+i) = *(dbuf+i);
560 }
561 break;
562 case UNSIGNED_CHAR_ID:
563 for (i=0; i<(max_x*max_y); i++) {
564 *(((unsigned char*) buf)+i) = (unsigned char) *(dbuf+i);
565 }
566 break;
567 default:
568 printf("Unexpected data type in convolve\n");
569 exit(0);
570 }
571 x_free(dbuf);
572 x_free(dbufin);
573 }
574
hf_convolveb(hf_struct_type * hf,gint radius,gdouble * matrix,gboolean if_tiles)575 void hf_convolveb (hf_struct_type *hf, gint radius, gdouble *matrix, gboolean if_tiles) {
576 // "Boolean" convolution filter
577 // matrix size = (2*radius+1)*(2*radius+1)
578 // If the current pixel is different from any neighbour, put it black
579 gint i,j,k,l, ii, x0,y0,width;
580 gboolean test;
581 gdouble weight=0.0;
582 hf_type *hf_buf, value;
583 gboolean if_weight = TRUE;
584 width = 2*radius + 1;
585 hf_buf = (hf_type *) x_malloc(sizeof(hf_type)*hf->max_x*hf->max_y, "hf_type (hf_buf in hf_convolveb)");
586 for (i=0;i<width*width;i++) {
587 weight += *(matrix+i);
588 }
589 if ((weight==0.0) || (weight==1.0))
590 if_weight = FALSE;
591 // printf("W: %s... %f\n",if_weight?"TRUE":"FALSE",weight);
592 for (i=0; i<hf->max_y; i++) {
593 ii = i*hf->max_x;
594 value = *(hf->hf_buf + ii+j) ;
595 for (j=0; j<hf->max_x; j++) {
596 test = TRUE;
597 y0 = i - radius;
598 x0 = j - radius;
599 // Add the values contained in the (2*radius)+1 *... square centered at current dot
600 for (k = 0; k < width; k++)
601 for (l = 0; l < width; l++) {
602 if ((*matrix+k*width+l)>0.0) {
603 if (if_tiles)
604 test = test && (value ==
605 *(hf->hf_buf +
606 WRAP2(y0+k,hf->max_y)*(hf->max_y) +
607 WRAP2(x0+l,hf->max_x)) ) ;
608 else
609 test = test && (value ==
610 *(hf->hf_buf +
611 MAX(0,MIN(y0+k,hf->max_y-1))*(hf->max_y) +
612 MAX(0,MIN(x0+l,hf->max_x-1)) ) );
613 }
614 }
615 if (!test)
616 *(hf_buf + ii+j) = 0 ;
617 else
618 *(hf_buf + ii+j) = *(hf->hf_buf+ii+j) ;
619 }
620 }
621 memcpy(hf->hf_buf,hf_buf,sizeof(hf_type)*hf->max_x*hf->max_y);
622 x_free(hf_buf);
623 }
624
625 // static gdouble matrix[] = {0.0625,0.125,0.0625,0.125,0.25,0.125,0.0625,0.125,0.0625}; // Blur
626 // "Sharpen" / add noise; without clamping and with sum = 0: B&W edge detection
627 // static gdouble matrix[] = {-0.5,-0.5,-0.5,-0.5,4.0,-0.5,-0.5,-0.5,-0.5};
628 // static gdouble matrix[] = {1,1,1,1,1,1,1,1,1}; // For min/max (erode/dilate)
629 // static gdouble matrix[] = {-1,-1,-1,-1,8.0,-1,-1,-1,-1};
630
erode_or_dilate(hf_type * hf_in,hf_type * hf_out,gint max_x,gint max_y,gdouble * matrix,gint radius,gint operation,gboolean wrap)631 void erode_or_dilate (hf_type *hf_in, hf_type *hf_out, gint max_x, gint max_y, gdouble *matrix, gint radius, gint operation, gboolean wrap) {
632
633 // Simple erosion with a modified min/max convolution filter
634 // matrix size = (2*radius+1)*(2*radius+1)
635 gint i,j,k,l, ii, x0,y0,width;
636 gdouble value,last,current, weight;
637 width = 2*radius + 1;
638 for (i=0; i<max_y; i++) {
639 ii = i*max_x;
640 for (j=0; j<max_x; j++) {
641 current = (gdouble) *(hf_in + ii+j) ;
642 value = current;
643 y0 = i - radius;
644 x0 = j - radius;
645 // Process the values contained in the (2*radius)+1 *... square centered at current dot
646 for (k = 0; k < width; k++)
647 for (l = 0; l < width; l++) {
648 weight = *(matrix+k*width+l);
649 if (weight>0.0)
650 {
651 // The minimum is normalized with matrix weight (typically distance)
652 if (wrap)
653 last = (gdouble) *(hf_in +
654 WRAP2(y0+k,max_y)*(max_y) +
655 WRAP2(x0+l,max_x) );
656 else
657 last = (gdouble) *(hf_in +
658 MAX(0,MIN(y0+k,max_y-1))*(max_y) +
659 MAX(0,MIN(x0+l,max_x-1)) );
660 if (operation==DILATE) {
661 if (last>current)
662 value = (1.0-weight)*current+weight*last;
663 }
664 else { // Erode
665 if (last<current)
666 value = (1.0-weight)*current+weight*last;
667 }
668 }
669 }
670 *(hf_out + ii+j) = (hf_type) value ;
671 } // j (x) loop
672 } // i (y) loop
673 }
674
dilation(hf_type * hf_in,hf_type * hf_out,gint max_x,gint max_y,gboolean wrap)675 void dilation (hf_type *hf_in, hf_type *hf_out, gint max_x, gint max_y, gboolean wrap) {
676 // Morphological dilation,
677 // done with max instead of boolean OR, for coping with non-binary values
678 static gdouble matrix[] = {0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0};
679 erode_or_dilate (hf_in, hf_out, max_x, max_y, matrix, 1, DILATE, wrap);
680 }
681
erosion(hf_type * hf_in,hf_type * hf_out,gint max_x,gint max_y,gboolean wrap)682 void erosion (hf_type *hf_in, hf_type *hf_out, gint max_x, gint max_y, gboolean wrap) {
683 // Morphological erosion,
684 // done with min instead of boolean AND, for coping with non-binary values
685 static gdouble matrix[] = {0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0};
686 erode_or_dilate (hf_in, hf_out, max_x, max_y, matrix, 1, ERODE, wrap);
687 }
688
min_max_erode(hf_type * hf_in,hf_type * hf_out,gint max_x,gint max_y,gint steps,gboolean wrap)689 void min_max_erode (hf_type *hf_in, hf_type *hf_out, gint max_x, gint max_y, gint steps, gboolean wrap) {
690
691 // Simple erosion with a modified min/max convolution filter
692 // matrix size = (2*radius+1)*(2*radius+1)
693 // static gdouble matrix[] = {0.7,1,0.7,1,1,1,0.7,1,0.7};
694 // static gdouble matrix[] = {0,1,0,1,1,1,0,1,0};
695 // static gdouble matrix[] = {1,1,1,1,1,1,1,1,1};
696 // gint i,j,k,l, ii, s, x0,y0;
697 // gdouble value,last,current, weight;
698 gint s, width, radius = 1;
699 width = 2*radius + 1;
700 hf_type *in,*out;
701 in = hf_in;
702 out = hf_out;
703 for (s=1; s<steps; s++) {
704 printf("STEP %d\n",s);
705 if (s>1)
706 swap_buffers ((gpointer*) &hf_in, (gpointer*) &hf_out);
707 erosion (hf_in, hf_out, max_x, max_y, wrap);
708 }
709 if (out==hf_in)
710 memcpy(out,hf_out,max_x*max_y*sizeof(hf_type));
711 }
712
min_max_spread(hf_type * hf_in,hf_type * hf_out,gint max_x,gint max_y,gint steps,gboolean wrap)713 void min_max_spread (hf_type *hf_in, hf_type *hf_out, gint max_x, gint max_y, gint steps, gboolean wrap) {
714
715 // Simple erosion (dark areas spread) derivated from min_max_erode
716 // matrix size = 2*2 of 1.0
717 // Use for expanding dark areas (lines) 1 pixel at a time, towards east and south
718 gint i,j, ii, s, x0,y0;
719 gdouble value,last,current;
720 hf_type *in,*out;
721 in = hf_in;
722 out = hf_out;
723 for (s=1; s<steps; s++) {
724 if (s>1)
725 swap_buffers ((gpointer*) &hf_in, (gpointer*) &hf_out);
726 for (i=0; i<max_y; i++) {
727 ii = i*max_x;
728 for (j=0; j<max_x; j++) {
729 current = (gdouble) *(hf_in + ii+j) ;
730 value = current;
731 if (wrap) {
732 y0 = WRAP2(i - 1,max_y);
733 x0 = WRAP2(j - 1,max_y);
734 }
735 else {
736 y0 = MAX(0,i - 1);
737 x0 = MAX(0,j - 1);
738 }
739 last = (gdouble) *(hf_in + VECTORIZE(x0,y0,max_x) );
740 if (last<current)
741 value = last;
742 y0 = i;
743 last = (gdouble) *(hf_in + VECTORIZE(x0,y0,max_x) );
744 if (last<current)
745 value = last;
746 if (wrap) {
747 y0 = WRAP2(i - 1,max_y);
748 }
749 else {
750 y0 = MAX(0,i - 1);
751 }
752 x0 = j;
753 last = (gdouble) *(hf_in + VECTORIZE(x0,y0,max_x) );
754 if (last<current)
755 value = last;
756 *(hf_out + ii+j) = (hf_type) value ;
757 } // j (x) loop
758 } // i (y) loop
759 } // s (steps) loop
760 if (out==hf_in)
761 memcpy(out,hf_out,max_x*max_y*sizeof(hf_type));
762 }
763
hf_erode(hf_struct_type * hf,gint steps)764 void hf_erode (hf_struct_type *hf, gint steps) {
765 hf_type *hf_in;
766 gint totsize;
767 totsize = sizeof(hf_type) * hf->max_x *hf->max_y;
768 hf_in = (hf_type *) x_malloc(totsize, "hf_type (hf_in in hf_erode)");
769 memcpy (hf_in,hf->hf_buf, totsize);
770 if (hf->if_tiles)
771 min_max_erode (hf_in, hf->hf_buf, hf->max_x, hf->max_y, steps, *hf->if_tiles);
772 else
773 min_max_erode (hf_in, hf->hf_buf, hf->max_x, hf->max_y, steps, TRUE);
774 x_free(hf_in);
775 }
776
init_blocks(gpointer buf,gint max_x,gint max_y,gint radius,gint nblocks,gfloat phase,gboolean if_tiles,gint data_type)777 gpointer init_blocks (gpointer buf, gint max_x, gint max_y, gint radius, gint nblocks,
778 gfloat phase, gboolean if_tiles, gint data_type) {
779 // Average the pixels of hf into a smaller nblocks * nblocks hf_buf
780 // Kind of "undersampling"; basis for processes as "smooth" and "cityscape"
781 gint i, j, k, l, x, y, size;
782 gpointer blocks;
783 gdouble value;
784 switch (data_type) {
785 case GINT_ID:
786 size = sizeof(gint);
787 break;
788 case HF_TYPE_ID:
789 size = sizeof(hf_type);
790 break;
791 case GDOUBLE_ID:
792 size = sizeof(gdouble);
793 break;
794 case UNSIGNED_CHAR_ID:
795 size = sizeof(unsigned char);
796 break;
797 default:
798 printf(_("Unexpected option in %s; contact the programmer!"),
799 "init_blocks");
800 printf("\n");
801 }
802 blocks = (gpointer) x_malloc(nblocks * nblocks * size, "gpointer (blocks in init_blocks)");
803 for (i=0; i<nblocks; i++) {
804 for (j=0; j<nblocks; j++) {
805 value = 0;
806 for (k=((i+phase)*radius); k<((i+1+phase)*radius); k++)
807 for (l=((j+phase)*radius); l<((j+1+phase)*radius); l++) {
808 if (if_tiles) {
809 x = (max_x+l)%max_x;
810 y = (max_y+k)%max_y;
811 }
812 else {
813 x = MIN(l, max_x-1);
814 y = MIN(k, max_y-1);
815 }
816 switch (data_type) {
817 case GINT_ID:
818 value = value + (gdouble) *( ((gint *) buf) + y*max_x + x);
819 break;
820 case HF_TYPE_ID:
821 value = value + (gdouble) *( ((hf_type *) buf) + y*max_x + x);
822 break;
823 case GDOUBLE_ID:
824 value = value + *( ((gdouble *) buf) + y*max_x + x);
825 break;
826 case UNSIGNED_CHAR_ID:
827 value = value + (gdouble) *( ((unsigned char *) buf) + y*max_x + x);
828 break;
829 default:
830 printf(_("Unexpected option in %s; contact the programmer!"),
831 "init_blocks");
832 printf("\n");
833 }
834 } // End k and l loops
835 value = value / (gdouble) (radius * radius);
836 switch (data_type) {
837 case GINT_ID:
838 *( ((gint *) blocks) + (i*nblocks)+j) = (gint) value;
839 break;
840 case HF_TYPE_ID:
841 *( ((hf_type *) blocks) + (i*nblocks)+j) = (hf_type) value;
842 break;
843 case GDOUBLE_ID:
844 *( ((gdouble *) blocks) + (i*nblocks)+j) = value;
845 break;
846 case UNSIGNED_CHAR_ID:
847 *( ((unsigned char *) blocks) + (i*nblocks)+j) = (unsigned char) value;
848 break;
849 default:
850 printf(_("Unexpected option in %s; contact the programmer!"),
851 "init_blocks");
852 printf("\n");
853 }
854 } // End j loop
855 } // End i loop
856 return blocks;
857 }
858
smooth_buf(gpointer buf,gint max_x,gint max_y,gint radius,dist_matrix_struct * dm,gboolean wrap,gdouble ** gauss_list,gint data_type)859 void smooth_buf (gpointer buf, gint max_x, gint max_y, gint radius, dist_matrix_struct *dm, gboolean wrap, gdouble **gauss_list, gint data_type) {
860 // Smoothing process, based on a sampling, followed by a rebuild
861 // based on a normal (gaussian) 3d function (used as "shape filters" elsewhere in Geomorph)
862 // For radius <= 36, we use a "separated" convolution in the spacial domain
863 //
864 // dist_matrix = the usual precompiled distances matrix, for managing gaussian filters
865 // wrap = TRUE => we reuse pixel 0 after pixel max_x or max_y
866
867 gint i, j, k, l, nblocks, box, index, indexf, offset, half_radius;
868 gpointer blocks;
869 filter_struct *filter;
870 gdouble *hf_weight,*hf_values, value;
871 gdouble *v;
872
873 if (!radius) // We are not supposed to be here!
874 return;
875
876 if (radius<=36) {
877
878 if (!gauss_list[radius])
879 gauss_list[radius] = normalized_bell_new(radius);
880 v = gauss_list[radius];
881
882 convolve_normalized_vector (buf,
883 buf,
884 max_x,
885 max_y,
886 wrap,
887 radius,
888 v,
889 data_type);
890 }
891 else {
892
893 // We add repeatedly a gaussian bell, at each sampled pixel
894 // 2003-02-14: finally there is only 1 level left, since a separable convolution is used
895 // Finally the best rebuilding box radius is equal to the specified radius
896 box = radius;
897 filter = sharp_filter_new(box<<1,dm); // The size must be even
898 half_radius = radius>>1;
899
900 // 1- Sample
901 nblocks = (gint) ((max_x-1+half_radius) / half_radius);
902 // printf("NBLOCKS: %d\n",nblocks);
903 blocks = init_blocks(buf, max_x, max_y, half_radius, nblocks, 0, wrap, data_type);
904
905 hf_values = (gdouble *) x_calloc(max_x*max_y, sizeof(gdouble), "gdouble (hf_values in smooth_buf)");
906 hf_weight = (gdouble *) x_calloc(max_x*max_y, sizeof(gdouble), "gdouble (hf_weight in smooth_buf)");
907
908 // 2- Rebuild
909 // The sampling / rebuild process introduces a translation, which we must correct
910 offset = half_radius/2;
911 for (i=0; i<nblocks; i++) {
912 for (j=0; j<nblocks; j++) {
913
914 for (k=0; k<filter->hf_size; k++)
915 for (l=0; l<filter->hf_size; l++) {
916
917 if (wrap)
918 index = WRAP2((i*half_radius+k-box+offset),max_y)*max_x
919 + WRAP2(j*half_radius +l-box+offset, max_x);
920 else
921 index = MAX(0,MIN(i*half_radius+k-box+offset,max_y-1))*max_x
922 + MAX(0,MIN(j*half_radius +l-box+offset, max_x-1)) ;
923 indexf = VECTORIZE(CENTER_XY(k,filter->hf_size),
924 CENTER_XY(l,filter->hf_size),filter->hf_size>>1) ;
925 // printf("INDEX dans filter->value: %d\n", indexf );
926 switch (data_type) {
927 case GINT_ID:
928 value = (gdouble) *(((gint *)blocks) + i*nblocks+j);
929 break;
930 case HF_TYPE_ID:
931 value = (gdouble) *(((hf_type *)blocks) + i*nblocks+j);
932 break;
933 case GDOUBLE_ID:
934 value = *(((gdouble *)blocks) + i*nblocks+j);
935 break;
936 case UNSIGNED_CHAR_ID:
937 value = (gdouble) *(((unsigned char *) blocks) + i*nblocks+j);
938 break;
939 default:
940 printf(_("Unexpected option in %s; contact the programmer!"),
941 "smooth_buf");
942 printf("\n");
943 }
944 *(hf_values + index) = *(hf_values + index) + value * *(filter->values + indexf) ;
945 *(hf_weight+index) = *(hf_weight+index) + *(filter->values+ indexf);
946 // printf("(i,j,k,l): (%d,%d,%d,%d); INDEX: %d; values: %f; w: %f\n",i,j,k,l,index,*(hf_values + index) ,*(hf_weight+index) );
947 }
948 }
949 }
950
951 for (i=0; i<max_x*max_y; i++) {
952
953 // Modified 2006-01-27 for different data types
954 switch (data_type) {
955 case GINT_ID:
956 *( ((gint *) buf) + i) = (gint) ( (*(hf_values+i))/(*(hf_weight+i)) ) ;
957 break;
958 case HF_TYPE_ID:
959 *( ((hf_type *) buf) + i) = (hf_type) ( (*(hf_values+i))/(*(hf_weight+i)) ) ;
960 break;
961 case GDOUBLE_ID:
962 *( ((gdouble *) buf) + i) = (gdouble) ( (*(hf_values+i))/(*(hf_weight+i)) ) ;
963 break;
964 case UNSIGNED_CHAR_ID:
965 *( ((unsigned char *) buf) + i) = (unsigned char) ( (*(hf_values+i))/(*(hf_weight+i)) ) ;
966 break;
967 default:
968 printf(_("Unexpected option in %s; contact the programmer!"),
969 "smooth_buf");
970 printf("\n");
971 }
972 }
973 x_free(blocks);
974 x_free(hf_values);
975 x_free(hf_weight);
976 filter_free(filter);
977 // We reapply the process for smoothing the artifacts of the rebuilt HF
978 smooth_buf(buf, max_x, max_y, 36, dm, wrap, gauss_list, data_type);
979 } // End else radius <= 36
980 }
981
hf_smooth(hf_struct_type * hf,gint radius,dist_matrix_struct * dm,gboolean wrap,gdouble ** gauss_list)982 void hf_smooth (hf_struct_type *hf, gint radius,
983 dist_matrix_struct *dm, gboolean wrap,gdouble **gauss_list) {
984 smooth_buf(hf->hf_buf, hf->max_x, hf->max_y, radius, dm, wrap, gauss_list, HF_TYPE_ID);
985 };
hf_mosaic(hf_struct_type * hf,gint radius,gint nblocks,gfloat phase)986 void hf_mosaic (hf_struct_type *hf, gint radius, gint nblocks, gfloat phase) {
987 gint i,j;
988 hf_type *blocks;
989 // Build the mosaic
990 blocks = (hf_type *) init_blocks(hf->tmp_buf, hf->max_x, hf->max_y, radius, nblocks, phase, FALSE, HF_TYPE_ID);
991 // Transfer the mosaic into the hf buffer
992 for (i=0; i<hf->max_y; i++)
993 for (j=0; j<hf->max_x; j++)
994 *(hf->hf_buf+WRAP(i+(gint) (radius*phase),hf->max_x)*hf->max_x+j+
995 WRAP((gint) (radius*phase),hf->max_y)) =
996 *(blocks+((gint) (i/radius) )*nblocks + (gint) (j/radius));
997 }
998
hf_honeycomb(hf_type * hf_in,hf_type * hf_out,gint max,gint radius,gint border)999 void hf_honeycomb (hf_type *hf_in, hf_type *hf_out, gint max, gint radius, gint border) {
1000 // Quantize with an hexagonal structure
1001 // Distance between cells centers = 2*radius
1002 // "border" = with of the line separating each hexagon, in pixels
1003 // This transformation does not preserve tiling
1004 // Important: radius must be even, otherwise we get black pixels (rounding problem)
1005 gint x,y,i,j, nbcellsh, nbcellsw, h, w, index, nrow, npix;
1006 gint baseh, basew, current_cell, block_h, xoffset, xstep;
1007 static gfloat tan60=1.73205, sin60=0.86602;
1008 unsigned long int sum;
1009 hf_type *cell_values;
1010 // Hexagons are not symmetrical in a square world...
1011 // I chose to put the longer axis on the width
1012 // h = height of each cell
1013 // w = width
1014 /* _ _
1015 / \_/ \
1016 \_/ \_/
1017 */
1018 radius = 2* (gint) ceil((gdouble) radius/2.0);
1019 h = (gint) ((gfloat) radius * sin60);
1020 w = 3 * radius;
1021 nbcellsh = (gint) ceil ( ((gdouble) 2*radius + max) / (gdouble) h);
1022 nbcellsw = (gint) ceil (((gdouble) 2*radius + max) / (gdouble) w);
1023 cell_values = (hf_type *) x_calloc(1+nbcellsw*nbcellsh, sizeof(hf_type), "hf_type (cell_values in hf_honeycomb)");
1024 // cell_values contains the average value for the current cell
1025 // We use hf->hf_buf to store the index of the cell to which the related pixel belongs
1026 for (i=0; i<max*max; i++)
1027 *(hf_out+i) = 0;
1028 // Processing the neighbourhood of each cell...
1029 // We "center" the structure
1030 baseh = (max%h)/2;
1031 // basew = (max%radius)/2;
1032 basew = 0;
1033 current_cell = 1;
1034 // Cell 0 value is 0, so we can detect glitches...
1035 block_h = radius*sin60;
1036 nrow=0;
1037 // printf("RADIUS: %d; CH: %d; CW: %d; WxH: %d; BASEH: %d; BASEW: %d; BLOCK_h: %d\n",radius, nbcellsh, nbcellsw, nbcellsh*nbcellsw, baseh, basew, block_h);
1038 for (y=baseh-h/2; y<max+h; y+=h) {
1039 xstep = 3*radius ;
1040 // The x offset is not the same depending if the row is even or odd
1041 if (nrow%2) {
1042 xoffset = basew-radius/2; // This is why radius must be even
1043 }
1044 else {
1045 xoffset = basew+radius;
1046 }
1047 // printf("ROW: %d; XOFFSET: %d; XSTEP: %d; CELL: %d\n",nrow,xoffset,xstep, current_cell);
1048 nrow++;
1049 for (x=xoffset; x<max+radius; x+=xstep) {
1050 // printf("X: %d; ",x);
1051 sum = 0;
1052 npix = 0;
1053 for (i=-block_h+border; i<block_h-border; i++) {
1054 if ((y+i)<0)
1055 continue;
1056 if ((y+i)>=max)
1057 break;
1058 for (j=-radius+border; j<radius-border; j++) {
1059 index = x + j;
1060 if (index<0)
1061 continue;
1062 if (index>=max)
1063 break;
1064 index = VECTORIZE(index,y+i,max);
1065 // Test: is it inside the hexagon?
1066 // If abs(x) (x is j) <=cos60*radius (<= radius/2), it's certainly OK
1067 // otherwise, we must check if abs(i) / (radius-abs(x)) < tan60
1068 if ( (ABS(j)<= (radius-border)/2) ||
1069 (( ((gfloat) ABS(i)) / (gfloat) (radius-border-ABS(j)) ) < tan60 ))
1070 {
1071 *(hf_out+index) = current_cell;
1072 npix++;
1073 sum = sum + *(hf_in+index);
1074 }
1075 }
1076 }
1077 // printf("Current_CELL: %d; NPIX: %d; SUM: %d\n",current_cell,npix,sum);
1078 if (npix) // Some cells are outside the image, for some radius
1079 *(cell_values+current_cell) = sum / npix;
1080 current_cell++;
1081 }
1082 }
1083 for (i=0; i<max*max; i++)
1084 if ((nbcellsh*nbcellsw)<*(hf_out+i)) // "Debugging" test - it's not supposed to happen...
1085 *(hf_out+i) = 0xFFFF;
1086 else
1087 *(hf_out+i) = *(cell_values+*(hf_out+i));
1088 x_free(cell_values);
1089 }
1090
hf_hexagon(hf_struct_type * hf,gint radius,gint border,gint smooth_radius,dist_matrix_struct * dist_matrix,gboolean apply_postprocess,gdouble ** gauss_list)1091 void hf_hexagon (hf_struct_type *hf, gint radius, gint border,
1092 gint smooth_radius,
1093 dist_matrix_struct *dist_matrix,
1094 gboolean apply_postprocess,
1095 gdouble **gauss_list) {
1096
1097 if (radius<4)
1098 return;
1099
1100 if (!hf->tmp_buf)
1101 hf_backup(hf);
1102 // We need a second temporary buffer for storing
1103 // the result of the quantization process before the postprocessing
1104 if (!hf->result_buf)
1105 hf->result_buf = (hf_type *) x_malloc(hf->max_x*hf->max_y*sizeof(hf_type), "hf_type (result_buf in hf_hexagon)");
1106
1107 hf_honeycomb (hf->tmp_buf, hf->result_buf, hf->max_x, radius, border);
1108
1109 if ( apply_postprocess && smooth_radius ) {
1110 // Smooth works on hf_buf, so we must swap it with result_buf
1111 swap_buffers ((gpointer*) &hf->result_buf, (gpointer*) &hf->hf_buf);
1112 if (smooth_radius)
1113 hf_smooth(hf, smooth_radius, dist_matrix, FALSE, gauss_list);
1114 swap_buffers ((gpointer*) &hf->result_buf, (gpointer*) &hf->hf_buf);
1115 } // end of "if (apply_postprocess)"
1116
1117 }
1118
1119
hf_city(hf_struct_type * hf,gint radius,gint skyscraper_var,gint streets_width)1120 void hf_city(hf_struct_type *hf, gint radius, gint skyscraper_var, gint streets_width) {
1121 // radius: skyscraper width, in pixels
1122 // skyscaper_var: upper limit of random width variation, in % of radius
1123 // streets_width: width of black lines between blocks simulating streets, in % of radius
1124 gint i,j,k,l,nblocks,str_w, var;
1125 hf_type *blocks;
1126 gint *rand_y, *rand_x;
1127
1128 if (!radius)
1129 return;
1130
1131 if (!hf->tmp_buf)
1132 hf_backup(hf);
1133 if (!hf->result_buf)
1134 hf->result_buf = (hf_type *) x_malloc(hf->max_x*hf->max_y*sizeof(hf_type), "hf_type (result_buf in hf_city)");
1135
1136 nblocks = (gint) ((hf->max_x-1+radius) / radius);
1137 // Sample the height of skyscapers
1138 blocks = (hf_type *) init_blocks(hf->tmp_buf, hf->max_x, hf->max_y, radius, nblocks, 0, FALSE, HF_TYPE_ID);
1139 // Transfer the "mosaic" into the hf buffer
1140 // Apply random width variation, bounded by skyscraper_var %, and draw streets
1141 str_w = (streets_width*radius)/200; // Divide by 200: a street extends one half on each block
1142 rand_x = (gint *) x_malloc((nblocks+2)*sizeof(gint), "gint (rand_x in hf_city)");
1143 rand_y = (gint *) x_malloc((nblocks+2)*sizeof(gint), "gint (rand_y in hf_city)");
1144 var = (skyscraper_var*radius)/100;
1145 *(rand_x) = 0;
1146 *(rand_y) = 0;
1147 if (var) {
1148 for (i=1; i<(nblocks+1); i++) {
1149 *(rand_y+i) = ((i-1)*radius) + ((rand()%var) - (var/2));
1150 *(rand_x+i) = ((i-1)*radius) + ((rand()%var) - (var/2));
1151 // printf("RAND_X(%d): %d; RAND_Y(%d): %d\n",i, *(rand_x+i),i,*(rand_y+i));
1152 }
1153 }
1154 else
1155 for (i=1; i<(nblocks+1); i++) {
1156 *(rand_y+i) = (i-1)*radius;
1157 *(rand_x+i) = (i-1)*radius;
1158 }
1159 *(rand_y+nblocks+1) = hf->max_y;
1160 *(rand_x+nblocks+1) = hf->max_x;
1161 for (i=0; i<(nblocks+1); i++) {
1162 for (j=0; j<(nblocks+1); j++) {
1163 for (k= *(rand_y+i); (k<*(rand_y+i+1)) ; k++)
1164 for (l= *(rand_x+j); (l<*(rand_x+j+1)) ; l++)
1165 if ((k<(*(rand_y+i) + str_w)) || (k>*(rand_y+i+1)-str_w) ||
1166 (l<(*(rand_x+j) + str_w)) || (l>*(rand_x+j+1)-str_w) )
1167 *(hf->result_buf+WRAP(k,hf->max_y)*hf->max_x+
1168 WRAP(l,hf->max_x)) = 0;
1169 else
1170 *(hf->result_buf+WRAP(k,hf->max_y)*hf->max_x+
1171 WRAP(l,hf->max_x)) =
1172 *(blocks+WRAP(i-1,nblocks)*nblocks + WRAP(j-1,nblocks)) ;
1173 }
1174 }
1175 if (rand_y)
1176 x_free(rand_y);
1177 if (rand_x)
1178 x_free(rand_x);
1179 }
1180
hf_stretch_v(hf_struct_type * hf)1181 void hf_stretch_v (hf_struct_type *hf) {
1182 // Stretch vertically (compress horizontally)
1183 // We simply divide the scale by 2 horizontally
1184 gint x,y;
1185 if (!hf->tmp_buf)
1186 hf_backup(hf);
1187 for (y=0; y<hf->max_y; y++) {
1188 for (x=0; x<hf->max_x; x++) {
1189 *(hf->hf_buf + VECTORIZE(x,y,hf->max_x)) =
1190 *(hf->tmp_buf + VECTORIZE(WRAP(x<<1,hf->max_x),y,hf->max_x));
1191 }
1192 }
1193 }
hf_stretch_h(hf_struct_type * hf)1194 void hf_stretch_h (hf_struct_type *hf) {
1195 // Stretch horizontally (compress vertically)
1196 // We simply divide the scale by 2 vertically
1197 gint x,y;
1198 if (!hf->tmp_buf)
1199 hf_backup(hf);
1200 for (y=0; y<hf->max_y; y++) {
1201 for (x=0; x<hf->max_x; x++) {
1202 *(hf->hf_buf + VECTORIZE(x,y,hf->max_x)) =
1203 *(hf->tmp_buf + VECTORIZE(x,WRAP((y<<1),hf->max_y),hf->max_x));
1204 }
1205 }
1206 }
1207
hf_math_fn(hf_struct_type * hf,gint math_op,gdouble parameter)1208 void hf_math_fn (hf_struct_type *hf, gint math_op, gdouble parameter) {
1209 // Some mathematical pixel transformations
1210 gint i;
1211 static gdouble PI2 = 6.283185307, maxd = (gdouble) 0xFFFF;
1212 gdouble *doublebuf, ratiolog;
1213 doublebuf = xalloc_double_hf(hf->max_x, "gdouble (doublebuf in hf_math_fn)");
1214 ratiolog = (10.0-parameter) * log1p(maxd) / maxd;
1215 // printf("MAXD: %f; LOG(MAX): %f; RATIOLOG: %f\n",maxd, log1p(maxd), ratiolog);
1216 if (!hf->tmp_buf)
1217 hf_backup(hf);
1218 // printf("MATH_OP: %d; PARAM: %f\n",math_op, parameter);
1219
1220 // Patch 2010-11-21 to avoid black result
1221 if (parameter==1.0)
1222 parameter += 0.001;
1223
1224 for (i=0; i<hf->max_x*hf->max_y; i++) {
1225
1226 switch (math_op) {
1227 case POWER_OP:
1228 *(doublebuf +i) = pow( (gdouble) *(hf->tmp_buf+i), parameter);
1229 break;
1230 case BASE_OP:
1231 *(doublebuf +i) = pow( parameter, ((gdouble) *(hf->tmp_buf+i)) / maxd );
1232 break;
1233 case LOG_OP:
1234 *(doublebuf+i) = ratiolog * *(hf->tmp_buf+i) + parameter * log1p((gdouble)*(hf->tmp_buf+i));
1235 break;
1236 case SINE_OP:
1237 *(doublebuf+i) = sin( parameter * PI2 *((gdouble) *(hf->tmp_buf+i)) / maxd );
1238 default:
1239 *(doublebuf+i) = (gdouble) *(hf->tmp_buf+i);
1240 }
1241 // if (!((i+(hf->max_x>>1))%hf->max_x))
1242 // printf("@i= %d; tmp_buf: %d; doublebuf: %f; \n", i,*(hf->tmp_buf+i),*(doublebuf+i) );
1243 }
1244 hf_min_max (hf);
1245 hf_double_clamp (hf, hf->min, hf->max, doublebuf);
1246 x_free(doublebuf);
1247 }
1248
hf_black_n_white(hf_type * buf,gint length,hf_type pivot)1249 void hf_black_n_white (hf_type *buf, gint length, hf_type pivot) {
1250 // Stretch the contrast of buf to 0x0000 and 0xFFFF,
1251 // using pivot as the value under which all is black
1252 // and over which all is white
1253 gint i;
1254 for (i=0; i<length; i++) {
1255 if (*(buf+i)<=pivot)
1256 *(buf+i) = 0;
1257 else
1258 *(buf+i) = 0xFFFF;
1259 }
1260 }
1261
merge_lift_edges(hf_struct_type * hf,gint level)1262 void merge_lift_edges (hf_struct_type *hf, gint level) {
1263 gint i;
1264 hf_type *hf1,*hf2,*out;
1265 gdouble ratio, glevel;
1266 glong abs_level = ((LIFT_EDGES_MAX>>1) * level) / 100;
1267 glevel = ((gdouble) level) / 100.0;
1268 hf1 = hf->tmp_buf;
1269 hf2 = hf->result_buf;
1270 out = hf->hf_buf;
1271 if ((hf->min>abs_level) || ((MAX_HF_VALUE-hf->max)>abs_level)) {
1272 ratio = ((gdouble) (hf->max - hf->min)) / ((gdouble) (hf->max - hf->min + 2*abs_level));
1273 for (i=0; i<(hf->max_x*hf->max_y); i++)
1274 *(out+i) = (hf_type) (ratio * (gdouble) ((glong) *(hf1+i)) + (glong) (glevel * (gdouble) *(hf2+i)) );
1275 }
1276 else
1277 for (i=0; i<(hf->max_x*hf->max_y); i++)
1278 *(out+i) = (hf_type) ((glong) *(hf1+i)) + (glong) (glevel * (gdouble) *(hf2+i));
1279 }
1280
lift_edges(hf_type * in,hf_type * out,gint max_x,gint max_y,gint radius,dist_matrix_struct * dm,gboolean wrap,gdouble ** gauss_list,gboolean use_black_point)1281 void lift_edges (hf_type *in, hf_type *out, gint max_x, gint max_y, gint radius,
1282 dist_matrix_struct *dm, gboolean wrap,gdouble **gauss_list, gboolean use_black_point) {
1283 // Writes in *out the differences to add to *in to lift the edges
1284 // of areas with black borders
1285 //
1286 // 1. Threshold the grid
1287 // (all pixels under the "black point" are black, others are white)
1288 // 2006-12-25: threshold only if use_black_point is TRUE
1289 // TRUE: better for surfaces with dark lines (e.g. hexagons)
1290 // FALSE: better for other quantized surfaces (e.g. terraces)
1291 // 2. Smooth
1292 // 3. Subtract the smooth result from the lower contrast source
1293 //
1294
1295 hf_type *tmp, min, max;
1296 gint i;
1297 tmp = (hf_type *) x_malloc(sizeof(hf_type)*max_x*max_y, "hf_type (tmp in lift_edges)");
1298 if ((!radius) || (!in) || (!out))
1299 return;
1300 min = max = *in;
1301 for (i=1; i<(max_x*max_y); i++)
1302 if (min>*(in+i))
1303 min = *(in+i);
1304 else
1305 if (max<*(in+i))
1306 max = *(in+i);
1307 if (use_black_point)
1308 hf_black_n_white (in, max_x*max_y, (max-min)>>3);
1309 memcpy(tmp,in,sizeof(hf_type)*max_x*max_y);
1310 smooth_buf (in, max_x, max_y, radius, dm, wrap, gauss_list, HF_TYPE_ID);
1311 hf_subtract (tmp, in, out, max_x*max_y, OVERFLOW_ZERO);
1312 hf_clamp_buffer (out, max_x * max_y, 0, LIFT_EDGES_MAX);
1313 x_free(tmp);
1314 }
1315
build_noise_map(gdouble wave_length_in_percent,gint max_x,gint max_y,dist_matrix_struct * dm,gdouble ** gauss_list)1316 hf_type *build_noise_map (gdouble wave_length_in_percent, gint max_x, gint max_y, dist_matrix_struct *dm, gdouble **gauss_list) {
1317 hf_type *output;
1318 gint nbpix, i;
1319 output = (hf_type *) x_calloc(sizeof(hf_type), max_x*max_y, "hf_type (output in build_noise_map)");
1320 nbpix = (gint) (100.0 / wave_length_in_percent);
1321 nbpix *= nbpix;
1322 nbpix *= 4;
1323 for (i=0; i<nbpix; i++) {
1324 *(output + (rand()%max_x) + (rand()%max_y)*max_x) = MAX_HF_VALUE;
1325 }
1326 smooth_buf (output, max_x, max_y, (gint) ((wave_length_in_percent / 100.0) * (gdouble) max_x), dm, TRUE, gauss_list, HF_TYPE_ID);
1327 hf_clamp_buffer(output, max_x*max_y, 0, MAX_HF_VALUE);
1328 // smooth_hf_buf (output, max_x, max_y, (gint) ((wave_length_in_percent / 100.0) * (gdouble) max_x), dm, TRUE, gauss_list);
1329 return output;
1330 }
1331
dfind_maximum(gdouble * in,hf_type * output,gint max_x,gint max_y,hf_type foreground_value,hf_type background_value)1332 void dfind_maximum (gdouble *in, hf_type *output, gint max_x, gint max_y, hf_type foreground_value, hf_type background_value) {
1333 // Find the maxima of the "in" map
1334 // Fill the maxima with a foreground value,
1335 // fill the other pixels with a background value
1336 // in the "output" result
1337 // "output" is already allocated
1338
1339 // If both foreground and background values are 0, we use 0 as the background
1340 // and the "in" value as the foreground
1341
1342 gdouble current;
1343 gint x,y,xx,yy,xxx,yyy;
1344 gboolean test;
1345
1346 test = (!foreground_value) && !background_value;
1347
1348 for (y=0; y<max_y; y++)
1349 for (x=0; x<max_x; x++) {
1350
1351 // We test the 4 axis to determine if the current pixel is
1352 // a local maximum
1353
1354 current = *(in + x + y*max_x);
1355
1356 if (test)
1357 foreground_value = (hf_type) current;
1358
1359 xx = WRAP2(x-1,max_x);
1360 yy = WRAP2(y-1,max_y);
1361 xxx = WRAP2(x+1,max_x);
1362 yyy = WRAP2(y+1,max_y);
1363
1364 if ((current>*(in+xxx+yyy*max_x)) &&
1365 (current>*(in+xx+yy*max_x))) {
1366 *(output+x+y*max_x) = foreground_value;
1367 continue;
1368 }
1369
1370 xx = x;
1371 yy = WRAP2(y-1,max_y);
1372 xxx = x;
1373 yyy = WRAP2(y+1,max_y);
1374
1375 if ((current>*(in+xxx+yyy*max_x)) &&
1376 (current>*(in+xx+yy*max_x))) {
1377 *(output+x+y*max_x) = foreground_value;
1378 continue;
1379 }
1380
1381 xx = WRAP2(x+1,max_x);
1382 yy = WRAP2(y-1,max_y);
1383 xxx = WRAP2(x-1,max_x);
1384 yyy = WRAP2(y+1,max_y);
1385
1386 if ((current>*(in+xxx+yyy*max_x)) &&
1387 (current>*(in+xx+yy*max_x))) {
1388 *(output+x+y*max_x) = foreground_value;
1389 continue;
1390 }
1391
1392 xx = WRAP2(x+1,max_x);
1393 yy = y;
1394 xxx = WRAP2(x-1,max_x);
1395 yyy = y;
1396
1397 if ((current>*(in+xxx+yyy*max_x)) &&
1398 (current>*(in+xx+yy*max_x))) {
1399 *(output+x+y*max_x) = foreground_value;
1400 continue;
1401 }
1402
1403 *(output+x+y*max_x) = background_value;
1404 }
1405 }
1406
find_maximum(hf_type * in,hf_type * output,gint max_x,gint max_y,hf_type foreground_value,hf_type background_value)1407 void find_maximum (hf_type *in, hf_type *output, gint max_x, gint max_y, hf_type foreground_value, hf_type background_value) {
1408 // Find the maxima of the "in" map
1409 // Fill the maxima with a foreground value, the other pixels a background value
1410 // in the "output" result
1411 // "output" is already allocated
1412
1413 // If both foreground and background values are 0, we use 0 as the background
1414 // and the "in" value as the foreground
1415
1416 hf_type current;
1417 gint x,y,xx,yy,xxx,yyy;
1418 gboolean test;
1419
1420 test = (!foreground_value) && !background_value;
1421
1422 for (y=0; y<max_y; y++)
1423 for (x=0; x<max_x; x++) {
1424
1425 // We test the 4 axis to determine if the current pixel is
1426 // a local maximum
1427
1428 current = *(in + x + y*max_x);
1429
1430 if (test)
1431 foreground_value = (hf_type) current;
1432
1433 xx = WRAP2(x-1,max_x);
1434 yy = WRAP2(y-1,max_y);
1435 xxx = WRAP2(x+1,max_x);
1436 yyy = WRAP2(y+1,max_y);
1437
1438 if ((current>*(in+xxx+yyy*max_x)) &&
1439 (current>*(in+xx+yy*max_x))) {
1440 *(output+x+y*max_x) = foreground_value;
1441 continue;
1442 }
1443
1444 xx = x;
1445 yy = WRAP2(y-1,max_y);
1446 xxx = x;
1447 yyy = WRAP2(y+1,max_y);
1448
1449 if ((current>*(in+xxx+yyy*max_x)) &&
1450 (current>*(in+xx+yy*max_x))) {
1451 *(output+x+y*max_x) = foreground_value;
1452 continue;
1453 }
1454
1455 xx = WRAP2(x+1,max_x);
1456 yy = WRAP2(y-1,max_y);
1457 xxx = WRAP2(x-1,max_x);
1458 yyy = WRAP2(y+1,max_y);
1459
1460 if ((current>*(in+xxx+yyy*max_x)) &&
1461 (current>*(in+xx+yy*max_x))) {
1462 *(output+x+y*max_x) = foreground_value;
1463 continue;
1464 }
1465
1466 xx = WRAP2(x+1,max_x);
1467 yy = y;
1468 xxx = WRAP2(x-1,max_x);
1469 yyy = y;
1470
1471 if ((current>*(in+xxx+yyy*max_x)) &&
1472 (current>*(in+xx+yy*max_x))) {
1473 *(output+x+y*max_x) = foreground_value;
1474 continue;
1475 }
1476
1477 *(output+x+y*max_x) = background_value;
1478 }
1479 }
1480
spread_over_max(hf_type * input,hf_type * output,gint max_x,gint max_y,gboolean wrap,gint min_width,gint max_width,hf_type min_val,hf_type max_val,gint function)1481 void spread_over_max (hf_type *input, hf_type *output, gint max_x, gint max_y, gboolean wrap, gint min_width, gint max_width, hf_type min_val, hf_type max_val, gint function) {
1482 // Widen the areas in input which are not null
1483 // "Spread" them depending on the pixel value
1484 // Typically: used for thickening lines with a variable width
1485 // We scale the values from 1 to max_width, depending on min_val and max_val
1486 // If both min_val and max_val are 0, or if max_val < min_val,
1487 // we recalculate them from *input
1488
1489 gint i, x, y, index_x, index_y, index_i;
1490 gdouble ratio;
1491 gint value;
1492 hf_type *tmp;
1493
1494 if (((!max_val) && (!min_val)) || (max_val<min_val)) {
1495 // Note: min_val is the lowest value greater than 0
1496 min_val = MAX_HF_VALUE;
1497 max_val = 0;
1498 for (i=0; i<(max_y*max_x); i++) {
1499 if (!*(input+i))
1500 continue;
1501 if (*(input+i)<min_val)
1502 min_val = *(input+i);
1503 if (*(input+i)>max_val)
1504 max_val = *(input+i);
1505 }
1506 }
1507 if (max_val==min_val) {
1508 printf("Nothing to spread!\n");
1509 return;
1510 }
1511
1512 ratio = ((gdouble) max_width) / (gdouble) max_val ;
1513
1514 // printf("MIN_VAL: %d; MAX_VAL: %d; MIN_WIDTH: %d; MAX_WIDTH: %d; RATIO: %7.3f\n",min_val, max_val, min_width, max_width, ratio);
1515
1516 tmp = (hf_type *) x_calloc(max_x*max_y,sizeof(hf_type), "hf_type (tmp in spread_over_max)");
1517
1518 // write_png ("output0.png", 16, input, max_x, max_y);
1519
1520 // write_png ("input.png", 16, input, max_x, max_y);
1521 // This is a separable algorithm, with a half-kernel of 'radius' width
1522 // When wrap == FALSE, we simply skip the outbounded indices
1523 // 1st pass
1524 for (y=0; y<max_y; y++) {
1525 for (x=0; x<max_x; x++) {
1526 for (i=x; i<(x+max_width); i++) {
1527 if (!wrap) {
1528 if ((i<0) || (i>=max_x))
1529 continue;
1530 }
1531 index_x = VECTORIZE(x,y,max_x);
1532 index_i = VECTORIZE(WRAP2(i,max_x),y,max_x);
1533 // When ABS(i-x) <= normalize(*input+index_i), we output normalize(*input+index_i) in output+index_x, if it's >
1534 value = (hf_type) (0.5 + ratio * ((gdouble) *(input+index_i)));
1535 if (value)
1536 value = MAX(value,min_width);
1537 else
1538 continue;
1539 if (ABS(i-x)<=value) {
1540 value = value - ABS(i-x); // TENT_FUNC
1541 if (value>*(tmp+index_x))
1542 *(tmp+index_x) = value;
1543 }
1544 }
1545 }
1546 }
1547
1548 // write_png ("output1.png", 16, tmp, max_x, max_y);
1549
1550 for (x=0; x<max_x; x++) {
1551 for (y=0; y<max_y; y++) {
1552 for (i=y; i<(y+max_width); i++) {
1553 if (!wrap) {
1554 if ((i<0) || (i>=max_y))
1555 continue;
1556 }
1557 index_y = VECTORIZE(x,y,max_x);
1558 index_i = VECTORIZE(x,WRAP2(i,max_y),max_x);
1559 value = *(tmp+index_i);
1560 if (!value)
1561 continue;
1562 if (ABS(i-y)<=value) {
1563 value = value - ABS(i-y); // TENT_FUNC
1564 if (value>*(output+index_y))
1565 *(output+index_y) = value;
1566 }
1567 }
1568 }
1569 }
1570
1571 // write_png ("output2.png", 16, output, max_x, max_y);
1572 // memcpy(output,tmp,sizeof(hf_type)*max_x*max_y);
1573 x_free(tmp);
1574 }
1575
improve_edges(hf_type * hf_in,hf_type * hf_out,gint max_x,gint max_y,gboolean wrap,hf_type threshold)1576 void improve_edges (hf_type *hf_in, hf_type *hf_out, gint max_x, gint max_y, gboolean wrap, hf_type threshold) {
1577
1578 // Improve edges (or lines) by smoothing against a threshold value
1579 // Works only with dark edges
1580 //
1581 // Pixel numbering:
1582 // 4 | 3 | 2
1583 // 5 | 0 | 1
1584 // 6 | 7 | 8
1585 //
1586 // Edges type == {DARK_EDGES, BRIGHT_EDGES}
1587 //
1588 gint i,j, ii, startx, starty, endx, endy, index;
1589 hf_type v0, v1, v2, v3, v4, v5, v6, v7, v8;
1590
1591 if (wrap) {
1592 startx = starty = 0;
1593 endx = max_x;
1594 endy = max_y;
1595 }
1596 else {
1597 startx = starty = 1;
1598 endx = max_x - 1;
1599 endy = max_y - 1;
1600 }
1601
1602 // We test each axis (N-S, E-W, diagonal) to replace v0 so that it is coherent with its neighbours
1603
1604 // 2006-04-04 2-pass version - the first one works best when input==output
1605
1606 // This pass fills "holes", for instance (B = background, F = foreground):
1607 // B F B
1608 // F B B
1609 // becomes:
1610 // B F B
1611 // F F B
1612 for (i=0; i<max_y; i++) {
1613 for (j=0; j<max_x; j++) {
1614 ii = i*max_x;
1615 index = ii + j;
1616 v0 = *(hf_in + index);
1617 v1 = *(hf_in + ii + WRAP2(j+1,max_x));
1618 v5 = *(hf_in + ii + WRAP2(j-1,max_x));
1619 ii = WRAP2(i-1,max_y)*max_x;
1620 v3 = *(hf_in + ii + j);
1621 v2 = *(hf_in + ii + WRAP2(j+1,max_x));
1622 v4 = *(hf_in + ii + WRAP2(j-1,max_x));
1623 ii = WRAP2(i+1,max_y)*max_x;
1624 v7 = *(hf_in + ii + j);
1625 v8 = *(hf_in + ii + WRAP2(j+1,max_x));
1626 v6 = *(hf_in + ii + WRAP2(j-1,max_x));
1627 if (v0>threshold) {
1628 if (((v3<=threshold) && (v1<=threshold) && (v2>threshold)) ) {
1629 *(hf_in+index) = hf_type_avrg(v1,v3);
1630 continue;
1631 }
1632 if ((v1<=threshold) && (v7<=threshold) && (v8>threshold)){
1633 *(hf_in+index) = hf_type_avrg(v7,v1);
1634 continue;
1635 }
1636 if ((v3<=threshold) && (v5<=threshold) && (v4>threshold)) {
1637 *(hf_in+index) = hf_type_avrg(v3,v5);
1638 continue;
1639 }
1640 if ((v5<=threshold) && (v7<=threshold) && (v6>threshold)){
1641 *(hf_in+index) = hf_type_avrg(v5,v7);
1642 continue;
1643 }
1644
1645 }
1646 } // j (x) loop
1647 } // i (y) loop
1648
1649 // This pass tests that kind of patterns (h or v):
1650 // 2
1651 // 1 1 1
1652 // Here we remove "2" - being the background
1653
1654 for (i=0; i<max_y; i++) {
1655 for (j=0; j<max_x; j++) {
1656 ii = i*max_x;
1657 index = ii + j;
1658 v0 = *(hf_in + index);
1659 v1 = *(hf_in + ii + WRAP2(j+1,max_x));
1660 v5 = *(hf_in + ii + WRAP2(j-1,max_x));
1661 ii = WRAP2(i-1,max_y)*max_x;
1662 v3 = *(hf_in + ii + j);
1663 v2 = *(hf_in + ii + WRAP2(j+1,max_x));
1664 v4 = *(hf_in + ii + WRAP2(j-1,max_x));
1665 ii = WRAP2(i+1,max_y)*max_x;
1666 v7 = *(hf_in + ii + j);
1667 v8 = *(hf_in + ii + WRAP2(j+1,max_x));
1668 v6 = *(hf_in + ii + WRAP2(j-1,max_x));
1669 *(hf_out+index) = v0;
1670 if (v0<=threshold) {
1671 if ((v3>threshold) && (v7>threshold)) {
1672 if (((v2>threshold) && (v1>threshold) && (v8>threshold)) ||((v4>threshold) && (v5>threshold) && (v6>threshold)) ) {
1673 *(hf_out+index) = hf_type_avrg(v3,v7);
1674 continue;
1675 }
1676 }
1677 if ((v5>threshold) && (v1>threshold)) {
1678 if (((v2>threshold) && (v3>threshold) && (v4>threshold)) ||((v6>threshold) && (v7>threshold) && (v8>threshold)) ) {
1679 *(hf_out+index) = hf_type_avrg(v5,v1);
1680 }
1681 continue;
1682 }
1683 } // Threshold test
1684 } // j (x) loop
1685 } // i (y) loop
1686
1687 // memcpy(hf_out,hf_in,max_x*max_y*sizeof(hf_type));
1688 }
1689
hf_integrate(hf_struct_type * hf,gint angle)1690 void hf_integrate (hf_struct_type *hf, gint angle) {
1691 // Integrates an image from an arbitrary angle, in degrees
1692 // Used basically for "simple view modelling"
1693 // - Building a HF from a terrain picture with lateral lighting
1694 // 1. Recalculate min, max, average
1695 // 2. Oversize the HF (for rotating it without loosing the edges)
1696 // 3. Convert the HF to floating point
1697 // 4. Rotate the HF
1698 // 5. Integrates
1699 // 6. Rotate back the HF
1700 // 7. Crop back the HF
1701 // 8. Convert back the HF to hf_type
1702
1703 gdouble *buffer_in, *buffer_out, a, avrg, dmin, dmax, ratio;
1704 gint i, j, max, offset, index;
1705
1706 if (!hf->tmp_buf)
1707 hf_backup(hf);
1708 hf_min_max (hf);
1709 avrg = (gdouble) hf->avrg;
1710
1711 max = (((gint) (1.7*hf->max_x))>>1)<<1; // Must be even
1712 offset = max>>2;
1713 buffer_in = xcalloc_double_hf(max, "gdouble (buffer_in in hf_integrate)");
1714 buffer_out = xcalloc_double_hf(max, "gdouble (buffer_out in hf_integrate)");
1715
1716 for (i=0; i<hf->max_x; i++)
1717 for (j=0; j<hf->max_y; j++)
1718 // The content is centered
1719 // The average is adjusted to 0
1720 *(buffer_in+VECTORIZE(i+offset,j+offset,max)) = (gdouble) *(hf->tmp_buf+VECTORIZE(i,j,hf->max_x));
1721
1722 a = PI * ((gdouble) angle) / 180.0;
1723 rotate (cos(a), sin(a), (gpointer) buffer_in, (gpointer) buffer_out, max, max, GDOUBLE_ID, OVERFLOW_ZERO);
1724
1725 /*****************************************/
1726 // Integration
1727 // We adjust the resulting values to 0-MAX_HF_VALUE
1728
1729 dmin = 0.0;
1730 dmax = 0.0;
1731
1732 for (i=0; i<max; i++) {
1733 if (*(buffer_out+i)<dmin)
1734 dmin = *(buffer_out+i);
1735 if (*(buffer_out+i)>dmax)
1736 dmax = *(buffer_out+i);
1737 for (j=1; j<max; j++) {
1738 index = VECTORIZE(i,j,max);
1739 *(buffer_out + index) += *(buffer_out+VECTORIZE(i,j-1,max)) / 1.05;
1740
1741 if (*(buffer_out+index)<dmin)
1742 dmin = *(buffer_out+index);
1743 if (*(buffer_out+index)>dmax)
1744 dmax = *(buffer_out+index);
1745 }
1746 }
1747
1748
1749 /*****************************************/
1750
1751 a = -PI * ((gdouble) angle) / 180.0;
1752 rotate (cos(a), sin(a), (gpointer) buffer_out, (gpointer) buffer_in, max, max, GDOUBLE_ID, OVERFLOW_ZERO);
1753
1754 ratio = ((gdouble) MAX_HF_VALUE ) / (dmax - dmin);
1755
1756 for (i=0; i<hf->max_x; i++)
1757 for (j=0; j<hf->max_y; j++)
1758 // The content is translated back to (0,0) (cropped)
1759 // and scaled back to 0-MAX_HF_VALUE
1760 *(hf->hf_buf+VECTORIZE(i,j,hf->max_x)) = (hf_type) (ratio * ( *(buffer_in+VECTORIZE(i+offset,j+offset,max)) - dmin ));
1761
1762 x_free(buffer_in);
1763 x_free(buffer_out);
1764 }
1765
hf_cut_graph(hf_struct_type * hf,gint index,gint axis)1766 void hf_cut_graph (hf_struct_type *hf, gint index, gint axis) {
1767
1768 // Draw a graph of line / column "index" on axis "axis"
1769
1770 guint i,j, shift;
1771 hf_type value;
1772
1773 if (!hf->tmp_buf)
1774 hf_backup(hf);
1775
1776 index = MAX(0,MIN(hf->max_x,index));
1777
1778 shift = 1 + log2i(MAX_HF_VALUE) - log2i(hf->max_x);
1779
1780 for (i=0; i<hf->max_x; i++) {
1781 if (axis==V_AXIS) // index is an X value
1782 value = *(hf->tmp_buf+VECTORIZE(index,i,hf->max_x));
1783 else // index is a Y value
1784 value = *(hf->tmp_buf+VECTORIZE(i,index,hf->max_x));
1785 if (shift)
1786 value = value>>shift;
1787 for (j=0; j<hf->max_x; j++) // Y output
1788 if ((hf->max_y-value)>j)
1789 *(hf->hf_buf+VECTORIZE(i,j,hf->max_x)) = 0;
1790 else
1791 *(hf->hf_buf+VECTORIZE(i,j,hf->max_x)) = MAX_HF_VALUE>>1;
1792 }
1793
1794 }
1795
hf_histogram(hf_struct_type * hf)1796 void hf_histogram (hf_struct_type *hf) {
1797
1798 // Draw a values histogram in hf->hf_buf
1799
1800 guint i,j, shift;
1801
1802 hf_type value;
1803 gulong *vector, max;
1804
1805 vector = (gulong *) x_calloc(sizeof(gulong),hf->max_x, "gulong (vector in hf_histogram)");
1806
1807 if (!hf->tmp_buf)
1808 hf_backup(hf);
1809
1810 shift = 1 + log2i(MAX_HF_VALUE) - log2i(hf->max_x);
1811
1812 for (i=0; i<(hf->max_x*hf->max_y); i++)
1813 *(vector+((*(hf->tmp_buf+i))>>shift)) += 1;
1814
1815 max = 0;
1816
1817 for (i=0; i<hf->max_x; i++)
1818 if (*(vector+i)>max)
1819 max = *(vector+i);
1820
1821 if (max>hf->max_y)
1822 shift = 1 + log2i(max) - log2i(hf->max_x);
1823
1824 // printf("MAX: %d; SHIFT: %d\n",max, shift);
1825
1826 for (i=0; i<hf->max_x; i++) {
1827 value = *(vector+i);
1828 if (shift)
1829 value = value>>shift;
1830 for (j=0; j<hf->max_y; j++) // Y output
1831 if ((hf->max_y-value)>j)
1832 *(hf->hf_buf+VECTORIZE(i,j,hf->max_x)) = 0;
1833 else
1834 *(hf->hf_buf+VECTORIZE(i,j,hf->max_x)) = MAX_HF_VALUE>>1;
1835 }
1836 x_free(vector);
1837 }
1838
histogram(hf_struct_type * hf,guint width,guint height,guchar * output)1839 gboolean histogram (hf_struct_type *hf, guint width, guint height, guchar *output) {
1840
1841 // Builds a histogram in output
1842
1843 gint i,j, k, value;
1844 gdouble ratio;
1845 gulong *vector, max;
1846
1847 if (!output)
1848 return FALSE;
1849
1850 vector = (gulong *) x_calloc(sizeof(gulong),width, "gulong (vector in histogram)");
1851
1852 ratio = ((gdouble) MAX_HF_VALUE) / (gdouble) width;
1853
1854 for (i=0; i<(hf->max_x*hf->max_y); i++) {
1855 k = (gint) (((gdouble)*(hf->hf_buf+i))/ratio);
1856 if (k>=width) {
1857 printf("WARNING! Vector index in histogram out of boundary: %d >= %d\n",k, width);
1858 continue;
1859 }
1860 *(vector+k) += 1;
1861 }
1862
1863 // Smooth the vector and calculate the max number of pixels
1864 /*
1865 last = *vector;
1866 for (i=0; i<width; i++) {
1867 llast = *(vector+i);
1868 *(vector+i) += last;
1869 last = llast;
1870 }
1871 */
1872 max = 0;
1873
1874 for (i=0; i<width; i++)
1875 if (*(vector+i)>max)
1876 max = *(vector+i);
1877
1878 ratio = ((gdouble) max) / (gdouble) height;
1879
1880 for (i=0; i<width; i++) {
1881 value = (guint) (((gdouble) *(vector+i)) / ratio);
1882 value = height - value;
1883 for (j=(height-1); j>value; j--) { // Y output
1884 k = VECTORIZE(i,j,width);
1885 if ((k<0) || (k>=(width*height))) {
1886 printf("WARNING! Output index in histogram out of boundary: %i >= %d\n",k, width*height);
1887 continue;
1888 }
1889 *(output+k) = 127;
1890 }
1891 }
1892
1893 x_free(vector);
1894 return TRUE;
1895 }
1896