1 /* hf_calc.c - hf transformations utilities
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 "hf.h"
23 #include "hf_calc.h"
24 #include "hf_filters.h"
25 #include "fill.h"
26
pow2_scale_grid(gpointer in,gpointer out,gint max_x,gint max_y,gint log_ratio,gint data_type_in,gint data_type_out)27 void pow2_scale_grid (gpointer in, gpointer out, gint max_x, gint max_y, gint log_ratio, gint data_type_in, gint data_type_out) {
28
29 // Scale input_grid in output_grid, given log_ratio
30 // log_ratio = -1 : reduce by 50%
31 // log_ratio = -2 : reduce by 75%
32 // log_ratio = 1 : increase size by 200%
33 // log_ratio = 2 : increase size by 400%...
34 // output_grid should be allocated
35 // max_x, max_y: size of the input grid
36 // m_x, m_y: deduced size of the output grid
37 gint i, j, indx, m_x, m_y, m_i;
38 m_x = LEFT_SHIFT(max_x,log_ratio);
39 m_y = LEFT_SHIFT(max_y,log_ratio);
40 // printf("MX: %d; MY: %d\n",m_x,m_y);
41 for (i = 0; i < m_x; i++)
42 for (j = 0; j < m_y; j++) {
43
44 indx = RIGHT_SHIFT(i,log_ratio) +
45 RIGHT_SHIFT(j,log_ratio)*max_x;
46 m_i = j*m_y + i;
47
48 switch (data_type_in) {
49 case GDOUBLE_ID:
50
51 switch (data_type_out) {
52 case GDOUBLE_ID:
53 *(((gdouble *) out) +m_i) = *(((gdouble *) in)+indx);
54 break;
55 case HF_TYPE_ID:
56 *(((hf_type *) out) +m_i) = (hf_type) *(((gdouble *) in)+indx);
57 break;
58 case GINT_ID:
59 *(((gint *) out) +m_i) = (gint) *(((gdouble *) in)+indx);
60 break;
61 default: // UNSIGNED_CHAR_ID:
62 *(((unsigned char *) out) +m_i) = (unsigned char) *(((gdouble *) in)+indx);
63 } // End switch (data_type_out)
64 break;
65
66 case HF_TYPE_ID:
67
68 switch (data_type_out) {
69 case GDOUBLE_ID:
70 *(((gdouble *) out) +m_i) = (gdouble) *(((hf_type *) in)+indx);
71 break;
72 case HF_TYPE_ID:
73 *(((hf_type *) out) +m_i) = *(((hf_type *) in)+indx);
74 break;
75 case GINT_ID:
76 *(((gint *) out) +m_i) = (gint) *(((hf_type *) in)+indx);
77 break;
78 default: // UNSIGNED_CHAR_ID:
79 *(((unsigned char *) out) +m_i) = *(((hf_type *) in)+indx) >> 8;
80
81 // if ((i==127) && !(j%64))
82 // printf("Value in: %d; Value out: %d\n",*(((hf_type *) in)+indx),*(((unsigned char *) out) +m_i));
83
84 } // End switch (data_type_out)
85 break;
86
87 case GINT_ID:
88
89 switch (data_type_out) {
90 case GDOUBLE_ID:
91 *(((gdouble *) out) +m_i) = (gdouble) *(((gint *) in)+indx);
92 break;
93 case HF_TYPE_ID:
94 *(((hf_type *) out) +m_i) = (hf_type) *(((gint *) in)+indx);
95 break;
96 case GINT_ID:
97 *(((gint *) out) +m_i) = *(((gint *) in)+indx);
98 break;
99 default: // UNSIGNED_CHAR_ID:
100 *(((unsigned char *) out) +m_i) = (unsigned char) *(((gint *) in)+indx) >> 8;
101 } // End switch (data_type_out)
102
103 break;
104 default: // UNSIGNED_CHAR_ID:
105
106 switch (data_type_out) {
107 case GDOUBLE_ID:
108 *(((gdouble *) out) +m_i) = (gdouble) *(((unsigned char *) in)+indx);
109 break;
110 case HF_TYPE_ID:
111 *(((hf_type *) out) +m_i) = (hf_type) *(((unsigned char *) in)+indx);
112 break;
113 case GINT_ID:
114 *(((gint *) out) +m_i) = (gint) *(((unsigned char *) in)+indx);
115 break;
116 default: // UNSIGNED_CHAR_ID:
117 *(((unsigned char *) out) +m_i) = *(((unsigned char *) in)+indx);
118 } // End switch (data_type_out)
119
120 } // End switch (data_type_in)
121 }
122
123 }
124
pow2_scale_uchar_grid_with_mask(gpointer in,unsigned char * out,gint max_x,gint max_y,gint log_ratio,gint data_type_in,gdouble * mask,gpointer mask_value_ptr)125 void pow2_scale_uchar_grid_with_mask (gpointer in, unsigned char *out, gint max_x, gint max_y, gint log_ratio, gint data_type_in, gdouble *mask, gpointer mask_value_ptr) {
126
127 // Scale input_grid in output_grid, given log_ratio
128 // Use a double mask, translated as a a RGB value, given mask_value_ptr
129 // Range of the mask: cut to (0.0 - 1.0)
130 // Input: any
131 // Output: (pix8_rgb *) == 3 * (unsigned char *)
132 //
133 // log_ratio = -1 : reduce by 50%
134 // log_ratio = -2 : reduce by 75%
135 // log_ratio = 1 : increase size by 200%
136 // log_ratio = 2 : increase size by 400%...
137 // output_grid should be allocated
138 // max_x, max_y: size of the input grid
139 // m_x, m_y: deduced size of the output grid
140 gint i, j, indx, m_x, m_y, m_i;
141 gdouble mask_level;
142 pix8_rgb *mask_value;
143 unsigned char value;
144 mask_value = (pix8_rgb *) mask_value_ptr;
145 m_x = LEFT_SHIFT(max_x,log_ratio);
146 m_y = LEFT_SHIFT(max_y,log_ratio);
147 // printf("MX: %d; MY: %d\n",m_x,m_y);
148 for (i = 0; i < m_x; i++)
149 for (j = 0; j < m_y; j++) {
150
151 indx = RIGHT_SHIFT(i,log_ratio) +
152 RIGHT_SHIFT(j,log_ratio)*max_x;
153 m_i = j*m_y + i;
154
155 switch (data_type_in) {
156 case GDOUBLE_ID:
157 value = (unsigned char) *(((gdouble *) in)+indx);
158 break;
159 case HF_TYPE_ID:
160 value = *(((hf_type *) in)+indx) >> 8;
161 break;
162 case GINT_ID:
163 value = *(((gint *) in)+indx);
164 break;
165 default: // UNSIGNED_CHAR_ID:
166 value = *(((unsigned char *) in)+indx);
167 }
168 if (mask) {
169 // When mask > 1.0, we revert the colour
170 if (1.0 < *(mask+indx)) {
171 mask_level = 1.0 - MIN(MAX(*(mask+indx)-1.0, 0.0),1.0);
172 // Red
173 *(out + m_i * 3) = (unsigned char) ((mask_level * (gdouble) value) + (1.0-mask_level) * (gdouble) ((~mask_value->r) & value));
174 // Green
175 *(out + 1 + m_i * 3) = (unsigned char) ((mask_level * (gdouble) value) + (1.0-mask_level) * (gdouble) ((~mask_value->g) & value));
176 // Blue
177 *(out + 2 + m_i * 3) = (unsigned char) ((mask_level * (gdouble) value) + (1.0-mask_level) * (gdouble) ((~mask_value->b) & value));
178 }
179 else {
180 mask_level = MIN(MAX(*(mask+indx),0.0),1.0);
181 // Red
182 *(out + m_i * 3) = (unsigned char) ((mask_level * (gdouble) value) + (1.0-mask_level) * (gdouble) (mask_value->r & value));
183 // Green
184 *(out + 1 + m_i * 3) = (unsigned char) ((mask_level * (gdouble) value) + (1.0-mask_level) * (gdouble) (mask_value->g & value));
185 // Blue
186 *(out + 2 + m_i * 3) = (unsigned char) ((mask_level * (gdouble) value) + (1.0-mask_level) * (gdouble) (mask_value->b & value));
187 }
188 }
189 else {
190 *(out + m_i * 3) = mask_value->r & value;
191 *(out + 1 + m_i * 3) = mask_value->g & value;
192 *(out + 2 + m_i * 3) = mask_value->b & value;
193 }
194 }
195 }
196
hf_scale(hf_struct_type * hf,gint log_ratio)197 hf_struct_type * hf_scale (hf_struct_type *hf, gint log_ratio) {
198 // Scale a HF...
199 // log_ratio = -1 : reduce by 50%
200 // log_ratio = -2 : reduce by 75%
201 // log_ratio = 1 : increase size by 200%
202 // log_ratio = 2 : increase size by 400%...
203 hf_struct_type *out;
204 gint i, j, max_x, max_y;
205 max_x = LEFT_SHIFT(hf->max_x,log_ratio);
206 max_y = LEFT_SHIFT(hf->max_y,log_ratio);
207 out = hf_new(max_x);
208 for (i = 0; i < max_x; i++)
209 for (j = 0; j < max_y; j++) {
210 *(out->hf_buf + j*max_y + i) =
211 *(hf->hf_buf + (RIGHT_SHIFT(i,log_ratio)) +
212 (RIGHT_SHIFT(j,log_ratio))*hf->max_x) ;
213 }
214 return out;
215 }
216
hf_merge_buf(hf_type * hf_pencil,gint pen_max_x,gint pen_max_y,hf_type * hf,gint max_x,gint max_y,int x,int y,gint pen_merge,gboolean pen_tiles,gfloat h_displacement,gboolean normalize)217 void hf_merge_buf (hf_type *hf_pencil, gint pen_max_x, gint pen_max_y,
218 hf_type *hf, gint max_x, gint max_y,
219 int x, int y,
220 gint pen_merge,
221 gboolean pen_tiles,
222 gfloat h_displacement,
223 gboolean normalize) {
224
225 // Writes hf_pencil (or any other HF) in hf centered at (x,y)
226 // For instance, when drawing, (x,y) is the mouse position
227
228 // h_displacement is the relative level of the added HF, from 0.0 to 1.0
229 // normalize = TRUE applies to pen_merge = ADD
230 // ... reduce the added HFs by half, to avoid "burning"
231
232 gint i,j, k, l, displ; // displ = displacement (we center the HF on the mouse position)
233 glong val;
234 hf_type max_value = 1; // Should be > 0, even if HF is black
235 if (pen_merge == MULTIPLY) {
236 for (i=0; i<(max_x*max_y); i++)
237 if (max_value<*(hf+i))
238 max_value = *(hf+i);
239 }
240 displ = pen_max_x >> 1; // Pencil is square!
241 x = x - displ;
242 y = y - displ;
243 for (i = 0; i < pen_max_x ; i++) {
244 // Calculate X index in hf world
245 k = x + i;
246 if (!pen_tiles) {
247 if (k>=max_x)
248 break;
249 else
250 if (k<0)
251 continue;
252 }
253 if (k < 0)
254 k = max_x + k;
255 else
256 k = k%max_x;
257
258 for (j = 0; j < pen_max_y; j++) {
259 // Calculate Y index in hf world
260 l = y + j;
261 if (!pen_tiles) {
262 if (l>=max_y)
263 break;
264 else
265 if (l<0)
266 continue;
267 }
268 if (l < 0)
269 l = max_y + l;
270 else
271 l = l%max_y;
272 if (pen_merge == SUBTRACT)
273 val = (*(hf + k + l*max_x)) - h_displacement *
274 (*(hf_pencil + i + j * pen_max_x));
275 else if (pen_merge == MULTIPLY) {
276 val = (glong) (((gfloat)*(hf + k + l*max_x)) *
277 ((gfloat) *(hf_pencil + i + j * pen_max_x)) / max_value );
278 val = h_displacement * val + (1.0-h_displacement) * *(hf + k + l*max_x);
279 }
280 else // pen_merge == ADD
281 if (normalize)
282 val = ((1.0-h_displacement)* *(hf + k + l*max_x)) + h_displacement *
283 (*(hf_pencil + i + j * pen_max_x));
284 else
285 val = (*(hf + k + l*max_x)) + h_displacement *
286 (*(hf_pencil + i + j * pen_max_x));
287 *(hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),0xFFFF);
288 }
289 }
290 }
291
292
hf_merge(hf_struct_type * hf_pencil,hf_struct_type * hf,int x,int y,gint pen_merge,gboolean pen_tiles,gfloat h_displacement,gboolean normalize)293 void hf_merge (hf_struct_type *hf_pencil, hf_struct_type *hf,
294 int x, int y,
295 gint pen_merge,
296 gboolean pen_tiles,
297 gfloat h_displacement,
298 gboolean normalize) {
299
300 // See hf_merge_buf
301
302 hf_merge_buf (hf_pencil->hf_buf, hf_pencil->max_x, hf_pencil->max_y,
303 hf->hf_buf, hf->max_x, hf->max_y, x, y,
304 pen_merge, pen_tiles, h_displacement, normalize);
305
306 }
307
generalized_merge(gpointer map,gint map_data_type,gint size_x,gint size_y,gpointer hf,gint hf_data_type,gint max_x,gint max_y,gint x,gint y,gint merge_type,gboolean wrap,gboolean square_symmetry)308 void generalized_merge (gpointer map, gint map_data_type, gint size_x, gint size_y,
309 gpointer hf, gint hf_data_type, gint max_x, gint max_y, gint x, gint y,
310 gint merge_type, gboolean wrap, gboolean square_symmetry) {
311
312 // 2004-01 Merge a map of size_x * size_y at center position (x,y) in hf
313 // Used for drawing a continuous line and craters
314 // [Eventually generalize to encompass hf_merge]
315 // data_type: type of map, can be any of
316 // GINT_ID, HF_TYPE_ID, GDOUBLE_ID, UNSIGNED_CHAR_ID
317 // merge_type can be any of : ADD, SUBTRACT, MULTIPLY, MULTIPLY2
318 // square_symmetry is true only when the map has a square symmetry,
319 // so that one quadrant is sufficient to deduce it
320 // Implies size_x, size_y are odd, and generally that size_x == size_y
321 // The size of a square map of radius R is (2*R+1) * (2*R+1)
322
323 gint i,j,ix,iy, k, l, radius_x, radius_y,lx;
324 glong val;
325 gdouble ratio, dval=0.0;
326 gboolean if_val = TRUE;
327 static gdouble dmax = (gdouble) MAX_HF_VALUE;
328
329 radius_x = size_x >>1;
330 radius_y = size_y >>1;
331
332 x = x - radius_x;
333 y = y - radius_y;
334
335 if (square_symmetry)
336 lx = radius_x + 1;
337 else
338 lx = size_x;
339
340 for (i = 0; i < size_x ; i++) {
341 // Calculate X index in hf world
342 k = x + i;
343 if (!wrap) {
344 if (k>=max_x)
345 break;
346 else
347 if (k<0)
348 continue;
349 }
350 if (k < 0)
351 k = max_x + k;
352 else
353 k = k%max_x;
354
355 for (j = 0; j < size_y; j++) {
356 // Calculate Y index in hf world
357 l = y + j;
358 if (!wrap) {
359 if (l>=max_y)
360 break;
361 else
362 if (l<0)
363 continue;
364 }
365 if (l < 0)
366 l = max_y + l;
367 else
368 l = l%max_y;
369
370 ix = i;
371 iy = j;
372 if (square_symmetry) {
373 if (i>radius_x)
374 ix = 2*radius_x - i;
375 if (j>radius_y)
376 iy = 2*radius_y - j;
377 }
378 // printf("(radius_y,radius_x): (%d,%d); (i,j): (%d,%d); (ix,iy): (%d,%d); lx: %d;\n",radius_x, radius_y,i,j, ix,iy, lx);
379 switch (map_data_type) {
380 case GINT_ID:
381 val = (glong) *(((gint *) map) + ix + iy * lx);
382 break;
383 case HF_TYPE_ID:
384 val = (glong) *(((hf_type *) map) + ix + iy * lx);
385 break;
386 case GDOUBLE_ID:
387 dval = *(((gdouble *) map) + ix + iy * lx);
388 if_val = FALSE;
389 break;
390 case UNSIGNED_CHAR_ID:
391 val = (glong) *(((unsigned char *) map) + ix + iy * lx);
392 break;
393 default:
394 printf(_("Unexpected option in %s; contact the programmer!"),"generalized_merge");
395 printf("\n");
396 return;
397 }
398
399 if (if_val && (hf_data_type==GDOUBLE_ID))
400 dval = (gdouble) val;
401
402 switch (merge_type) {
403 case ADD:
404
405 switch (hf_data_type) {
406 case GINT_ID:
407 val = *((gint *) hf + k + l*max_x) + val;
408 *((gint *)hf + k + l*max_x) = (gint) MIN(MAX(val,0),MAX_HF_VALUE);
409 break;
410 case HF_TYPE_ID:
411 val = *((hf_type *) hf + k + l*max_x) + val;
412 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
413 break;
414 case GDOUBLE_ID:
415 dval = *((gdouble *) hf + k + l*max_x) + dval;
416 *((gdouble *)hf + k + l*max_x) = (gdouble) MIN(MAX(dval,0.0),dmax);
417 break;
418 case UNSIGNED_CHAR_ID:
419 val = *((unsigned char *) hf + k + l*max_x) + val;
420 *((unsigned char *)hf + k + l*max_x) = (unsigned char) MIN(MAX(val,0),MAX_HF_VALUE);
421 break;
422 default:
423 val = *((hf_type *) hf + k + l*max_x) + val;
424 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
425 }
426 break;
427 case SUBTRACT:
428 switch (hf_data_type) {
429 case GINT_ID:
430 val = *((gint *) hf + k + l*max_x) - val;
431 *((gint *)hf + k + l*max_x) = (gint) MIN(MAX(val,0),MAX_HF_VALUE);
432 break;
433 case HF_TYPE_ID:
434 val = *((hf_type *) hf + k + l*max_x) - val;
435 // printf("hf: %d\n", *((hf_type *) hf + k + l*max_x));
436 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
437 break;
438 case GDOUBLE_ID:
439 dval = *((gdouble *) hf + k + l*max_x) - dval;
440 *((gdouble *)hf + k + l*max_x) = (gdouble) MIN(MAX(dval,0.0),dmax);
441 break;
442 case UNSIGNED_CHAR_ID:
443 val = *((unsigned char *) hf + k + l*max_x) - val;
444 *((unsigned char *)hf + k + l*max_x) = (unsigned char) MIN(MAX(val,0),MAX_HF_VALUE);
445 break;
446 default:
447 val = *((hf_type *) hf + k + l*max_x) - val;
448 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
449 }
450 break;
451 case MULTIPLY:
452 case MULTIPLY2:
453 if (merge_type == MULTIPLY)
454 ratio = ((gdouble) val) / dmax;
455 else
456 ratio = 1.0 + 0.5 * ((gdouble) val) / dmax;
457 switch (hf_data_type) {
458 case GINT_ID:
459 val = (glong) (ratio * (gdouble) *((gint *) hf + k + l*max_x));
460 *((gint *)hf + k + l*max_x) = (gint) MIN(MAX(val,0),MAX_HF_VALUE);
461 break;
462 case HF_TYPE_ID:
463 val =(glong) (ratio * (gdouble) *((hf_type *) hf + k + l*max_x));
464 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
465 break;
466 case GDOUBLE_ID:
467 if (merge_type == MULTIPLY)
468 ratio = dval / dmax;
469 else
470 ratio = 1.0 + 0.5 * dval / dmax;
471 dval = (ratio * *((gdouble *) hf + k + l*max_x));
472 *((gdouble *)hf + k + l*max_x) = (gdouble) MIN(MAX(dval,0.0),dmax);
473 break;
474 case UNSIGNED_CHAR_ID:
475 val = (glong) (ratio * (gdouble) *((unsigned char *) hf + k + l*max_x));
476 *((unsigned char *)hf + k + l*max_x) = (unsigned char) MIN(MAX(val,0),MAX_HF_VALUE);
477 break;
478 default:
479 val =(glong) (ratio * (gdouble) *((hf_type *) hf + k + l*max_x));
480 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
481 }
482 break;
483 case MAX_MERGE:
484 switch (hf_data_type) {
485 case GINT_ID:
486 val = MAX(*((gint *) hf + k + l*max_x),val);
487 *((gint *)hf + k + l*max_x) = (gint) MIN(MAX(val,0),MAX_HF_VALUE);
488 break;
489 case HF_TYPE_ID:
490 val = MAX(*((hf_type *) hf + k + l*max_x),val);
491 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
492 break;
493 case GDOUBLE_ID:
494 dval = MAX(*((gdouble *) hf + k + l*max_x),dval);
495 *((gdouble *)hf + k + l*max_x) = (gdouble) MIN(MAX(dval,0.0),dmax);
496 break;
497 case UNSIGNED_CHAR_ID:
498 val = MAX(*((unsigned char *) hf + k + l*max_x),val);
499 *((unsigned char *)hf + k + l*max_x) = (unsigned char) MIN(MAX(val,0),MAX_HF_VALUE);
500 break;
501 default:
502 val = MAX(*((hf_type *) hf + k + l*max_x),val);
503 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
504 }
505 break;
506 case MIN_MERGE:
507 switch (hf_data_type) {
508 case GINT_ID:
509 val = MIN(*((gint *) hf + k + l*max_x),val);
510 *((gint *)hf + k + l*max_x) = (gint) MIN(MAX(val,0),MAX_HF_VALUE);
511 break;
512 case HF_TYPE_ID:
513 val = MIN(*((hf_type *) hf + k + l*max_x),val);
514 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
515 break;
516 case GDOUBLE_ID:
517 dval = MIN(*((gdouble *) hf + k + l*max_x),dval);
518 *((gdouble *)hf + k + l*max_x) = (gdouble) MIN(MAX(dval,0.0),dmax);
519 break;
520 case UNSIGNED_CHAR_ID:
521 val = MIN(*((unsigned char *) hf + k + l*max_x),val);
522 *((unsigned char *)hf + k + l*max_x) = (unsigned char) MIN(MAX(val,0),MAX_HF_VALUE);
523 break;
524 default:
525 val = MIN(*((hf_type *) hf + k + l*max_x),val);
526 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
527 }
528 break;
529 default:
530 printf(_("Unexpected option in %s; contact the programmer!"),"generalized_merge");
531 printf("\n");
532 return;
533 }
534 }
535 }
536 }
537
interpolated_merge(gpointer map,gint map_data_type,gint size_x,gint size_y,gpointer hf,gint hf_data_type,gint max_x,gint max_y,gdouble x,gdouble y,gint merge_type,gboolean wrap,gboolean square_symmetry)538 void interpolated_merge (gpointer map,
539 gint map_data_type, gint size_x, gint size_y,
540 gpointer hf, gint hf_data_type, gint max_x, gint max_y,
541 gdouble x, gdouble y,
542 gint merge_type, gboolean wrap, gboolean square_symmetry) {
543
544 // 2004-01 Merge similar to generalized merge
545 // (x,y) are real coordinates in HF
546 // We interpolate the pixels in the map to merge
547
548 gint i,j, k, l, radius_x, radius_y,lx, return_gint = 0;
549 unsigned char return_uns_char = 0;
550 glong val;
551 gboolean if_val=TRUE;
552 hf_type return_hf_type = 0;
553 gdouble ratio, dx, dy, ix, iy, return_gdouble = 0.0, dval = 0.0;
554 static gdouble dmax = (gdouble) MAX_HF_VALUE;
555
556 radius_x = size_x >>1;
557 radius_y = size_y >>1;
558
559 dx = x - floor(x);
560 dy = y - floor(y);
561 // printf("INTERPOLATED MERGE AT: (%5.2f, %5.2f); (DX,DY): (%5.2f,%5.2f)\n",x,y,dx,dy);
562 x = floor(x) - (gdouble) radius_x ;
563 y = floor(y) - (gdouble) radius_y;
564
565 if (square_symmetry)
566 lx = radius_x + 1;
567 else
568 lx = size_x;
569
570 for (i = 0; i < (size_x+1) ; i++) {
571 // Calculate X index in hf world
572 k = i + (gint) x;
573 if (!wrap) {
574 if (k>=max_x)
575 break;
576 else
577 if (k<0)
578 continue;
579 }
580 if (k < 0)
581 k = max_x + k;
582 else
583 k = k%max_x;
584
585 for (j = 0; j < (size_y+1); j++) {
586 // Calculate Y index in hf world
587 l = j + (gint) y;
588 if (!wrap) {
589 if (l>=max_y)
590 break;
591 else
592 if (l<0)
593 continue;
594 }
595 if (l < 0)
596 l = max_y + l;
597 else
598 l = l%max_y;
599
600 ix = -dx + (gdouble) i;
601 iy = -dy + (gdouble) j;
602 if (square_symmetry) {
603 if (ix>(gdouble) radius_x)
604 ix = 2.0*(gdouble) radius_x - ix;
605 if (iy>(gdouble) radius_y)
606 iy = 2.0*(gdouble) radius_y - iy;
607 }
608 // printf("(radius_y,radius_x): (%d,%d); (i,j): (%d,%d); (ix,iy): (%5.2f,%5.2f); lx: %d;\n",radius_x, radius_y,i,j, ix,iy, lx);
609 switch (map_data_type) {
610 case GINT_ID:
611 interpolate (ix, iy, map, lx, lx, &return_gint,
612 GINT_ID, OVERFLOW_ZERO);
613 val = (glong) return_gint;
614 break;
615 case HF_TYPE_ID:
616 interpolate (ix, iy, map, lx, lx, &return_hf_type,
617 HF_TYPE_ID, OVERFLOW_ZERO);
618 val = (glong) return_hf_type;
619 break;
620 case GDOUBLE_ID:
621 interpolate (ix, iy, map, lx, lx, &return_gdouble,
622 GDOUBLE_ID, OVERFLOW_ZERO);
623 // val = (glong) return_gdouble;
624 dval = return_gdouble;
625 if_val = FALSE;
626 break;
627 case UNSIGNED_CHAR_ID:
628 interpolate (ix, iy, map, lx, lx, &return_uns_char,
629 UNSIGNED_CHAR_ID, OVERFLOW_ZERO);
630 val = (glong) return_uns_char;
631 break;
632 default:
633 printf(_("Unexpected option in %s; contact the programmer!"),"interpolated_merge");
634 printf("\n");
635 return;
636 }
637
638 if (if_val && (hf_data_type==GDOUBLE_ID))
639 dval = (gdouble) val;
640
641 switch (merge_type) {
642 case ADD:
643 switch (hf_data_type) {
644 case GINT_ID:
645 val = *((gint *) hf + k + l*max_x) + val;
646 *((gint *)hf + k + l*max_x) = (gint) MIN(MAX(val,0),MAX_HF_VALUE);
647 break;
648 case HF_TYPE_ID:
649 val = *((hf_type *) hf + k + l*max_x) + val;
650 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
651 break;
652 case GDOUBLE_ID:
653 dval = *((gdouble *) hf + k + l*max_x) + dval;
654 *((gdouble *)hf + k + l*max_x) = (gdouble) MIN(MAX(dval,0.0),dmax);
655 break;
656 case UNSIGNED_CHAR_ID:
657 val = *((unsigned char *) hf + k + l*max_x) + val;
658 *((unsigned char *)hf + k + l*max_x) = (unsigned char) MIN(MAX(val,0),MAX_HF_VALUE);
659 break;
660 default:
661 val = *((hf_type *) hf + k + l*max_x) + val;
662 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
663 }
664 break;
665 case SUBTRACT:
666 switch (hf_data_type) {
667 case GINT_ID:
668 val = *((gint *) hf + k + l*max_x) - val;
669 *((gint *)hf + k + l*max_x) = (gint) MIN(MAX(val,0),MAX_HF_VALUE);
670 break;
671 case HF_TYPE_ID:
672 val = *((hf_type *) hf + k + l*max_x) - val;
673 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
674 break;
675 case GDOUBLE_ID:
676 dval = *((gdouble *) hf + k + l*max_x) - dval;
677 *((gdouble *)hf + k + l*max_x) = (gdouble) MIN(MAX(dval,0.0),dmax);
678 break;
679 case UNSIGNED_CHAR_ID:
680 val = *((unsigned char *) hf + k + l*max_x) - val;
681 *((unsigned char *)hf + k + l*max_x) = (unsigned char) MIN(MAX(val,0),MAX_HF_VALUE);
682 break;
683 default:
684 val = *((hf_type *) hf + k + l*max_x) - val;
685 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
686 }
687 break;
688 case MULTIPLY:
689 case MULTIPLY2:
690 if (merge_type == MULTIPLY)
691 ratio = ((gdouble) val) / dmax;
692 else
693 ratio = 1.0 + 0.5 * ((gdouble) val) / dmax;
694 switch (hf_data_type) {
695 case GINT_ID:
696 val = (glong) (ratio * (gdouble) *((gint *) hf + k + l*max_x));
697 *((gint *)hf + k + l*max_x) = (gint) MIN(MAX(val,0),MAX_HF_VALUE);
698 break;
699 case HF_TYPE_ID:
700 val =(glong) (ratio * (gdouble) *((hf_type *) hf + k + l*max_x));
701 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
702 break;
703 case GDOUBLE_ID:
704 if (merge_type == MULTIPLY)
705 ratio = dval / dmax;
706 else
707 ratio = 1.0 + 0.5 * dval / dmax;
708 dval = (ratio * *((gdouble *) hf + k + l*max_x));
709 *((gdouble *)hf + k + l*max_x) = (gdouble) MIN(MAX(dval,0.0),dmax);
710 break;
711 case UNSIGNED_CHAR_ID:
712 val = (glong) (ratio * (gdouble) *((unsigned char *) hf + k + l*max_x));
713 *((unsigned char *)hf + k + l*max_x) = (unsigned char) MIN(MAX(val,0),MAX_HF_VALUE);
714 break;
715 default:
716 val =(glong) (ratio * (gdouble) *((hf_type *) hf + k + l*max_x));
717 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
718 }
719 case MAX_MERGE:
720 switch (hf_data_type) {
721 case GINT_ID:
722 val = MAX(*((gint *) hf + k + l*max_x),val);
723 *((gint *)hf + k + l*max_x) = (gint) MIN(MAX(val,0),MAX_HF_VALUE);
724 break;
725 case HF_TYPE_ID:
726 val = MAX(*((hf_type *) hf + k + l*max_x),val);
727 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
728 break;
729 case GDOUBLE_ID:
730 dval = MAX(*((gdouble *) hf + k + l*max_x),dval);
731 *((gdouble *)hf + k + l*max_x) = (gdouble) MIN(MAX(dval,0.0),dmax);
732 break;
733 case UNSIGNED_CHAR_ID:
734 val = MAX(*((unsigned char *) hf + k + l*max_x),val);
735 *((unsigned char *)hf + k + l*max_x) = (unsigned char) MIN(MAX(val,0),MAX_HF_VALUE);
736 break;
737 default:
738 val = MAX(*((hf_type *) hf + k + l*max_x),val);
739 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
740 }
741 break;
742 case MIN_MERGE:
743 switch (hf_data_type) {
744 case GINT_ID:
745 val = MIN(*((gint *) hf + k + l*max_x),val);
746 *((gint *)hf + k + l*max_x) = (gint) MIN(MAX(val,0),MAX_HF_VALUE);
747 break;
748 case HF_TYPE_ID:
749 val = MIN(*((hf_type *) hf + k + l*max_x),val);
750 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
751 break;
752 case GDOUBLE_ID:
753 dval = MIN(*((gdouble *) hf + k + l*max_x),dval);
754 *((gdouble *)hf + k + l*max_x) = (gdouble) MIN(MAX(dval,0.0),dmax);
755 break;
756 case UNSIGNED_CHAR_ID:
757 val = MIN(*((unsigned char *) hf + k + l*max_x),val);
758 *((unsigned char *)hf + k + l*max_x) = (unsigned char) MIN(MAX(val,0),MAX_HF_VALUE);
759 break;
760 default:
761 val = MIN(*((hf_type *) hf + k + l*max_x),val);
762 *((hf_type *)hf + k + l*max_x) = (hf_type) MIN(MAX(val,0),MAX_HF_VALUE);
763 }
764 break;
765 default:
766 printf(_("Unexpected option in %s; contact the programmer!"),"interpolated_merge");
767 printf("\n");
768 return;
769 }
770 }
771 }
772 }
773
774 // There are some explanations in hf.h!
775
dist_matrix_init(dist_matrix_struct * dist,gint hf_size)776 void dist_matrix_init(dist_matrix_struct *dist, gint hf_size) {
777 // Initializes / reinit the distances matrix with a new hf_size
778 gint i, j;
779 if (dist->hf_size == hf_size) // Nothing to do!
780 return;
781 dist->hf_size = hf_size;
782 dist->size = (hf_size>>1) * (hf_size>>1);
783 if (dist->distances)
784 x_free(dist->distances);
785 dist->distances = (gfloat *) x_malloc(sizeof(gfloat) * dist->size, "gfloat (dist->distances)");
786 for (i=0; i<(hf_size>>1); i++) {
787 for (j=0; j<(hf_size>>1); j++) {
788 (*(dist->distances+VECTORIZE(i,j,hf_size>>1))) = (gfloat) sqrt(i*i + j*j);
789 // printf("I: %d; J: %d; VECTORIZE: %d; DIST: %-10.2f\n",i, j, VECTORIZE(i,j,hf_size>>1), (gfloat) sqrt(i*i + j*j));
790 }
791 }
792 }
793
dist_matrix_new(gint hf_size)794 dist_matrix_struct *dist_matrix_new(gint hf_size) {
795
796 dist_matrix_struct *dist;
797 dist = (dist_matrix_struct *) x_malloc(sizeof(dist_matrix_struct), "dist_matrix_struct");
798 dist->distances = NULL;
799 dist->hf_size = 0;
800 if (hf_size)
801 dist_matrix_init(dist, hf_size);
802 return dist;
803 }
804
805
hf_translate(hf_type * source,hf_type * result,gint max_x,gint max_y,gint dx,gint dy)806 void hf_translate(hf_type *source, hf_type *result,
807 gint max_x, gint max_y, gint dx, gint dy) {
808 gint i,j;
809 for (i=0; i<max_y; i++) {
810 for (j=0; j<max_x; j++) {
811 *(result + VECTORIZE(i,j,max_x)) =
812 *(source + VECTORIZE(WRAP(i+dy,max_y),
813 WRAP2(j+dx,max_x), max_x));
814 }
815 }
816 }
817
hf_slide(hf_struct_type * hf,gint slidev,gint slideh)818 void hf_slide(hf_struct_type *hf, gint slidev, gint slideh) {
819 // Translate a HF vertically or horizontally from 0 to 100 %
820 // Suitable for "wrappable" HF only...
821 gint dx, dy;
822 dx = (slideh * hf->max_x) / 100;
823 dy = - (slidev * hf->max_y) / 100;
824 if (!hf->tmp_buf)
825 hf_backup(hf);
826 hf_translate (hf->tmp_buf, hf->hf_buf, hf->max_x, hf->max_y, dx, dy);
827 }
828
hf_revert(hf_struct_type * hf)829 void hf_revert (hf_struct_type *hf) {
830 gint i;
831 if (hf->tmp_buf) {
832 // We free the temp buffer because we don't need it and we want to avoid
833 // inconsistencies in other actions... "revert" is totally reversible
834 x_free(hf->tmp_buf);
835 hf->tmp_buf = NULL;
836 }
837
838 for (i=0; i<hf->max_x*hf->max_y; i++) {
839 *(hf->hf_buf + i) = ~*(hf->hf_buf +i);
840 }
841 }
842
hf_vertical_mirror(hf_struct_type * hf)843 void hf_vertical_mirror(hf_struct_type *hf) {
844 gint i,j;
845 // We need a fresh temporary buffer here
846 if (!hf->tmp_buf)
847 hf_backup(hf);
848 for (i=0; i<hf->max_y; i++)
849 for (j=0; j<hf->max_x; j++)
850 *(hf->hf_buf+VECTORIZE(i,hf->max_x-j-1, hf->max_x)) =
851 *(hf->tmp_buf+VECTORIZE(i,j, hf->max_x));
852 }
853
hf_horizontal_mirror(hf_struct_type * hf)854 void hf_horizontal_mirror(hf_struct_type *hf) {
855 gint i,j;
856 // We need a fresh temporary buffer here
857 if (!hf->tmp_buf)
858 hf_backup(hf);
859 for (i=0; i<hf->max_y; i++)
860 for (j=0; j<hf->max_x; j++)
861 *(hf->hf_buf+VECTORIZE(hf->max_y-i-1,j, hf->max_x)) =
862 *(hf->tmp_buf+VECTORIZE(i,j, hf->max_x));
863 }
864
interpolate_get_xy(gdouble dbx,gdouble dby,hf_type * hf,gint x_size,gint y_size,gdouble (* get_x)(gdouble,gpointer),gdouble (* get_y)(gdouble,gpointer),gpointer arg_get_x,gpointer arg_get_y)865 hf_type interpolate_get_xy (gdouble dbx, gdouble dby, hf_type *hf,
866 gint x_size, gint y_size,
867 gdouble (*get_x) (gdouble, gpointer), gdouble (*get_y) (gdouble, gpointer),
868 gpointer arg_get_x, gpointer arg_get_y) {
869 // See "interpolate"
870 // "get_x" and "get_y" are functions for getting the coordinates of x and y,
871 // when a special transformation is required (then dbx and dby are not used)
872 // "get_x" and "get_y" could be NULL
873 if (get_x)
874 dbx = (gdouble) (*get_x) (dbx, arg_get_x);
875 if (get_y)
876 dby = (gdouble) (*get_y) (dby, arg_get_y);
877 return interpolate2 (dbx,dby, hf, x_size, y_size);
878 }
879
interpolate2(gdouble dbx,gdouble dby,hf_type * hf,gint x_size,gint y_size)880 hf_type interpolate2 (gdouble dbx, gdouble dby, hf_type *hf,
881 gint x_size, gint y_size) {
882
883 // Find an interpolated value for (dbx,dby) in *hf
884 // Interpolation from distance
885 // Used for rotation and similar functions
886 // "hf" is not required to be square
887
888 // 2004-01-02 Renamed interpolate2, kept for documentation
889 // Use "interpolate" instead
890
891 gdouble dx,dy, v1, v2, v3, v4,d1,d2,d3,d4, value;
892 long int xa, ya, xb, yb;
893 xa = (long int) floor(dbx);
894 ya = (long int) floor(dby);
895 xb = WRAP2(xa ,x_size);
896 yb = WRAP2(ya ,y_size);
897 dx = ABS(dbx - floor(dbx)); // dx,dy = float part, used for interpolation
898 dy = ABS(dby - floor(dby));
899 // printf("(x,y): (%5.2f,%5.2f); (dx,dy): (%f,%f); (xa,ya): (%d,%d); (xb,yb): (%d,%d)\n",dbx,dby,dx,dy,xa,ya,xb,yb);
900 // Simple interpolation from distance
901 d1 = MAX(0.0,1.0 - DIST2(0.0,0.0,dx,dy));
902 d2 = MAX(0.0,1.0 - DIST2(0.0,1.0,dx,dy));
903 d3 = MAX(0.0,1.0 - DIST2(1.0, 0.0, dx, dy));
904 d4 = MAX(0.0,1.0- DIST2(1.0,1.0, dx,dy));
905 /* Too simple bi-linear interpolation... produces some aliasing
906 d1 = 1.0 - dx*dy;
907 d2 = 1.0 - dx * (1.0-dy);
908 d3 = 1.0 - (1.0-dx) * dy;
909 d4 = 1.0 - (1.0-dx) * (1.0-dy);
910 */
911 // printf("(d1,d2,d3,d4): (%f,%f,%f,%f)\n",d1,d2,d3,d4);
912 v1 = d1 * (gdouble) *(hf+VECTORIZE(xb,yb,x_size)) ;
913 v2 = d2 * (gdouble) *(hf+VECTORIZE(xb,WRAP(yb+1,y_size),x_size));
914 v3 = d3 * (gdouble) *(hf+VECTORIZE(WRAP(xb+1,x_size),yb,x_size));
915 v4 = d4 * (gdouble) *(hf+VECTORIZE(WRAP(xb+1,x_size),WRAP(yb+1,y_size),x_size));
916 value = ((v1+v2+v3+v4) / (d1+d2+d3+d4));
917 return (hf_type) value;
918 }
919
translate_real_forward_mapping(gpointer source_grid,gpointer output_grid,gint data_type,gint x_size,gint y_size,gint x,gint y,gdouble ox,gdouble oy)920 void translate_real_forward_mapping (gpointer source_grid,
921 gpointer output_grid,
922 gint data_type,
923 gint x_size, gint y_size,
924 gint x,gint y,
925 gdouble ox, gdouble oy) {
926
927 // Translate the value source_grid(x,y) to output_grid(ox,oy),
928 // given that (ox,oy) are real values, using forward mapping
929 // Opposite of interpolating...
930 // We "spread" the value over the 4 pixels surrounding the target (ox,oy),
931 // using the distance between (ox,oy) and the integer coordinates
932 // We assume that overflow / tiling control is not relevant
933
934 gdouble dx, dy, d1,d2,d3,d4, tot, value;
935 gint ix, iy,x1,y1, x2, y2, i, j1, j2, j3, j4;
936 gboolean t1=FALSE,t2=FALSE,t3=FALSE,t4=FALSE;
937
938 ix = (gint) floor(ox);
939 iy = (gint) floor(oy);
940 x1 = x + ix;
941 y1 = y + iy;
942 if ((x1<-1) || (y1<-1) || (x1>x_size) || (y1>y_size)) // Target trivially out of boundaries
943 return;
944 x2 = x1+1;
945 y2 = y1+1;
946 dx = ox - floor(ox);
947 dy = oy - floor(oy);
948
949 d1 = MAX(0.0,1.0 - DIST2(0.0,0.0,dx,dy));
950 d2 = MAX(0.0,1.0 - DIST2(0.0,1.0,dx,dy));
951 d3 = MAX(0.0,1.0 - DIST2(1.0, 0.0, dx, dy));
952 d4 = MAX(0.0,1.0- DIST2(1.0,1.0, dx,dy));
953 tot = d1+d2+d3+d4;
954 // printf("D1: %5.2f; D2: %5.2f; D3: %5.2f; D4: %5.2f; tot: %5.2f; ",d1,d2,d3,d4, tot);
955 i = VECTORIZE(x, y, x_size);
956 j1= VECTORIZE(x1,y1, x_size); // (x1,y1)
957 j2= j1 + x_size; // (x1,y2)
958 j3 = j1 + 1; // (x2,y1)
959 j4 = j2 + 1; // (x2,y2)
960
961 switch (data_type) {
962 case GINT_ID:
963 value = (gdouble) *( ((gint *) source_grid) + i);
964 if ( (x1>=0) && (y1>=0) && (x1<x_size) && (y1<y_size))
965 *( ((gint *) output_grid) + j1) += (gint) (d1 * value / tot);
966 if ( (x1>=0) && (y2>=0) && (x1<x_size) && (y2<y_size))
967 *( ((gint *) output_grid) + j2) += (gint) (d2 * value / tot);
968 if ( (x2>=0) && (y1>=0) && (x2<x_size) && (y1<y_size))
969 *( ((gint *) output_grid) + j3) += (gint) (d3 * value / tot);
970 if ( (x2>=0) && (y2>=0) && (x2<x_size) && (y2<y_size))
971 *( ((gint *) output_grid) + j4) += (gint) (d4 * value / tot);
972 break;
973 case HF_TYPE_ID:
974 value = (gdouble) *( ((hf_type *) source_grid) + i);
975
976
977 // printf("VALUE (%d, %d): %5.2f; (ix,iy): (%d, %d); j1: %d; j2: %d; j3: %d; j4: %d\n",x,y,value, ix, iy, j1, j2, j3, j4);
978 if ( (x1>=0) && (y1>=0) && (x1<x_size) && (y1<y_size)) {
979 *( ((hf_type *) output_grid) + j1) += (hf_type) (d1 * value / tot);
980 t1 = TRUE;
981 }
982 if ( (x1>=0) && (y2>=0) && (x1<x_size) && (y2<y_size)) {
983 *( ((hf_type *) output_grid) + j2) += (hf_type) (d2 * value / tot);
984 t2 = TRUE;
985 }
986 if ( (x2>=0) && (y1>=0) && (x2<x_size) && (y1<y_size)) {
987 *( ((hf_type *) output_grid) + j3) += (hf_type) (d3 * value / tot);
988 t3 = TRUE;
989 }
990 if ( (x2>=0) && (y2>=0) && (x2<x_size) && (y2<y_size)) {
991 *( ((hf_type *) output_grid) + j4) += (hf_type) (d4 * value / tot);
992 t4 = TRUE;
993 }
994 // if ( (x1<1) || (x1>=(x_size-1)) || (y1<1) || (y1>=(y_size-1)))
995 // printf ("VALUE (%d, %d): %5.2f; (x1,y1): (%d, %d); (x2,y2): (%d, %d); t1,t2,t3,t4: %d,%d,%d,%d; j1: %d; j2: %d; j3: %d; j4: %d; maxi: %d\n",x,y,value,x1,y1,x2,y2,t1,t2,t3,t4,j1,j2,j3,j4,VECTORIZE(x_size-1,y_size-1,x_size));
996 break;
997 case GDOUBLE_ID:
998 value = *( ((gdouble *) source_grid) + i);
999 if ( (x1>=0) && (y1>=0) && (x1<x_size) && (y1<y_size))
1000 *( ((gdouble *) output_grid) + j1) += (gdouble) (d1 * value / tot);
1001 if ( (x1>=0) && (y2>=0) && (x1<x_size) && (y2<y_size))
1002 *( ((gdouble *) output_grid) + j2) += (gdouble) (d2 * value / tot);
1003 if ( (x2>=0) && (y1>=0) && (x2<x_size) && (y1<y_size))
1004 *( ((gdouble *) output_grid) + j3) += (gdouble) (d3 * value / tot);
1005 if ( (x2>=0) && (y2>=0) && (x2<x_size) && (y2<y_size))
1006 *( ((gdouble *) output_grid) + j4) += (gdouble) (d4 * value / tot);
1007
1008 break;
1009 case UNSIGNED_CHAR_ID:
1010 value = *( ((unsigned char *) source_grid) + i);
1011 if ( (x1>=0) && (y1>=0) && (x1<x_size) && (y1<y_size))
1012 *( ((unsigned char *) output_grid) + j1) += (unsigned char) (d1 * value / tot);
1013 if ( (x1>=0) && (y2>=0) && (x1<x_size) && (y2<y_size))
1014 *( ((unsigned char *) output_grid) + j2) += (unsigned char) (d2 * value / tot);
1015 if ( (x2>=0) && (y1>=0) && (x2<x_size) && (y1<y_size))
1016 *( ((unsigned char *) output_grid) + j3) += (unsigned char) (d3 * value / tot);
1017 if ( (x2>=0) && (y2>=0) && (x2<x_size) && (y2<y_size))
1018 *( ((unsigned char *) output_grid) + j4) += (unsigned char) (d4 * value / tot);
1019 break;
1020 default:
1021 printf(_("Unexpected option in %s; contact the programmer!"),
1022 "translate_real_forward_mapping");
1023 printf("\n");
1024 }
1025 }
1026
translate_forward_mapping(gpointer source_grid,gpointer output_grid,gint data_type,gint x_size,gint y_size,gint x,gint y,gint ix,gint iy)1027 void translate_forward_mapping (gpointer source_grid,
1028 gpointer output_grid,
1029 gint data_type,
1030 gint x_size, gint y_size,
1031 gint x,gint y,
1032 gint ix, gint iy) {
1033
1034 // Translate the value source_grid(x,y) to output_grid(ix,iy),
1035 // given that (ix,iy) are integer values, using forward mapping
1036 // We assume that overflow / tiling control is not relevant
1037 // Typically used instead of translate_real_forward_mapping, in preview modes
1038
1039 gint x1,y1, i, j1;
1040
1041 x1 = x + ix;
1042 y1 = y + iy;
1043 if ((x1<0) || (y1<0) || (x1>=x_size) || (y1>=y_size)) // Target out of boundaries
1044 return;
1045
1046 i = VECTORIZE(x, y, x_size);
1047 j1= VECTORIZE(x1,y1, x_size);
1048
1049 switch (data_type) {
1050 case GINT_ID:
1051 *( ((gint *) output_grid) + j1) = *( ((gint *) source_grid) + i);
1052 break;
1053 case HF_TYPE_ID:
1054 *( ((hf_type *) output_grid) + j1) = *( ((hf_type *) source_grid) + i);
1055 break;
1056 case GDOUBLE_ID:
1057 *( ((gdouble *) output_grid) + j1) = *( ((gdouble *) source_grid) + i);
1058 break;
1059 case UNSIGNED_CHAR_ID:
1060 *( ((unsigned char *) output_grid) + j1) = *( ((unsigned char *) source_grid) + i);
1061 break;
1062 default:
1063 printf(_("Unexpected option in %s; contact the programmer!"),
1064 "translate_forward_mapping");
1065 printf("\n");
1066 }
1067 }
interpolate(gdouble dbx,gdouble dby,gpointer grid,gint x_size,gint y_size,gpointer return_value_ptr,gint data_type,gint overflow)1068 void interpolate (gdouble dbx, gdouble dby, gpointer grid,
1069 gint x_size, gint y_size, gpointer return_value_ptr, gint data_type,
1070 gint overflow) {
1071
1072 // Find an interpolated value for (dbx,dby) in "grid"
1073 // Interpolation from distance
1074 // Used for rotation and similar functions
1075 // (dbx,dby): real coordinates in "grid"
1076 // grid: any table of type "data_type"; not required to be square
1077 // x_size * y_size: size of the grid, in "data_type" units
1078 // return_value_ptr: pointer in which the value is returned - allocated by the calling function
1079 // data_type: any of GINT_ID, HF_TYPE_ID,
1080 // GDOUBLE_ID, UNSIGNED_CHAR_ID (defined in hf.h)
1081 // overflow: any of OVERFLOW_WRAP, ..._REBOUND, ..._ZERO, ..._IDLE
1082
1083 gdouble dx,dy, v1=0.0, v2=0.0, v3=0.0, v4=0.0,d1,d2,d3,d4;
1084 gboolean tx1=TRUE,tx2=TRUE,ty1=TRUE,ty2=TRUE;
1085 glong xa, ya, i1, i2, i3, i4, x1, x2, y1, y2;
1086
1087 // printf("INTERPOLATING (%5.2f,%5.2f); ", dbx,dby);
1088 xa = (glong) floor(dbx);
1089 ya = (glong) floor(dby);
1090 dx = ABS(dbx - floor(dbx)); // dx,dy = float part, used for interpolation
1091 dy = ABS(dby - floor(dby));
1092
1093 // Simple interpolation from distance
1094 d1 = MAX(0.0,1.0 - DIST2(0.0,0.0,dx,dy));
1095 d2 = MAX(0.0,1.0 - DIST2(0.0,1.0,dx,dy));
1096 d3 = MAX(0.0,1.0 - DIST2(1.0, 0.0, dx, dy));
1097 d4 = MAX(0.0,1.0- DIST2(1.0,1.0, dx,dy));
1098
1099 // printf("(d1,d2,d3,d4): (%f,%f,%f,%f)\n",d1,d2,d3,d4);
1100 switch (overflow) {
1101 case OVERFLOW_WRAP:
1102 x1 = WRAP2(xa,x_size);
1103 x2 = WRAP2(xa+1,x_size);
1104 y1 = WRAP2(ya,y_size);
1105 y2 = WRAP2(ya+1,y_size);
1106 break;
1107 case OVERFLOW_REBOUND:
1108 x1 = REBOUND(xa,x_size);
1109 x2 = REBOUND(xa+1,x_size);
1110 y1 = REBOUND(ya,y_size);
1111 y2 = REBOUND(ya+1,y_size);
1112 break;
1113 case OVERFLOW_IDLE:
1114 x1 = MAX(0,MIN(xa,x_size-1));
1115 x2 = MAX(0,MIN(xa+1,x_size-1));
1116 y1 = MAX(0,MIN(ya,y_size-1));
1117 y2 = MAX(0,MIN(ya+1,y_size-1));
1118 break;
1119 case OVERFLOW_ZERO:
1120 x1 = xa;
1121 x2 = xa+1;
1122 y1 = ya;
1123 y2 = ya+1;
1124 if ( (x1<0) || (x1>=x_size) ) {
1125 tx1 = FALSE;
1126 x1 = MAX(0,MIN(x1,x_size-1)); // We initialize to a valid index, for avoiding errors
1127 }
1128 if ( (x2>=x_size) || (x2<0) ) {
1129 tx2 = FALSE;
1130 x2 = MAX(0,MIN(x2,x_size-1));
1131 }
1132 if ( (y1<0) || (y1>=y_size) ) {
1133 ty1 = FALSE;
1134 y1 = MAX(0,MIN(y1,y_size-1));
1135 }
1136 if ( (y2>=y_size) || (y2<0) ) {
1137 ty2 = FALSE;
1138 y2 = MAX(0,MIN(y2,y_size-1));
1139 }
1140 break;
1141 default:
1142 printf(_("Unexpected option in %s; contact the programmer!"),"interpolate");
1143 printf("\n");
1144 return;
1145 }
1146
1147 i1 = VECTORIZE(x1,y1,x_size);
1148 i2 = VECTORIZE(x1,y2,x_size);
1149 i3 = VECTORIZE(x2,y1,x_size);
1150 i4 = VECTORIZE(x2,y2,x_size);
1151
1152 switch (data_type) {
1153 case GINT_ID:
1154 if (tx1 && ty1)
1155 v1 = d1 * (gdouble) *( ((gint*) grid) + i1) ;
1156 if (tx1 && ty2)
1157 v2 = d2 * (gdouble) *( ((gint*) grid) + i2);
1158 if (tx2 && ty1)
1159 v3 = d3 * (gdouble) *( ((gint*) grid) + i3);
1160 if (tx2 && ty2)
1161 v4 = d4 * (gdouble) *( ((gint*) grid) + i4);
1162 *((gint *) return_value_ptr) = (gint) ((v1+v2+v3+v4) / (d1+d2+d3+d4)) ;
1163 break;
1164 case HF_TYPE_ID:
1165 if (tx1 && ty1)
1166 v1 = d1 * (gdouble) *( ((hf_type*) grid) + i1) ;
1167 if (tx1 && ty2)
1168 v2 = d2 * (gdouble) *( ((hf_type*) grid) + i2);
1169 if (tx2 && ty1)
1170 v3 = d3 * (gdouble) *( ((hf_type*) grid) + i3);
1171 if (tx2 && ty2)
1172 v4 = d4 * (gdouble) *( ((hf_type*) grid) + i4);
1173 *((hf_type *) return_value_ptr) = (hf_type) ((v1+v2+v3+v4) / (d1+d2+d3+d4)) ;
1174 // printf("VALUE: %d; (x1,y1): (%d,%d); (x2,y2): (%d,%d); d1: %5.2f; d2: %5.2f; d3: %5.2f; d4: %5.2f; v1: %5.2f; v2: %5.2f; v3: %5.2f; v4: %5.2f\n",* (hf_type*) return_value_ptr,x1,y1,x2,y2,d1,d2,d3,d4,v1,v2,v3,v4);
1175 break;
1176 case GDOUBLE_ID:
1177 if (tx1 && ty1)
1178 v1 = d1 * (gdouble) *( ((gdouble*) grid) + i1) ;
1179 if (tx1 && ty2)
1180 v2 = d2 * (gdouble) *( ((gdouble*) grid) + i2);
1181 if (tx2 && ty1)
1182 v3 = d3 * (gdouble) *( ((gdouble*) grid) + i3);
1183 if (tx2 && ty2)
1184 v4 = d4 * (gdouble) *( ((gdouble*) grid) + i4);
1185 *((gdouble *) return_value_ptr) = (gdouble) ((v1+v2+v3+v4) / (d1+d2+d3+d4)) ;
1186 break;
1187 case UNSIGNED_CHAR_ID:
1188 if (tx1 && ty1)
1189 v1 = d1 * (gdouble) *( ((unsigned char*) grid) + i1) ;
1190 if (tx1 && ty2)
1191 v2 = d2 * (gdouble) *( ((unsigned char*) grid) + i2);
1192 if (tx2 && ty1)
1193 v3 = d3 * (gdouble) *( ((unsigned char*) grid) + i3);
1194 if (tx2 && ty2)
1195 v4 = d4 * (gdouble) *( ((unsigned char*) grid) + i4);
1196 *((unsigned char *) return_value_ptr) = (unsigned char) ((v1+v2+v3+v4) / (d1+d2+d3+d4)) ;
1197 break;
1198 default:
1199 printf(_("Unexpected option in %s; contact the programmer!"),"interpolate");
1200 printf("\n");
1201 }
1202 }
1203
vector_interpolate(gdouble dbx,gpointer vector,gint x_size,gpointer return_value_ptr,gint data_type,gint overflow)1204 void vector_interpolate (gdouble dbx, gpointer vector, gint x_size,
1205 gpointer return_value_ptr, gint data_type, gint overflow) {
1206 // Same as "interpolate", for a vector instead of a grid
1207 glong xa, x1, x2;
1208 gdouble fx, d1, d2, v1=0.0, v2=0.0;
1209 gboolean tx1=TRUE, tx2=TRUE;
1210
1211 fx = floor(dbx);
1212 xa = (glong) fx;
1213 d2 = ABS(dbx - fx);
1214 d1 = 1.0 - d2;
1215
1216 switch (overflow) {
1217 case OVERFLOW_WRAP:
1218 x1 = WRAP2(xa,x_size);
1219 x2 = WRAP2(xa+1,x_size);
1220 break;
1221 case OVERFLOW_REBOUND:
1222 x1 = REBOUND(xa,x_size);
1223 x2 = REBOUND(xa+1,x_size);
1224 break;
1225 case OVERFLOW_IDLE:
1226 x1 = MAX(0,MIN(xa,x_size));
1227 x2 = MAX(0,MIN(xa+1,x_size));
1228 break;
1229 case OVERFLOW_ZERO:
1230 x1 = xa;
1231 x2 = xa+1;
1232 if ( (x1<0) || (x1>=x_size) ) {
1233 tx1 = FALSE;
1234 x1 = MAX(0,MIN(x1,x_size)); // We initialize to a valid index, for avoiding errors
1235 }
1236 if ( (x2>=x_size) || (x2<0) ) {
1237 tx2 = FALSE;
1238 x2 = MAX(0,MIN(x2,x_size));
1239 }
1240 break;
1241 default:
1242 printf(_("Unexpected option in %s; contact the programmer!"),"vector_interpolate");
1243 printf("\n");
1244 return;
1245 }
1246 switch (data_type) {
1247 case GINT_ID:
1248 if (tx1)
1249 v1 = d1 * (gdouble) *( ((gint*) vector) + x1) ;
1250 if (tx2)
1251 v2 = d2 * (gdouble) *( ((gint*) vector) + x2);
1252 *((gint *) return_value_ptr) = (gint) ((v1+v2) / (d1+d2)) ;
1253 break;
1254 case HF_TYPE_ID:
1255 if (tx1)
1256 v1 = d1 * (gdouble) *( ((hf_type*) vector) + x1) ;
1257 if (tx2)
1258 v2 = d2 * (gdouble) *( ((hf_type*) vector) + x2);
1259 *((hf_type *) return_value_ptr) = (hf_type) ((v1+v2) / (d1+d2)) ;
1260
1261 // printf("VALUE: %d; x1: %d; x2: %d; d1: %5.2f; d2: %5.2f; v1: %5.2f; v2: %5.2f;\n",* (hf_type*) return_value_ptr,x1,x2,d1,d2,v1,v2);
1262 break;
1263 case GDOUBLE_ID:
1264 if (tx1)
1265 v1 = d1 * (gdouble) *( ((gdouble*) vector) + x1) ;
1266 if (tx2)
1267 v2 = d2 * (gdouble) *( ((gdouble*) vector) + x2);
1268 *((gdouble *) return_value_ptr) = (gdouble) ((v1+v2) / (d1+d2)) ;
1269 // printf("VALUE: 5.2f; dbx: %5.2f; x1: %d; x2: %d; d1: %5.2f; d2: %5.2f; v1: %5.2f; v2: %5.2f;\n",* (gdouble*) return_value_ptr,dbx,x1,x2,d1,d2,v1,v2);
1270 break;
1271 case UNSIGNED_CHAR_ID:
1272 if (tx1)
1273 v1 = d1 * (gdouble) *( ((unsigned char*) vector) + x1) ;
1274 if (tx2)
1275 v2 = d2 * (gdouble) *( ((unsigned char*) vector) + x2);
1276 *((unsigned char *) return_value_ptr) = (unsigned char) ((v1+v2) / (d1+d2)) ;
1277 break;
1278 default:
1279 printf(_("Unexpected option in %s; contact the programmer!"),"vector_interpolate");
1280 printf("\n");
1281 }
1282 }
1283
hf_circular_projection(hf_type * hf_in,hf_type * hf_out,gint hf_size)1284 void hf_circular_projection (hf_type *hf_in, hf_type *hf_out, gint hf_size) {
1285 // Compresses the corners of the HF so that it fits into a circle
1286 // Completes the corner by wrapping the coordinates
1287 // The HF should be tileable for a seamless result
1288 // 2003-11-29
1289 // Intended uses: (1) curiosity; (2) minimize the "square effect" of some transformations
1290 // 2003-11-30 A bit disappointing - won't implement
1291 // We use inverse mapping
1292 gint x,y, half;
1293 gdouble dist, dbx, dby, in_x, in_y;
1294 half = hf_size / 2; // We center the projection
1295 for (y=0; y<hf_size; y++) {
1296 dby = (gdouble) (y-half);
1297 for (x=0; x<hf_size; x++) {
1298 dbx = (gdouble) (x-half);
1299 if ((x==0) || (y==0)) // Nothing to do
1300 *(hf_out+VECTORIZE(x,y,hf_size)) = *(hf_in+VECTORIZE(x,y,hf_size));
1301 else {
1302 dist = DIST2(0.0,0.0, dbx, dby); // Distance from center
1303 if (ABS(dbx)>ABS(dby)) {
1304 in_x = (dbx<0?-1:1)*dist;
1305 in_y = in_x*dby/dbx;
1306 }
1307 else {
1308 in_y = (dby<0?-1:1)*dist;
1309 in_x = in_y*dbx/dby;
1310 } // if (test abs) ... else
1311
1312 // printf("(x,y): (%d,%d); (dbx,dby): (%5.1f,%5.1f); dist: %7.3f; (in_x,in_y): (%d,%d)\n",x,y,dbx,dby,dist,half +(gint) in_x,half + (gint) in_y);
1313 *(hf_out+VECTORIZE(x,y,hf_size)) =
1314 *(hf_in+VECTORIZE(WRAP2(half + (gint) in_x,hf_size),
1315 WRAP2(half + (gint) in_y,hf_size),hf_size));
1316 } // if (test 0)...
1317 } // for y
1318 } // for x
1319 }
1320
rotate_sin_cos(gdouble sin,gdouble cos,gpointer map_in,gpointer map_out,gint hsize,gint vsize,gint data_type,gint overflow)1321 void rotate_sin_cos (gdouble sin, gdouble cos, gpointer map_in, gpointer map_out, gint hsize, gint vsize, gint data_type, gint overflow) {
1322
1323 // Rotate map_in in map_out
1324 // Both arrays are of type data_type (GINT_ID, HF_TYPE_ID,
1325 // GDOUBLE_ID, UNSIGNED_CHAR_ID)
1326 // Both of size hsize * vsize
1327 // Sine and cosine are given by dx, dy (no explicit computation of angle)
1328 // map_out is supposed to be allocated to the right size
1329 // Reverse mapping process
1330
1331 gdouble x_rot, y_rot, xx, yy, dradius_y, dradius_x;
1332 gint x, y;
1333 hf_type return_hf_type;
1334 unsigned char return_uns_char;
1335 gdouble return_gdouble;
1336 gint return_gint;
1337 sin = -sin ; // Reverse mapping
1338 dradius_y = 0.5 * ((gdouble) vsize - 1.0);
1339 dradius_x = 0.5 * ((gdouble) hsize - 1.0);
1340
1341 for (y=0; y<vsize; y++) {
1342 for (x=0; x<hsize; x++) {
1343 xx = ((gdouble) x) - dradius_x; // Rotate around the center, not (0,0)
1344 yy = ((gdouble) y) - dradius_y;
1345 x_rot = dradius_x + xx*cos - yy*sin;
1346 y_rot = dradius_y + yy*cos + xx*sin;
1347
1348 switch (data_type) {
1349 case GINT_ID:
1350 interpolate (x_rot, y_rot, map_in, hsize, vsize, &return_gint,
1351 GINT_ID, overflow);
1352 *(((gint *) map_out)+VECTORIZE(x,y,hsize)) = return_gint;
1353 break;
1354 case HF_TYPE_ID:
1355 // The next lines help to test the effect of "interpolate"
1356 // return_hf_type = *(((hf_type *) map_in)+VECTORIZE(WRAP2((gint) x_rot,hsize), WRAP2((gint) y_rot, vsize), hsize));
1357 // return_hf_type = *(((hf_type *) map_in)+
1358 // VECTORIZE(REBOUND((gint) x_rot,hsize), REBOUND((gint) y_rot, vsize), hsize));
1359 interpolate (x_rot, y_rot, map_in, hsize, vsize, &return_hf_type,
1360 HF_TYPE_ID, overflow);
1361 *(((hf_type *) map_out)+VECTORIZE(x,y,hsize)) = return_hf_type;
1362 break;
1363 case GDOUBLE_ID:
1364 interpolate (x_rot, y_rot, map_in, hsize, vsize, &return_gdouble,
1365 GDOUBLE_ID, overflow);
1366 *(((gdouble *) map_out)+VECTORIZE(x,y,hsize)) = return_gdouble;
1367 break;
1368 case UNSIGNED_CHAR_ID:
1369 interpolate (x_rot, y_rot, map_in, hsize, vsize, &return_uns_char,
1370 UNSIGNED_CHAR_ID, overflow);
1371 *(((unsigned char *) map_out)+VECTORIZE(x,y,hsize)) = return_uns_char;
1372 break;
1373 default:
1374 printf(_("Unexpected option in %s; contact the programmer!"),"rotate_sin_cos");
1375 printf("\n");
1376 return;
1377 } // switch
1378 } // for (x)
1379 } // for (y)
1380 }
1381
rotate(gdouble dx,gdouble dy,gpointer map_in,gpointer map_out,gint hsize,gint vsize,gint data_type,gint overflow)1382 void rotate (gdouble dx, gdouble dy, gpointer map_in, gpointer map_out,
1383 gint hsize, gint vsize, gint data_type, gint overflow) {
1384
1385 // Rotate map_in in map_out
1386 // Both arrays are of type data_type (GINT_ID, HF_TYPE_ID,
1387 // GDOUBLE_ID, UNSIGNED_CHAR_ID)
1388 // Both of size hsize * vsize
1389 // Sine and cosine are given by dx, dy (no explicit computation of angle)
1390 // map_out is supposed to be allocated to the right size
1391 // Reverse mapping process
1392
1393 gdouble ddist, cos, sin;
1394 ddist = sqrt( pow( dx ,2.0) + pow(dy, 2.0) ) ;
1395 cos = dx / ddist;
1396 sin = dy / ddist;
1397 rotate_sin_cos (sin, cos, map_in, map_out, hsize, vsize, data_type, overflow);
1398 }
1399
hf_rotate(hf_type * hf_in,hf_type * hf_out,gint hf_size,gint angle,gboolean wrap,gint overflow)1400 void hf_rotate (hf_type *hf_in, hf_type *hf_out, gint hf_size,
1401 gint angle, gboolean wrap, gint overflow) {
1402 // Rotation of a square HF
1403 // Angle is in degrees
1404 gdouble a;
1405 a = PI * ((gdouble) angle) / 180.0;
1406 rotate_sin_cos (sin(a), cos(a), (gpointer) hf_in, (gpointer) hf_out, hf_size, hf_size, HF_TYPE_ID,
1407 overflow);
1408 }
1409
hf_fast_rotate(hf_type * hf_in,hf_type * hf_out,gint hf_size,gint angle)1410 void hf_fast_rotate (hf_type *hf_in, hf_type *hf_out, gint hf_size, gint angle) {
1411 // Rotation of a square HF, preview style, with aliasing
1412 // "Wraps"
1413 // Angle is in degrees
1414 gdouble a, s, c, x_rot, y_rot;
1415 long int xa, ya, x, y;
1416 a = PI * ((gdouble) angle) / 180;
1417 // Reverse mapping
1418 s = sin(-a);
1419 c = cos(-a);
1420 for (y=0; y<hf_size; y++)
1421 for (x=0; x<hf_size; x++) {
1422 // Nearest neighbor interpolation
1423 x_rot = x*c - y*s ;
1424 y_rot = y*c + x*s ;
1425 xa = WRAP2((long int) (x_rot+0.5),hf_size);
1426 ya = WRAP2((long int) (y_rot+0.5),hf_size);
1427 *(hf_out+VECTORIZE(x,y,hf_size)) = *(hf_in+VECTORIZE(xa,ya,hf_size)) ;
1428 }
1429 }
1430
norm(vector * v)1431 void norm(vector *v) {
1432 // Normalizes a vector, so that it's length is 1
1433 v->l = sqrt(v->x*v->x + v->y*v->y + v->z*v->z);
1434 v->x = v->x / v->l;
1435 v->y = v->y / v->l;
1436 v->z = v->z / v->l;
1437 }
1438
normalized_bell_new(gint radius)1439 gdouble* normalized_bell_new(gint radius) {
1440 // Allocates and initializes half a gaussian bell (plus the 0 value)
1441 // It's normalized, so that the sum of the whole bell is 1.0
1442 gdouble *result, sum=0.0, value;
1443 gint i;
1444 result = (gdouble *) x_malloc(sizeof(gdouble)* (radius+1), "gdouble (result in normalized_bell_new)");
1445 if (!radius) {
1446 *result = 1.0;
1447 return result;
1448 }
1449 *result = BELLD(CONST_E,2.0,1,radius);
1450 sum = *result;
1451 for (i=1; i<=radius; i++) {
1452 value = BELLD(CONST_E,2.0,i+1,radius);
1453 *(result+i) = value;
1454 sum+=2*value; // once for the positive part, once for the negative
1455 }
1456 // Normalize
1457 value = 1.0 / sum;
1458 for (i=0; i<=radius; i++) {
1459 *(result+i) = *(result+i) * value;
1460 }
1461 return result;
1462 }
1463
convolve_normalized_vector(gpointer in,gpointer out,gint max_x,gint max_y,gboolean wrap,gint radius,gdouble * vector,gint data_type)1464 void convolve_normalized_vector (gpointer in,gpointer out, gint max_x, gint max_y, gboolean wrap,
1465 gint radius, gdouble *vector, gint data_type) {
1466 // Convolve with a normalized "half-vector" on x and y
1467 // (the vector is symmetric around 0 and the sum of the whole is 1)
1468 // "Separated" algorithm
1469 // Actually it's the same as convolving with a matrix
1470 // equal to the cross product of the vector
1471 // This function can take a rectangular image
1472 // Accepts 4 data types:
1473 // data_type = GINT_ID for gint / int
1474 // data_type = HF_TYPE_ID for hf_type (unsigned short int)
1475 // data_type = GDOUBLE_ID for gdouble / double
1476 // data_type = UNSIGNED_CHAR_ID for unsigned char
1477 gdouble *buf, value, *d_in, *d_out;
1478 gint x,y,i, idx, idx_in,max2x, max2y, *i_in, *i_out;
1479 hf_type *hf_in, *hf_out;
1480 unsigned char *us_in, *us_out;
1481 i_in = in; i_out = out;
1482 d_in = in; d_out = out;
1483 hf_in = in; hf_out = out;
1484 us_in = in; us_out = out;
1485 max2x = 2*max_x-1;
1486 max2y = 2*max_y-1;
1487 buf = (gdouble *) x_malloc(sizeof(gdouble)*max_x*max_y, "gdouble (buf in convolve_normalized_vector)");
1488 // Summarize on the x axis
1489 for (y=0; y<max_y; y++) {
1490 for (x=0; x<max_x; x++) {
1491 idx = VECTORIZE(x,y,max_x);
1492 *(buf+idx) = 0.0;
1493 for (i=-radius; i<=radius; i++) {
1494 if (wrap)
1495 idx_in = VECTORIZE(WRAP2(x+i,max_x),y, max_x);
1496 else { // "rebounds"
1497 idx_in = ABS(x+i);
1498 idx_in = (idx_in>=max_x)?(max2x-idx_in):idx_in;
1499 idx_in = VECTORIZE( idx_in, y, max_x);
1500 }
1501 switch (data_type) {
1502 case GINT_ID:
1503 *(buf+idx) += *(vector+ABS(i)) * (gdouble) *(i_in+idx_in);
1504 break;
1505 case HF_TYPE_ID:
1506 *(buf+idx) += *(vector+ABS(i)) * (gdouble) *(hf_in+idx_in);
1507 break;
1508 case GDOUBLE_ID:
1509 *(buf+idx) += *(vector+ABS(i)) * (gdouble) *(d_in+idx_in);
1510 break;
1511 case UNSIGNED_CHAR_ID:
1512 *(buf+idx) += *(vector+ABS(i)) * (gdouble) *(us_in+idx_in);
1513 break;
1514 default:
1515 *(buf+idx) += *(vector+ABS(i)) * (gdouble) *(d_in+idx_in);
1516 }
1517 }
1518 }
1519 }
1520 // Summarize on the y axis
1521 for (y=0; y<max_y; y++) {
1522 for (x=0; x<max_x; x++) {
1523 idx = VECTORIZE(x,y,max_x);
1524 value = 0.0 ;
1525
1526 for (i=-radius; i<=radius; i++) {
1527 if (wrap)
1528 idx_in = VECTORIZE(x,WRAP2(y+i,max_y), max_y);
1529 else { // "rebounds"
1530 idx_in = ABS(y+i);
1531 idx_in = (idx_in>=max_y)?(max2y-idx_in):idx_in;
1532 idx_in = VECTORIZE( x, idx_in, max_y);
1533 }
1534 value += *(vector+ABS(i)) * *(buf+idx_in);
1535 }
1536 switch (data_type) {
1537 case GINT_ID:
1538 *(i_out+idx) = (gint) value;
1539 break;
1540 case HF_TYPE_ID:
1541 *(hf_out+idx) = (hf_type) value;
1542 break;
1543 case GDOUBLE_ID:
1544 *(d_out+idx) = (gdouble) value;
1545 break;
1546 case UNSIGNED_CHAR_ID:
1547 *(us_out+idx) = (unsigned char) value;
1548 break;
1549 default:
1550 *(d_out+idx) = value;
1551 }
1552 }
1553 }
1554
1555 x_free(buf);
1556 }
1557
map_convolve(hf_type * map,gint map_max_x,gint map_max_y,hf_type * background,gint max_x,gint max_y,gint cx,gint cy,gboolean wrap,gint level,gdouble ** gauss_list,gboolean square_symmetry)1558 void map_convolve (hf_type *map, gint map_max_x, gint map_max_y,
1559 hf_type *background, gint max_x, gint max_y, gint cx, gint cy,
1560 gboolean wrap, gint level, gdouble **gauss_list, gboolean square_symmetry) {
1561
1562 // Smoothing (normalized convolution) of a square subset of the image "background"
1563 // The map is typically a stroke dot from the pen
1564 // whose values are used as relative smoothing radii
1565 // The map is centered at (cx,cy) in the background
1566 // The map and the background don't have to be square
1567 // Level is % of maximum radius
1568 // The result is "local smoothing"
1569 // The gaussian vectors generated are cached in gauss_list, for performance reasons
1570
1571 gint x,y,i, idxbuf, idxmap, idxin, idxout, max2x, max2y,
1572 sx,sy, radius, radius_x, radius_y, mx, my, lx, ix, iy;
1573 guint buf_HEAP_index=0xFFFF, background_HEAP_index=0xFFFF;
1574 gdouble *buf,value,*vector;
1575 gfloat ratio; // Ratio used to transform map values to radii
1576 ratio = (gfloat) ( ((gfloat) level) * ((gfloat) GAUSS_LIST_LENGTH) / ((gfloat) MAX_HF_VALUE) / 100.0);
1577 // printf("MAP_CONVOLVE_RATIO: %8.6f\n",ratio);
1578 // We prepare a buffer of the region to smooth
1579 // The buffer must range from (x or y)-(map_size>>1) to (x or y)+(map_size>>1)
1580 buf = (gdouble *) x_malloc(sizeof(gdouble)*map_max_x*map_max_y, "gdouble (buf in map_convolve)");
1581 max2x = 2*map_max_x-1;
1582 max2y = 2*map_max_y-1;
1583 // Control of the square symmetry (when the map contains only one quadrant)
1584 // lx: actual length of one line of the map
1585 radius_x = (map_max_x>>1);
1586 radius_y = (map_max_y>>1);
1587 if (square_symmetry)
1588 lx = radius_x + 1;
1589 else
1590 lx = map_max_x;
1591 // Start point in the background
1592 sy = cy - (map_max_y>>1);
1593 sx = cx - (map_max_x>>1);
1594 // Summarize along the x axis
1595 // (x,y) are indexes in the map
1596 // From the map standpoint, the convolution matrix always "rebounds" on the boundaries
1597 // (the input - index is "idxin")
1598 // From the background standpoint, the map wraps over the boundaries
1599 // when "wrap" is TRUE, or the map overflow is lost in a blackhole if "wrap" is FALSE
1600 // (the output - index is "idxout")
1601
1602 for (y=0; y<map_max_y; y++) {
1603 for (x=0; x<map_max_x; x++) {
1604 idxbuf = VECTORIZE(x,y,map_max_x);
1605 *(buf+idxbuf) = 0.0;
1606 if (square_symmetry) { // Only one quadrant is defined in the map
1607 if (x>radius_x)
1608 ix = 2* radius_x - x;
1609 else
1610 ix = x;
1611 if (y>radius_y)
1612 iy = 2* radius_y - y;
1613 else
1614 iy = y;
1615 idxmap = VECTORIZE(ix,iy,lx);
1616 }
1617 else
1618 idxmap = idxbuf;
1619 // Find the current radius
1620 radius = (gint) (ratio * *(map+idxmap));
1621 // if (!radius)
1622 // continue;
1623 // printf("RADIUS for (%d,%d,%d): %d\n",x,y,*(map+idxmap),radius);
1624 if (!radius) {
1625 *(buf+idxbuf) = (gdouble)*(background+VECTORIZE(WRAP2(sx+ABS(x),max_x),
1626 WRAP2(sy+y,max_y),max_x));
1627 continue;
1628 }
1629 if (!gauss_list[radius])
1630 gauss_list[radius] = normalized_bell_new(radius);
1631 vector = gauss_list[radius];
1632 for (i=-radius; i<=radius; i++) {
1633 // Find the index of the value to add at this step
1634 // The index is always in the "window" delimited by the map
1635 // Overflows rebound
1636 // (mx,y) is the relative pixel to use in the map "window"
1637 // (sx+mx, sy+y) would be the absolute address of the pixel, in the background
1638 // In this set of loops, we run along the X axis, so Y won't overflow here
1639 mx = ABS(x+i);
1640 mx = (mx>=map_max_x)?(max2x-mx):mx;
1641 idxin = VECTORIZE(WRAP2(sx+mx,max_x),WRAP2(sy+y,max_y),max_x);
1642 *(buf+idxbuf) += *(vector+ABS(i)) * (gdouble) *(background+idxin);
1643 }
1644 }
1645 }
1646
1647 // Summarize along the y axis
1648 for (y=0; y<map_max_y; y++) {
1649 for (x=0; x<map_max_x; x++) {
1650 if (square_symmetry) { // Only one quadrant is defined in the map
1651 if (x>radius_x)
1652 ix = 2* radius_x - x;
1653 else
1654 ix = x;
1655 if (y>radius_y)
1656 iy = 2* radius_y - y;
1657 else
1658 iy = y;
1659 idxmap = VECTORIZE(ix,iy,lx);
1660 }
1661 else
1662 idxmap = VECTORIZE(x,y,map_max_x);
1663 value = 0.0 ;
1664 // Find the current radius
1665 radius = (gint) (ratio * *(map+idxmap));
1666 if (!radius) {
1667 value = *(buf+VECTORIZE(x,y,map_max_y));
1668 }
1669 else {
1670 // Obviously we don't need to check the existence of the vector for this radius!
1671 vector = gauss_list[radius];
1672 for (i=-radius; i<=radius; i++) {
1673 // (x,my) is the relative pixel to use in the map "window"
1674 // (sx+x, sy+my) would be the absolute address of the pixel, in the background
1675 my = ABS(y+i);
1676 my = (my>=map_max_y)?(max2y-my):my;
1677 // We take the buffer as input now
1678 idxin = VECTORIZE(x,my, map_max_y);
1679 value += *(vector+ABS(i)) * *(buf+idxin);
1680 }
1681 }
1682 // Here we drop outbound values in the nothingness if wrap is FALSE
1683 if (wrap)
1684 idxout = VECTORIZE(WRAP2(sx+x,max_x),WRAP2(sy+y, max_y), max_y);
1685 else {
1686 if ( ((sx+x)<0) || ((sy+y)<0) || ((sx+x)>=max_x) || ((sy+y)>=max_y) )
1687 continue;
1688 else
1689 idxout = VECTORIZE(sx+x,sy+y,max_y);
1690 }
1691 *(background+idxout) = (hf_type) value;
1692 } // end for(x...)
1693 } // end for(y...)
1694 x_free(buf);
1695 }
1696
intersect_rect(gint xmin,gint ymin,gint xmax,gint ymax,gint x0,gint y0,gint x1,gint y1,gdouble * startx,gdouble * starty,gdouble * endx,gdouble * endy)1697 gboolean intersect_rect (gint xmin, gint ymin, gint xmax, gint ymax,
1698 gint x0, gint y0, gint x1, gint y1,
1699 gdouble *startx, gdouble *starty, gdouble *endx, gdouble *endy) {
1700
1701 // Extend the line (x0,y0) -> (x1,y1) so that it intersects
1702 // the rectangle (xmin,ymin) -> (xmax,ymax)
1703 // Put the result in (startx, starty) and (endx, endy)
1704 // Adapted from Liang-Barsky clipping line algorithm (Foley - van Dam)
1705
1706 // IMPORTANT: xmax and ymax are not sizes, but actual indexes in base 0
1707 gdouble t[4], tn, tp, dx, dy;
1708 gint i=3;
1709 static gdouble epsilon=0.1,MAXPOS=1000000.0, MAXNEG=-1000000.0 ;
1710
1711 tn = MAXNEG;
1712 tp = MAXPOS;
1713
1714 // P0 and P1 must be in the rectangle
1715 if ( (x0<0) || (x0>xmax) || (y0<0) || (y0>ymax) ||
1716 (x1<0) || (x1>xmax) || (y1<0) || (y1>ymax) )
1717 return FALSE;
1718 // printf ("RECT: (%d,%d) - (%d,%d); LINE: (%d,%d) - (%d,%d); \n",
1719 // xmin,ymin,xmax,ymax, x0, y0, x1, y1);
1720
1721 dx = (gdouble) (x1-x0);
1722 dy = (gdouble) (y1-y0);
1723
1724 if (ABS(dx)<epsilon) {
1725 t[0] = MAXNEG;
1726 t[1] = MAXPOS;
1727 }
1728 else {
1729 t[0] = - ((gdouble) (x0 - xmin)) / dx; // Left
1730 t[1] = ((gdouble) (x0 - xmax)) / -dx; // Right
1731 }
1732 if (ABS(dy)<epsilon) {
1733 t[2] = MAXNEG;
1734 t[3] = MAXPOS;
1735 }
1736 else {
1737 t[2] = - ((gdouble) (y0-ymin)) / dy; // Bottom
1738 t[3] = ((gdouble) (y0-ymax)) / -dy; // Top
1739 }
1740
1741 // printf("T[0]: %5.2f; T[1]: %5.2f; T[2]: %5.2f; T[3]: %5.2f; \n",t[0],t[1],t[2],t[3]);
1742 // We retain the smallest negative T and the smallest positive (nearest from P0)
1743 // then we compute the points from the parametric formula of P0...P1
1744 // (since P0 and P1 are always inside the rectangle)
1745
1746 for (i=0; i<4; i++) {
1747 if (t[i]>0) {
1748 if (tp>t[i])
1749 tp = t[i];
1750 }
1751 else {
1752 if (tn<t[i])
1753 tn = t[i];
1754 }
1755 }
1756 // printf("TN: %5.2f; TP: %5.2f\n",tn,tp);
1757
1758 if ((tn==MAXNEG) || (tp==MAXPOS))
1759 return FALSE; // It shoudn't happen!
1760
1761 *startx = ((gdouble) x0) + dx *tn;
1762 *starty = ((gdouble) y0) + dy * tn;
1763 *endx = ((gdouble) x0) + dx * tp;
1764 *endy = ((gdouble) y0) + dy * tp;
1765
1766 return TRUE;
1767 }
1768
hexagonal_row_sampling_with_type(gpointer hf,gint hf_size,gboolean wrap,gboolean even,gboolean shift_right,gint data_type,gint axis)1769 gpointer hexagonal_row_sampling_with_type (gpointer hf, gint hf_size,
1770 gboolean wrap, gboolean even, gboolean shift_right, gint data_type, gint axis) {
1771 // Creates a new hf map, translating each odd or even row by 0.5 pixel
1772 // This gives a rough approximation of a hexagonal tiling
1773 // This could improve some transformations, like the erosion
1774 // We use a very rough approximation, for now (average)
1775 // If (even), we modify the even rows, otherwise we modify the odd rows
1776 // If (shift_right), the row is translated to the right
1777 // (we put in the current pixel half of its value plus half
1778 // of the value of the preceding pixel)
1779 // If (!shift_right), the row is translated to the left
1780 hf_type *hf_out;
1781 gint x, y, i, j, shift;
1782 gboolean test;
1783 if (shift_right)
1784 shift = -1;
1785 else
1786 shift = 1;
1787 switch (data_type) {
1788 case GINT_ID:
1789 hf_out = (gpointer) x_malloc(sizeof(gint)*hf_size*hf_size, "gpointer (hf_out in hexagonal_row_sampling_with_type)");
1790 break;
1791 case HF_TYPE_ID:
1792 hf_out = (gpointer) x_malloc(sizeof(hf_type)*hf_size*hf_size, "gpointer (hf_out in hexagonal_row_sampling_with_type)");
1793 break;
1794 case GDOUBLE_ID:
1795 hf_out = (gpointer) x_malloc(sizeof(gdouble)*hf_size*hf_size, "gpointer (hf_out in hexagonal_row_sampling_with_type)");
1796 break;
1797 case UNSIGNED_CHAR_ID:
1798 hf_out = (gpointer) x_malloc(sizeof(unsigned char)*hf_size*hf_size, "gpointer (hf_out in hexagonal_row_sampling_with_type)");
1799 break;
1800 default:
1801 printf("Unexpected data type in hexagonal_row_sampling_with_type\n");
1802 exit(0);
1803 }
1804 for (y=0; y<hf_size; y++) {
1805 test = !(y & (gint) 1); // test = TRUE if value is odd
1806 for (x=0; x<hf_size; x++) {
1807 if (axis==H_AXIS)
1808 // i = VECTORIZE(x+shift, y, hf_size);
1809 i = VECTORIZE( WRAP2( x + shift, hf_size), y, hf_size);
1810 else
1811 // i = VECTORIZE(y, x+shift, hf_size);
1812 i = VECTORIZE( y, (WRAP2( x + shift, hf_size)), hf_size);
1813 if ( (test && even) || ((!test) && (!even)) ) { // XOR
1814 switch (data_type) {
1815 case GINT_ID:
1816 *(((gint *) hf_out)+i) = *(((gint *) hf)+i);
1817 break;
1818 case HF_TYPE_ID:
1819 *(((hf_type *) hf_out)+i) = *(((hf_type *) hf)+i);
1820 break;
1821 case GDOUBLE_ID:
1822 *(((gdouble *) hf_out)+i) = *(((gdouble *) hf)+i);
1823 break;
1824 case UNSIGNED_CHAR_ID:
1825 *(((unsigned char *) hf_out)+i) = *(((unsigned char *) hf)+i);
1826 break;
1827 default:
1828 printf("Unexpected data type in hexagonal_row_sampling_with_type\n");
1829 exit(0);
1830 }
1831 continue;
1832 }
1833 if (wrap) {
1834 if (axis==H_AXIS)
1835 j = VECTORIZE( WRAP2( x + shift, hf_size), y, hf_size);
1836 else
1837 j = VECTORIZE( y, (WRAP2( x + shift, hf_size)), hf_size);
1838 }
1839 else { // Repeat the pixel
1840 if (axis==H_AXIS)
1841 j = VECTORIZE(MAX( MIN(hf_size, x + shift),0), y, hf_size);
1842 else
1843 j = VECTORIZE(y, MAX( MIN(hf_size, x + shift),0), hf_size);
1844 }
1845
1846 switch (data_type) {
1847 case GINT_ID:
1848 *(((gint *) hf_out)+i) = (gint) ((((glong) *(((gint *) hf)+i) ) + (glong) *(((gint *) hf)+j) )>>1 );
1849 break;
1850 case HF_TYPE_ID:
1851 *(((hf_type *) hf_out) + i) = (hf_type) ((((glong) *(((hf_type *) hf)+i) ) + (glong) *(((hf_type *) hf)+j) )>>1 );
1852 break;
1853 case GDOUBLE_ID:
1854 *(((gdouble *) hf_out)+i) = ( *(((gdouble *) hf)+i) + *(((gdouble *) hf)+j) ) / 2.0 ;
1855 break;
1856 case UNSIGNED_CHAR_ID:
1857 *(((unsigned char *) hf_out)+i) = (unsigned char) ((((glong) *(((unsigned char *) hf)+i) ) + (glong) *(((unsigned char *) hf)+j) )>>1 );
1858 break;
1859 default:
1860 printf("Unexpected data type in hexagonal_row_sampling_with_type\n");
1861 exit(0);
1862 }
1863 // printf("SHIFTING (%d, %d) = i: %d; j: %d\n",x,y,i,j);
1864 }
1865 }
1866 return hf_out;
1867 }
1868
hexagonal_row_sampling(hf_type * hf,gint hf_size,gboolean wrap,gboolean even,gboolean shift_right)1869 hf_type *hexagonal_row_sampling (hf_type *hf, gint hf_size,
1870 gboolean wrap, gboolean even, gboolean shift_right) {
1871
1872 return (hf_type *) hexagonal_row_sampling_with_type (hf, hf_size, wrap, even, shift_right, HF_TYPE_ID, V_AXIS);
1873 }
1874
add_spread_3x3(hf_type * hf,gint max,gint x,gint y,gint value,gboolean wrap)1875 void add_spread_3x3 (hf_type *hf, gint max, gint x, gint y, gint value, gboolean wrap) {
1876 // Add value in hf at point(x,y), spreading the value given "matrix"
1877 // (3x3 convolution)
1878 static gdouble matrix[3][3]={ {0.5, 1.0, 0.5}, {1.0, 8.0, 1.0}, {0.5, 1.0, 0.5}};
1879 // static gdouble matrix[3][3]={{1.0, 2.0, 1.0}, {2.0, 4.0, 2.0}, {1.0, 2.0, 1.0}};
1880 static gdouble sum = 14.0; // Sum of the matrix
1881 gint i, m, n, xi, yi;
1882 gdouble qty;
1883
1884 for (m=0; m<3; m++)
1885 for (n=0; n<3; n++) {
1886 if (wrap)
1887 i = VECTORIZE(WRAP2(x-1+n,max),WRAP2(y-1+m,max),max);
1888 else {
1889 xi = x-1+n;
1890 yi = y-1+m;
1891 if ( (xi<0) || (yi<0) || (xi>=max) || (yi>=max) )
1892 continue;
1893 else
1894 i = VECTORIZE(xi,yi,max);
1895 }
1896 qty = (((gdouble)value) * matrix[m][n]) / sum;
1897 *(hf+i) = (hf_type) MIN(MAX_HF_VALUE,MAX(0,((gint) *(hf+i)) + (gint) qty));
1898 }
1899 }
1900
hf_clamp_buffer(hf_type * buffer,gint length,hf_type newmin,hf_type newmax)1901 void hf_clamp_buffer (hf_type *buffer, gint length, hf_type newmin, hf_type newmax) {
1902 // "Clamp" values between newmin and newmax
1903 // buffer is considered as a 1D array
1904 gfloat ratio;
1905 gint i;
1906 hf_type vmin, vmax;
1907 vmin = *buffer;
1908 vmax = vmin;
1909 // We first need to know the max and min value of the HF
1910 for (i=0; i<length; i++) {
1911 if (vmin>*(buffer+i) )
1912 vmin = *(buffer+i);
1913 else
1914 if (vmax<*(buffer+i))
1915 vmax = *(buffer+i);
1916 }
1917 if (vmin == vmax)
1918 return;
1919 ratio = ((gfloat) ABS(newmax-newmin)) / (gfloat) (vmax - vmin);
1920 for (i=0; i<length; i++) {
1921 *(buffer+i) = newmin + (hf_type) (ratio * (gfloat) (*(buffer+i) - vmin) ) ;
1922 }
1923 }
1924
hf_clamp(hf_struct_type * hf,hf_type newmin,hf_type newmax)1925 void hf_clamp (hf_struct_type *hf, hf_type newmin, hf_type newmax) {
1926 // "Clamp" values between newmin and newmax
1927 // Version using the 16 bits display buffer
1928 hf_clamp_buffer(hf->hf_buf, hf->max_x*hf->max_y, newmin, newmax);
1929 }
1930
double_clamp(hf_type * hf_buf,gint max_x,gint max_y,hf_type newmin,hf_type newmax,gdouble * dbuf)1931 void double_clamp (hf_type *hf_buf, gint max_x, gint max_y, hf_type newmin, hf_type newmax, gdouble *dbuf) {
1932 gdouble ratio, min, max=0.0, newmind;
1933 static gdouble maxd = (gdouble) MAX_HF_VALUE;
1934 gint i;
1935 newmind = (gdouble) newmin;
1936 min = maxd;
1937 for (i=0; i<max_x*max_y; i++) {
1938 if (*(dbuf+i)<min)
1939 min = *(dbuf+i);
1940 else
1941 if (*(dbuf+i)>max)
1942 max = *(dbuf+i);
1943 // if (!((i+(max_x>>1))%max_x))
1944 // printf("%p; @i= %d :: min=%f;max=%f\n",dbuf,i,min,max);
1945 }
1946 if (min == max)
1947 ratio = 0.0;
1948 else
1949 ratio = ((gdouble) ABS(newmax-newmin)) / (max - min);
1950 // printf("NEWMIND: %f; MIN: %f; MAX: %f; RATIO: %f\n",newmind, min,max,ratio);
1951 for (i=0; i< max_x*max_y; i++) {
1952 *(hf_buf+i) = (hf_type) ( (newmind + ratio * ( ( (gdouble) *(dbuf+i)) - min)));
1953 }
1954 }
1955
hf_double_clamp(hf_struct_type * hf,hf_type newmin,hf_type newmax,gdouble * dbuf)1956 void hf_double_clamp (hf_struct_type *hf, hf_type newmin, hf_type newmax, gdouble *dbuf) {
1957 // Same as "hf clamp", but with a double float buffer as input
1958 double_clamp (hf->hf_buf, hf->max_x, hf->max_y, newmin, newmax, dbuf);
1959 }
1960
hf_type_to_double_scaled(hf_type * hf_buf,gint max_x,gint max_y,gdouble * doublebuf)1961 void hf_type_to_double_scaled (hf_type *hf_buf, gint max_x, gint max_y, gdouble *doublebuf) {
1962 // Convert hf_buf to double precision data
1963 // Values are scaled from 0 - MAX_HF_VALUE to 0.0 - 1.0
1964 // doublebuf MUST be allocated
1965 static gdouble maxd = (gdouble) MAX_HF_VALUE;
1966 gint i;
1967 for (i=0; i<max_x*max_y; i++) {
1968 *(doublebuf+i) = ((gdouble) *(hf_buf+i)) / maxd;
1969 }
1970 }
hf_type_to_double(hf_type * hf_buf,gint max_x,gint max_y,gdouble * doublebuf)1971 void hf_type_to_double (hf_type *hf_buf, gint max_x, gint max_y, gdouble *doublebuf) {
1972 // Convert hf_buf to double precision data
1973 // doublebuf MUST be allocated
1974 gint i;
1975 for (i=0; i<max_x*max_y; i++) {
1976 *(doublebuf+i) = (gdouble) *(hf_buf+i);
1977 }
1978 }
1979
hf_subtract(hf_type * hf1,hf_type * hf2,hf_type * result,gint length,gint behaviour)1980 void hf_subtract(hf_type *hf1, hf_type *hf2, hf_type *result, gint length, gint behaviour) {
1981 // Simple subtraction of hf2 from hf1, in result
1982 // Already allocated buffers
1983 // behaviour is OVERFLOW_ZERO, OVERFLOW_REBOUND (absolute value) or
1984 // OVERFLOW_OFFSET (set the negative minimum value to 0)
1985 // OVERFLOW_SCALE (set the negative minimum value to 0, and scale down
1986 // if required)
1987 gint i;
1988 glong min, max, *array = NULL;
1989 gdouble ratio;
1990 if ((!hf1) || (!hf2) || (!result))
1991 return;
1992 if ( (behaviour==OVERFLOW_OFFSET) || (behaviour==OVERFLOW_SCALE) ) {
1993 // Calculate the minimum and maximum
1994 array = (glong*) x_malloc(sizeof(glong)*length, "glong (array in hf_subtract)");
1995 min = max = (*hf1) - *hf2;
1996 for (i=1; i<length; i++) {
1997 *(array+i) = ((glong) *(hf1+i)) - (glong) *(hf2+i);
1998 if (*(array+i)<min)
1999 min = *(array+i);
2000 else
2001 if (*(array+i)>max)
2002 max = *(array+i);
2003 }
2004 if ((max-min) < (glong) MAX_HF_VALUE)
2005 ratio = 1.0;
2006 else
2007 ratio = ((gdouble) MAX_HF_VALUE) / (gdouble) (max-min);
2008 // printf("MIN: %ld; MAX: %ld; Ratio: %7.4f\n",min,max,ratio);
2009 }
2010 for (i=0; i<length; i++) {
2011 switch (behaviour) {
2012 case OVERFLOW_ZERO:
2013 *(result+i) = (hf_type) MAX(0 , ((glong) *(hf1+i)) - (glong) *(hf2+i));
2014 break;
2015 case OVERFLOW_REBOUND:
2016 *(result+i) = (hf_type) ABS(((glong) *(hf1+i)) - (glong) *(hf2+i));
2017 break;
2018 case OVERFLOW_OFFSET:
2019 *(result+i) = (hf_type) (*(array+i) - min);
2020 break;
2021 case OVERFLOW_SCALE:
2022 *(result+i) = (hf_type) (((gdouble) (*(array+i) - min)) * ratio);
2023 }
2024 }
2025 if (array)
2026 x_free (array);
2027 }
2028
2029