1 /* hf_calc.h - headers for hf transformation utilities 2 * 3 * Copyright (C) 2001-2010 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 #ifndef _HF_CALC 22 #define _HF_CALC 1 23 24 #include "globals.h" 25 #include "hf.h" 26 #include "merge.h" 27 #include <math.h> 28 29 // TEST... looking for a way to display the 16 bits lower heights 30 // in a 8 bits scale... See hf_16_to_8 in hf.c 31 #ifndef BLACK_THRESHOLD 32 #define BLACK_THRESHOLD 5.0 33 #endif 34 35 // Data types id 36 #ifndef GDOUBLE_ID 37 #define GDOUBLE_ID 0 38 #endif 39 40 #ifndef HF_TYPE_ID 41 #define HF_TYPE_ID 1 42 #endif 43 44 // Same as HF_TYPE 45 #ifndef UNSIGNED_SHORT_ID 46 #define UNSIGNED_SHORT_ID 1 47 #endif 48 49 #ifndef GINT_ID 50 #define GINT_ID 2 51 #endif 52 53 #ifndef UNSIGNED_CHAR_ID 54 #define UNSIGNED_CHAR_ID 3 55 #endif 56 57 #ifndef PIX8_RGB_ID 58 #define PIX8_RGB_ID 4 59 #endif 60 61 #ifndef PIX16_RGB_ID 62 #define PIX16_RGB_ID 5 63 #endif 64 65 // "E" 66 #ifndef CONST_E 67 #define CONST_E 2.718281828 68 #endif 69 70 #define hf_type_avrg(v1,v2) (hf_type)((((glong)(v1))+(glong)(v2))/2) 71 // Kind of bell curve, for using in "shape filters" and so on 72 73 // A cleaned version of the formula: 74 // 1 / pow( base, pow(2*index/size ,exp) / exp) 75 // The derived exponent to "base" is divided by "exp" 76 // When base == CONST_E and exp == 2, 77 // the formula gives a centered but not reduced gaussian (normal) bell 78 // One has to divide the result by pow(PI,0.5) to get a centered and reduced bell 79 // Actually we divide the independent variable (index) by half the radius 80 // so that the standard deviation is equal to half the radius, and we cover 95% of 81 // the values with the whole radius (given exp==2) 82 // The non reduced formula gives 1.0 at x==0. 83 84 // How to control the functions: 85 // exponent = "hardness" of the edges - greater = steepier 86 // size = relative width of the bell... start with something like half the actual radius (exponent = 2) 87 // it's probably not necessary to fiddle with "base" (check filters!) 88 89 #define BELL(base,exp,index,size) (gfloat)(1.0/ pow((gfloat)(base), pow(2.0*((gfloat)(index))/(gfloat)(size), (exp) ) / (exp) ) ) 90 91 #define BELLD(base,exp,index,size) (gdouble)(1.0/ pow((gdouble)(base), pow(2.0*((gdouble)(index))/(gdouble)(size), (exp) ) / (exp) ) ) 92 93 // Wrapping... ("tiling") 94 #define WRAP(x,max) (((x)+(max))%(max)) 95 96 // Case where x could be negative 97 #define WRAP2(x,max) ( ((x)<0) ? ((max)-WRAP(ABS(x),(max))) : WRAP((x),(max)) ) 98 99 // DWRAP is a float version of WRAP2 100 #define DOUBLEWRAP(x,max) ( ((x)<(max)) ? (x) : (x)-( ((gdouble) (max)) * floor((x)/(max)) ) ) 101 #define DWRAP(x,max) ( ((x)<0) ? ((max)-DOUBLEWRAP(ABS(x),(max))) : DOUBLEWRAP((x),(max)) ) 102 103 // An alternative to WRAP for processing overflowed indexes... 104 // X should be in base 0, and ranges from 0 to size-1 105 // Old version is valid only between -(size-1) and (2*size-1) 106 // #define REBOUND(x,size) ( ((x)>=size) ? (2*size-(x)-1) : ABS(x) ) 107 // New version 2004-04-11, valid for any value (indexes follow a "triangle" wave) 108 #define REBOUND(x,size) ( ((ABS(x)%(2*size))<size) ? (ABS(x)%size) : (2*size - (ABS(x)%(2*size)) - 1)) 109 110 // This macro returns the length of the intersection between two segments of an axis 111 // Allows calculating the area of the intersection between two rectangles 112 // No intersection = negative result - so we MAX the result to 0 113 #define INTERSECT(x1,size1,x2,size2) (MAX(0,MIN((x1)+(size1),(x2)+(size2))-MAX((x1),(x2)))) 114 115 // Redefinition of shift allowing negative values 116 // - interpreted as the reverse shift 117 #define RIGHT_SHIFT(value,shift) ( ((shift)>=0) ? (value)>>(shift) : (value)<<(-(shift)) ) 118 #define LEFT_SHIFT(value,shift) ( ((shift)>=0) ? (value)<<(shift) : (value)>>(-(shift)) ) 119 120 // As the name says... 121 #define FILL_MAP(map,size,value) for(;size-1;size--) *(map+size-1)=value 122 123 #define HF_BLACK 0x0000 124 #define HF_WHITE 0xFFFF 125 126 // Merge operations, for filters and paint strokes 127 128 // MULTIPLY: transform the multiplier as a ratio between 0.0 and 1.0 129 #ifndef MULTIPLY 130 #define MULTIPLY 1 131 #endif 132 133 // MULTIPLY2: transform the multiplier as a ratio between 1.0 and 2.0 134 #ifndef MULTIPLY2 135 #define MULTIPLY2 2 136 #endif 137 138 #ifndef ADD 139 #define ADD 3 140 #endif 141 #ifndef SUBTRACT 142 #define SUBTRACT 4 143 #endif 144 #ifndef MIN_MERGE 145 #define MIN_MERGE 5 146 #endif 147 #ifndef MAX_MERGE 148 #define MAX_MERGE 6 149 #endif 150 #ifndef XOR_MERGE 151 #define XOR_MERGE 7 152 #endif 153 154 #ifndef POWER_OP 155 #define POWER_OP 8 156 #endif 157 #ifndef LOG_OP 158 #define LOG_OP 9 159 #endif 160 #ifndef EXP_OP 161 #define EXP_OP 10 162 #endif 163 #ifndef BASE_OP 164 #define BASE_OP 11 165 #endif 166 #ifndef SINE_OP 167 #define SINE_OP 12 168 #endif 169 #ifndef SMOOTH_OP 170 #define SMOOTH_OP 13 171 #endif 172 #ifndef ALT_MERGE 173 #define ALT_MERGE 14 174 #endif 175 176 #ifndef MERGE_HARDNESS 177 #define MERGE_HARDNESS 0.15 178 #endif 179 180 // Cache of normalized gaussian vectors, generated by normalized_bell_new 181 #ifndef GAUSS_LIST_LENGTH 182 #define GAUSS_LIST_LENGTH 37 183 #endif 184 185 #ifndef TILING_AUTO 186 #define TILING_AUTO 0 187 #endif 188 #ifndef TILING_YES 189 #define TILING_YES 1 190 #endif 191 #ifndef TILING_NO 192 #define TILING_NO 2 193 #endif 194 195 // Parameters for overflow index processing 196 // in functions like interpolate and similar 197 198 // WRAP: max+1 -> 0, max+2 ->1...; -1->max, -2->max-1 199 // REBOUND: max+1 -> max-1, -1->1, -2->2... 200 // ZERO: outbound values set to 0 201 // IDLE: max+1, max+2... -> max 202 203 #ifndef OVERFLOW_WRAP 204 #define OVERFLOW_WRAP 1 205 #endif 206 #ifndef OVERFLOW_REBOUND 207 #define OVERFLOW_REBOUND 2 208 #endif 209 #ifndef OVERFLOW_ZERO 210 #define OVERFLOW_ZERO 3 211 #endif 212 #ifndef OVERFLOW_IDLE 213 #define OVERFLOW_IDLE 4 214 #endif 215 #ifndef OVERFLOW_OFFSET 216 #define OVERFLOW_OFFSET 5 217 #endif 218 #ifndef OVERFLOW_SCALE 219 #define OVERFLOW_SCALE 6 220 #endif 221 222 #ifndef NORTH 223 #define NORTH 0 224 #endif 225 226 #ifndef SOUTH 227 #define SOUTH 1 228 #endif 229 230 #ifndef EAST 231 #define EAST 2 232 #endif 233 234 #ifndef WEST 235 #define WEST 3 236 #endif 237 238 239 // The distance matrix (lookup table) 240 // We avoid recalculating distances from hf center for each use 241 // by storing them. Could become big for, say, 4096x4096 HF! 242 // We store only one quadrant and translate the coordinates with CENTER_XY! 243 244 // The most important rule: 245 // At any moment, hf_size for which the distance matrix is computed should be >= 246 // than the hf_size for which the current filter is computed, which should be >= 247 // than the hf_size 248 249 typedef struct { 250 gint hf_size; // HF size the current structure is able to process 251 gint size; // Vector size = square of the quadrant size 252 // Example: 256 for a 16 pixels wide square for a 32 pixels wide HF 253 gfloat *distances; // malloc(size*sizeof(gfloat)) 254 } dist_matrix_struct; 255 256 // This macro translates x,y coordinates relative to 0,0 257 // into absolute coordinates relative to the image center, for calculating distances 258 // The result is a X or Y index in a dist_matrix_struct 259 // Since we have an even number of pixels, the center is 2 pixels wide 260 #define CENTER_XY(xy,max) (((xy)>=((max)>>1))?(xy)-((max)>>1):((max)>>1)-(xy)-1) 261 262 // This macro "vectorizes" x,y coordinates, allowing the use of a * pointer instead of a [x,y] table 263 #define VECTORIZE(x,y,max) (((y)*(max))+(x)) 264 265 // Absolute value 266 #ifndef ABS 267 #define ABS(v) (((v)<0)?-(v):(v)) 268 #endif 269 270 // Getting the distance between 2 arbitrary points 271 #define DIST(dm,x1,y1,x2,y2) (*((dm)->distances+VECTORIZE(ABS((x2)-(x1)),ABS((y2)-(y1)),(dm)->hf_size>>1))) 272 // Straight version, without lookup table 273 #define POWDIST(x1,y1,x2,y2) ((x2)-(x1))*((x2)-(x1))+((y2)-(y1))*((y2)-(y1)) 274 #define DIST2(x1,y1,x2,y2) sqrt(POWDIST(x1,y1,x2,y2)) 275 276 // Scalar product of 2 vectors given by dots P1, P2, P3 277 // Vector 1 = P1-P2; Vector 2 = P3-P2 278 // Typically: for calculating the angle of a path 279 // If scalar product < 0, angle > 90 280 // If scalar product == 0, angle = 90 281 #define SCALAR_PRODUCT(x1,y1,x2,y2,x3,y3) (((x1)-(x2))*((x3)-(x2))+((y1)-(y2))*((y3)-(y2))) 282 // Same, with normalized vectors 283 #define NORM_SCALAR_PRODUCT(x1,y1,x2,y2,x3,y3) (((gdouble)SCALAR_PRODUCT(x1,y1,x2,y2,x3,y3))/((gdouble)(DIST2(x1,y1,x2,y2)*DIST2(x2,y2,x3,y3)))) 284 285 // Get the absolute slope value from degrees, for a given HF (1 pixel wide) 286 #ifndef dget_absolute_slope 287 #define dget_absolute_slope(deg,max) ((gdouble)(tan(PI*((gdouble)deg/180.0))*((gdouble)MAX_HF_VALUE)/(gdouble)max)) 288 #endif 289 290 #ifndef iget_absolute_slope 291 #define iget_absolute_slope(deg,max) ((gint)dget_absolute_slope(deg,max)) 292 #endif 293 294 #ifndef alloc_double_hf 295 #define alloc_double_hf(max) ((gdouble *)malloc(max*max*sizeof(gdouble))) 296 #endif 297 298 #ifndef calloc_double_hf 299 #define calloc_double_hf(max) ((gdouble *)calloc(max*max,sizeof(gdouble))) 300 #endif 301 302 #ifndef alloc_hf_type 303 #define alloc_hf_type(max) ((hf_type *)malloc(max*max*sizeof(hf_type))) 304 #endif 305 306 #ifndef calloc_hf_type 307 #define calloc_hf_type(max) ((hf_type *)calloc(max*max,sizeof(hf_type))) 308 #endif 309 310 #ifndef xalloc_double_hf 311 #define xalloc_double_hf(max,string) ((gdouble *)x_malloc(max*max*sizeof(gdouble),string)) 312 #endif 313 314 #ifndef xcalloc_double_hf 315 #define xcalloc_double_hf(max,string) ((gdouble *)x_calloc(max*max,sizeof(gdouble),string)) 316 #endif 317 318 #ifndef xalloc_hf_type 319 #define xalloc_hf_type(max,string) ((hf_type *)x_malloc(max*max*sizeof(hf_type),string)) 320 #endif 321 322 #ifndef xcalloc_hf_type 323 #define xcalloc_hf_type(max,string) ((hf_type *)x_calloc(max*max,sizeof(hf_type),string)) 324 #endif 325 326 #ifndef hf_min_merge 327 328 #define hf_min_merge(out,in1,in2,length) gint _i;for (_i=0;_i<length;_i++) *(out+_i)=MIN(*(in1+_i),*(in2+_i)) 329 330 #endif 331 332 #ifndef H_AXIS 333 #define H_AXIS 0 334 #endif 335 336 #ifndef V_AXIS 337 #define V_AXIS 1 338 #endif 339 340 // Standard vector with room for length 341 typedef struct { 342 gfloat x; 343 gfloat y; 344 gfloat z; 345 gfloat l; 346 } vector; 347 348 // Prototypes 349 350 void hf_merge( hf_struct_type *from_hf, 351 hf_struct_type *into_hf, 352 gint x, gint y, // "from_hf" is centered at "into_hf(x,y)" 353 gint merge_action, // ADD, SUBTRACT 354 gboolean if_tiles, 355 gfloat hf_displacement, 356 gboolean normalize); 357 358 // Same as hf_merge, decomposing the hf_struct_type* of the pen and the background 359 void hf_merge_buf (hf_type *hf_pencil, gint pen_max_x, gint pen_max_y, 360 hf_type *hf, gint max_x, gint max_y, 361 int x, int y, 362 gint pen_merge, 363 gboolean pen_tiles, 364 gfloat h_displacement, 365 gboolean normalize) ; 366 367 void generalized_merge( gpointer map, 368 gint map_data_type, 369 gint size_x, gint size_y, 370 gpointer hf, 371 gint hf_data_type, 372 gint max_x, gint max_y, 373 gint x, gint y, 374 gint merge_type, 375 gboolean wrap, 376 gboolean square_symmetry); 377 378 void interpolated_merge (gpointer map, 379 gint map_data_type, 380 gint size_x, gint size_y, 381 gpointer hf, 382 gint hf_data_type, 383 gint max_x, gint max_y, 384 gdouble x, 385 gdouble y, 386 gint merge_type, 387 gboolean wrap, 388 gboolean square_symmetry); 389 390 void hf_slide(hf_struct_type *hf, gint slideh, gint slidev) ; 391 392 void interpolate (gdouble dbx, gdouble dby, gpointer grid, 393 gint x_size, gint y_size, gpointer return_value_ptr, gint data_type, 394 gint overflow); 395 396 void translate_real_forward_mapping (gpointer source_grid, 397 gpointer output_grid, 398 gint data_type, 399 gint x_size, gint y_size, 400 gint x,gint y, 401 gdouble ox, gdouble oy) ; 402 403 void translate_forward_mapping (gpointer source_grid, 404 gpointer output_grid, 405 gint data_type, 406 gint x_size, gint y_size, 407 gint x,gint y, 408 gint ix, gint iy) ; 409 410 // Old version - use "interpolate" instead 411 hf_type interpolate2 (gdouble dbx, gdouble dby, hf_type *hf, gint x_size, gint y_size); 412 413 hf_type interpolate_get_xy (gdouble x, gdouble y, hf_type *hf, 414 gint x_size, gint y_size, 415 gdouble (*get_x) (gdouble, gpointer), gdouble (*get_y) (gdouble, gpointer), 416 gpointer arg_get_x, gpointer arg_get_y); 417 418 void vector_interpolate (gdouble dbx, gpointer vector, gint x_size, 419 gpointer return_value_ptr, gint data_type, gint overflow) ; 420 421 void rotate (gdouble dx, gdouble dy, gpointer map_in, gpointer map_out, 422 gint hsize, gint vsize, gint data_type, gint overflow); 423 424 void rotate_sin_cos (gdouble sin, gdouble cos, gpointer map_in, gpointer map_out, gint hsize, gint vsize, gint data_type, gint overflow); 425 426 void hf_rotate (hf_type *hf_in, hf_type *hf_out, gint hf_size, gint angle, 427 gboolean wrap, gint overflow); 428 void hf_fast_rotate (hf_type *hf_in, hf_type *hf_out, gint hf_size, gint angle); 429 430 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); 431 432 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); 433 434 hf_struct_type * hf_scale (hf_struct_type *hf_in, gint log_ratio); 435 436 void hf_horizontal_mirror(hf_struct_type *hf); 437 void hf_vertical_mirror(hf_struct_type *hf); 438 void hf_circular_projection(hf_type *hf_in, hf_type *hf_out, gint hf_size); 439 void hf_revert(hf_struct_type *hf); 440 441 void norm(vector *v); 442 443 // void hf_save(hf_struct_type *, gchar *); 444 445 dist_matrix_struct *dist_matrix_new(gint hf_size); 446 void dist_matrix_init(dist_matrix_struct *, gint hf_size); 447 gdouble *normalized_bell_new(gint radius); 448 void convolve_normalized_vector (gpointer in,gpointer out, gint max_x, gint max_y, gboolean wrap, 449 gint radius, gdouble *vector, gint data_type) ; 450 void map_convolve (hf_type *map, gint map_max_x, gint map_max_y, 451 hf_type *background, gint max_x, gint max_y, gint cx, gint cy, 452 gboolean wrap, gint level, gdouble **gauss_list, gboolean square_symmetry) ; 453 454 gboolean intersect_rect (gint startx_rect, gint starty_rect, gint endx_rect, gint endy_rect, 455 gint x0, gint y0, gint x1, gint y1, 456 gdouble *startx, gdouble *starty, gdouble *endx, gdouble *endy); 457 458 void fill_area (hf_type *buffer_in, hf_type *buffer_out, gint xmax, gint ymax, 459 hf_type filling_value, hf_type range, gint x, gint y); 460 461 gpointer hexagonal_row_sampling_with_type (gpointer hf, gint hf_size, 462 gboolean wrap, gboolean even, gboolean shift_right, gint data_type, gint axis) ; 463 464 hf_type *hexagonal_row_sampling (hf_type *hf, gint hf_size, 465 gboolean wrap, gboolean even, gboolean shift_right) ; 466 467 void add_spread_3x3 (hf_type *hf, gint max, gint x, gint y, gint value, gboolean wrap); 468 469 void hf_clamp (hf_struct_type *hf, hf_type min, hf_type max); 470 void hf_clamp_buffer (hf_type *buffer, gint length, 471 hf_type newmin, hf_type newmax); 472 void hf_double_clamp (hf_struct_type *hf, hf_type min, hf_type max, gdouble *buffer); 473 void double_clamp (hf_type *buf, gint max_x, gint max_y, hf_type newmin, hf_type newmax, gdouble *dbuf); 474 void hf_type_to_double (hf_type *hf_buf, gint max_x, gint max_y, gdouble *doublebuf); 475 void hf_type_to_double_scaled (hf_type *hf_buf, gint max_x, gint max_y, gdouble *doublebuf); 476 477 void hf_subtract(hf_type *hf1, hf_type *hf2, hf_type *result, gint length, gint behaviour); 478 479 #endif // _HF_CALC 480