1 /********************************************************************* 2 tile -- work with tesselations over a host dataset. 3 This is part of GNU Astronomy Utilities (Gnuastro) package. 4 5 Original author: 6 Mohammad Akhlaghi <mohammad@akhlaghi.org> 7 Contributing author(s): 8 Copyright (C) 2017-2021, Free Software Foundation, Inc. 9 10 Gnuastro is free software: you can redistribute it and/or modify it 11 under the terms of the GNU General Public License as published by the 12 Free Software Foundation, either version 3 of the License, or (at your 13 option) any later version. 14 15 Gnuastro is distributed in the hope that it will be useful, but 16 WITHOUT ANY WARRANTY; without even the implied warranty of 17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 18 General Public License for more details. 19 20 You should have received a copy of the GNU General Public License 21 along with Gnuastro. If not, see <http://www.gnu.org/licenses/>. 22 **********************************************************************/ 23 #ifndef __GAL_TILE_H__ 24 #define __GAL_TILE_H__ 25 26 /* Include other headers if necessary here. Note that other header files 27 must be included before the C++ preparations below */ 28 #include <gnuastro/data.h> 29 #include <gnuastro/fits.h> 30 #include <gnuastro/dimension.h> 31 32 /* C++ Preparations */ 33 #undef __BEGIN_C_DECLS 34 #undef __END_C_DECLS 35 #ifdef __cplusplus 36 # define __BEGIN_C_DECLS extern "C" { 37 # define __END_C_DECLS } 38 #else 39 # define __BEGIN_C_DECLS /* empty */ 40 # define __END_C_DECLS /* empty */ 41 #endif 42 /* End of C++ preparations */ 43 44 /* Actual header contants (the above were for the Pre-processor). */ 45 __BEGIN_C_DECLS /* From C++ preparations */ 46 47 48 49 50 51 /***********************************************************************/ 52 /************** Single tile ******************/ 53 /***********************************************************************/ 54 void 55 gal_tile_start_coord(gal_data_t *tile, size_t *start_coord); 56 57 void 58 gal_tile_start_end_coord(gal_data_t *tile, size_t *start_end, int rel_block); 59 60 void * 61 gal_tile_start_end_ind_inclusive(gal_data_t *tile, gal_data_t *work, 62 size_t *start_end_inc); 63 64 65 66 67 /***********************************************************************/ 68 /************** Series of tiles ******************/ 69 /***********************************************************************/ 70 gal_data_t * 71 gal_tile_series_from_minmax(gal_data_t *block, size_t *minmax, size_t number); 72 73 74 75 76 77 /***********************************************************************/ 78 /************** Allocated block **********************/ 79 /***********************************************************************/ 80 gal_data_t * 81 gal_tile_block(gal_data_t *tile); 82 83 size_t 84 gal_tile_block_increment(gal_data_t *block, size_t *tsize, 85 size_t num_increment, size_t *coord); 86 87 gal_data_t * 88 gal_tile_block_write_const_value(gal_data_t *tilevalues, gal_data_t *tilesll, 89 int withblank, int initialize); 90 91 gal_data_t * 92 gal_tile_block_check_tiles(gal_data_t *tiles); 93 94 void * 95 gal_tile_block_relative_to_other(gal_data_t *tile, gal_data_t *other); 96 97 void 98 gal_tile_block_blank_flag(gal_data_t *tile_ll, size_t numthreads); 99 100 101 102 103 104 /***********************************************************************/ 105 /************** Tile full dataset ********************/ 106 /***********************************************************************/ 107 108 struct gal_tile_two_layer_params 109 { 110 /* Inputs */ 111 size_t *tilesize; /* Tile size along each dim. (C order). */ 112 size_t *numchannels; /* Channel no. along each dim. (C order). */ 113 float remainderfrac; /* Frac. of remainers in each dim to cut. */ 114 uint8_t workoverch; /* Convolve over channel borders. */ 115 uint8_t checktiles; /* Tile IDs in an img, the size of input. */ 116 uint8_t oneelempertile; /* Only use one element for each tile. */ 117 118 /* Internal parameters. */ 119 size_t ndim; /* The number of dimensions. */ 120 size_t tottiles; /* Total number of tiles in all dim. */ 121 size_t tottilesinch; /* Number of tiles in one channel. */ 122 size_t totchannels; /* Total number of channels in all dim. */ 123 size_t *channelsize; /* Size of channels along each dimension. */ 124 size_t *numtiles; /* Tile no. in each dim. over-all. */ 125 size_t *numtilesinch; /* Tile no. in each dim. on one channel. */ 126 char *tilecheckname; /* Name of file to check tiles. */ 127 size_t *permutation; /* Tile pos. in memory --> pos. overall. */ 128 size_t *firsttsize; /* See 'gal_tile_full_regular_first'. */ 129 130 /* Actual tile and channel data structures. */ 131 gal_data_t *tiles; /* Tiles array (also linked with 'next'). */ 132 gal_data_t *channels; /* Channels array (linked with 'next'). */ 133 }; 134 135 136 size_t * 137 gal_tile_full(gal_data_t *input, size_t *regular, 138 float remainderfrac, gal_data_t **out, size_t multiple, 139 size_t **firsttsize); 140 141 void 142 gal_tile_full_sanity_check(char *filename, char *hdu, gal_data_t *input, 143 struct gal_tile_two_layer_params *tl); 144 145 void 146 gal_tile_full_two_layers(gal_data_t *input, 147 struct gal_tile_two_layer_params *tl); 148 149 void 150 gal_tile_full_permutation(struct gal_tile_two_layer_params *tl); 151 152 void 153 gal_tile_full_values_write(gal_data_t *tilevalues, 154 struct gal_tile_two_layer_params *tl, 155 int withblank, char *filename, 156 gal_fits_list_key_t *keys, char *program_string); 157 158 gal_data_t * 159 gal_tile_full_values_smooth(gal_data_t *tilevalues, 160 struct gal_tile_two_layer_params *tl, 161 size_t width, size_t numthreads); 162 163 size_t 164 gal_tile_full_id_from_coord(struct gal_tile_two_layer_params *tl, 165 size_t *coord); 166 167 void 168 gal_tile_full_blank_flag(gal_data_t *tile_ll, size_t numthreads); 169 170 void 171 gal_tile_full_free_contents(struct gal_tile_two_layer_params *tl); 172 173 174 175 176 177 /***********************************************************************/ 178 /************** Function-like macros ******************/ 179 /***********************************************************************/ 180 /* Useful when the input and other types are already known. We want this to 181 be self-sufficient (and be possible to call it independent of 182 'GAL_TILE_PARSE_OPERATE'), so some variables (basic definitions) that 183 are already defined in 'GAL_TILE_PARSE_OPERATE' re-defined here. */ 184 #define GAL_TILE_PO_OISET(IT, OT, IN, OTHER, PARSE_OTHER, CHECK_BLANK, OP) { \ 185 IT *i=IN->array; \ 186 gal_data_t *tpo_other=OTHER; /* 'OTHER' may be NULL. */ \ 187 gal_data_t *tpo_oblock = OTHER ? gal_tile_block(OTHER) : NULL; \ 188 \ 189 size_t tpo_s_e_i_junk[2]={0,0}; \ 190 IT b, *tpo_st=NULL, *tpo_f=i+IN->size; \ 191 size_t tpo_i_increment=0, tpo_num_i_inc=1; \ 192 size_t tpo_o_increment=0, tpo_num_o_inc=1; \ 193 int tpo_parse_other=(OTHER && PARSE_OTHER); \ 194 gal_data_t *tpo_iblock = gal_tile_block(IN); \ 195 OT *tpo_ost=NULL, *o = tpo_other ? tpo_other->array : NULL; \ 196 int tpo_hasblank = CHECK_BLANK ? gal_blank_present(IN, 0) : 0; \ 197 size_t tpo_s_e_i[2]={0,tpo_iblock->size-1}; /* -1: this is INCLUSIVE */ \ 198 \ 199 \ 200 /* A small sanity check: if 'OTHER' is given, and it is a block, */ \ 201 /* then it must have the same size as 'IN's block. On the other */ \ 202 /* hand, when 'OTHER' is a tile, its must have 'IN's size. */ \ 203 if( tpo_parse_other ) \ 204 { \ 205 if( OTHER==tpo_oblock ) /* 'OTHER' is a block. */ \ 206 { \ 207 if( gal_dimension_is_different(tpo_iblock, tpo_oblock) ) \ 208 { \ 209 /* 'error' function, is a GNU extension, see above. */ \ 210 fprintf(stderr, "GAL_TILE_PO_OISET: when " \ 211 "'PARSE_OTHER' is non-zero, the allocated " \ 212 "block size of 'IN' and 'OTHER' must be " \ 213 "equal, but they are not: %zu and %zu " \ 214 "elements respectively\n", tpo_iblock->size, \ 215 tpo_oblock->size); \ 216 exit(EXIT_FAILURE); \ 217 } \ 218 } \ 219 else \ 220 if( gal_dimension_is_different(IN, OTHER) ) \ 221 { \ 222 /* The 'error' function, is a GNU extension and this */ \ 223 /* is a header, not a library which the user has to */ \ 224 /* compile every time (on their own system). */ \ 225 fprintf(stderr, "GAL_TILE_PO_OISET: when " \ 226 "'PARSE_OTHER' is non-zero, the sizes of 'IN' " \ 227 "and 'OTHER' must be equal (in all " \ 228 "dimensions), but they are not: %zu and %zu " \ 229 "elements respectively\n", IN->size, \ 230 tpo_other->size); \ 231 exit(EXIT_FAILURE); \ 232 } \ 233 } \ 234 \ 235 \ 236 /* Write the blank value for the input type into 'b'. */ \ 237 gal_blank_write(&b, tpo_iblock->type); \ 238 \ 239 \ 240 /* If this is a tile, not a full block, then we need to set the */ \ 241 /* starting pointers ('tpo_st' and 'tpo_ost'). The latter needs */ \ 242 /* special attention: if it is a block, then we will use the */ \ 243 /* the same starting element as the input tile. If 'OTHER' is a */ \ 244 /* tile, then use its own starting position (recall that we have */ \ 245 /* already made sure that 'IN' and 'OTHER' have the same size. */ \ 246 if(IN!=tpo_iblock) \ 247 { \ 248 tpo_st = gal_tile_start_end_ind_inclusive(IN, tpo_iblock, \ 249 tpo_s_e_i); \ 250 if( tpo_parse_other ) \ 251 tpo_ost = ( OTHER==tpo_oblock \ 252 ? ( (OT *)(tpo_oblock->array) \ 253 + ( tpo_st - (IT *)(tpo_iblock->array) ) ) \ 254 : gal_tile_start_end_ind_inclusive(tpo_other, \ 255 tpo_oblock, \ 256 tpo_s_e_i_junk) ); \ 257 } \ 258 \ 259 \ 260 /* Go over contiguous patches of memory. */ \ 261 while( tpo_s_e_i[0] + tpo_i_increment <= tpo_s_e_i[1] ) \ 262 { \ 263 /* If we are on a tile, reset 'i' and 'o'. */ \ 264 if(IN!=tpo_iblock) \ 265 { \ 266 tpo_f = ( ( i = tpo_st + tpo_i_increment ) \ 267 + IN->dsize[IN->ndim-1] ); \ 268 if(tpo_parse_other) o = tpo_ost + tpo_o_increment; \ 269 } \ 270 \ 271 /* Do the operation depending the nature of the blank value. */ \ 272 /* Recall that for integer types, the blank value must be */ \ 273 /* checked with '=='. But for floats, the blank value can be */ \ 274 /* a NaN. Recall that a NAN will fail any comparison */ \ 275 /* including '=='. So when the blank value is not equal to */ \ 276 /* itself, then it is floating point and is a NAN. In that */ \ 277 /* case, the only way to check if a data element is blank or */ \ 278 /* not is to see if the value of each element is equal to */ \ 279 /* itself or not. */ \ 280 if(tpo_hasblank) \ 281 { \ 282 if(b==b) \ 283 do{if(*i!=b) {OP;} if(tpo_parse_other) ++o;} while(++i<tpo_f);\ 284 else \ 285 do{if(*i==*i) {OP;} if(tpo_parse_other) ++o;} while(++i<tpo_f);\ 286 } \ 287 else \ 288 do { {OP;} if(tpo_parse_other) ++o;} while(++i<tpo_f);\ 289 \ 290 /* Set the incrementation. On a fully allocated iblock (when */ \ 291 /* 'IN==tpo_iblock'), we have already gone through the whole */ \ 292 /* array, so we'll set the incrementation to the size of the */ \ 293 /* whole block. This will stop the 'while' loop above. On a */ \ 294 /* tile, we need to increment to the next contiguous patch */ \ 295 /* of memory to continue parsing this tile. */ \ 296 tpo_i_increment += ( IN==tpo_iblock \ 297 ? tpo_iblock->size \ 298 : gal_tile_block_increment(tpo_iblock, \ 299 IN->dsize, \ 300 tpo_num_i_inc++, \ 301 NULL) ); \ 302 \ 303 /* Similarly, increment the other array if necessary. Like */ \ 304 /* the above, when 'OTHER' is a full block, we'll just use */ \ 305 /* the same increment as 'IN'. Otherwise, when 'OTHER' is a */ \ 306 /* tile, calculate its increment based on its own block. */ \ 307 if(tpo_parse_other) \ 308 { \ 309 if(OTHER==tpo_oblock) tpo_o_increment=tpo_i_increment; \ 310 else \ 311 tpo_o_increment += gal_tile_block_increment(tpo_oblock, \ 312 tpo_other->dsize, \ 313 tpo_num_o_inc++, \ 314 NULL); \ 315 } \ 316 } \ 317 \ 318 /* This is done in case the caller doesn't need 'o' to avoid */ \ 319 /* compiler warnings. */ \ 320 o = o ? o+0 : NULL; \ 321 } 322 323 324 #define GAL_TILE_PO_OSET(OT, IN, OTHER, PARSE_OTHER, CHECK_BLANK, OP) { \ 325 switch(tpo_iblock->type) \ 326 { \ 327 case GAL_TYPE_UINT8: \ 328 GAL_TILE_PO_OISET(uint8_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 329 break; \ 330 case GAL_TYPE_INT8: \ 331 GAL_TILE_PO_OISET(int8_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 332 break; \ 333 case GAL_TYPE_UINT16: \ 334 GAL_TILE_PO_OISET(uint16_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 335 break; \ 336 case GAL_TYPE_INT16: \ 337 GAL_TILE_PO_OISET(int16_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 338 break; \ 339 case GAL_TYPE_UINT32: \ 340 GAL_TILE_PO_OISET(uint32_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 341 break; \ 342 case GAL_TYPE_INT32: \ 343 GAL_TILE_PO_OISET(int32_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 344 break; \ 345 case GAL_TYPE_UINT64: \ 346 GAL_TILE_PO_OISET(uint64_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 347 break; \ 348 case GAL_TYPE_INT64: \ 349 GAL_TILE_PO_OISET(int64_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 350 break; \ 351 case GAL_TYPE_FLOAT32: \ 352 GAL_TILE_PO_OISET(float, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 353 break; \ 354 case GAL_TYPE_FLOAT64: \ 355 GAL_TILE_PO_OISET(double, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 356 break; \ 357 default: \ 358 { /* 'error' function might not be available for the user. */ \ 359 fprintf(stderr, "GAL_TILE_PO_OSET: type code %d not recognized",\ 360 tpo_iblock->type); \ 361 exit(EXIT_FAILURE); \ 362 } \ 363 } \ 364 } 365 366 367 /* Parse over a region of memory (can be an n-dimensional tile or a fully 368 allocated block of memory) and do a certain operation. If 'OTHER' is not 369 NULL, this macro will also parse it at the same time . Note that OTHER 370 must either have only one element (for the whole input) or have exactly 371 the same number of elements as the input (one value for one 372 pixel/element of the input). See the documentation for more on this 373 macro and some examples. */ 374 #define GAL_TILE_PARSE_OPERATE(IN, OTHER, PARSE_OTHER, CHECK_BLANK, OP) { \ 375 gal_data_t *tpo_iblock = gal_tile_block(IN); \ 376 gal_data_t *tpo_oblock = OTHER ? gal_tile_block(OTHER) : NULL; \ 377 \ 378 /* First set the OTHER type. */ \ 379 if(OTHER) \ 380 switch(tpo_oblock->type) \ 381 { \ 382 case GAL_TYPE_UINT8: \ 383 GAL_TILE_PO_OSET(uint8_t, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 384 break; \ 385 case GAL_TYPE_INT8: \ 386 GAL_TILE_PO_OSET(int8_t, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 387 break; \ 388 case GAL_TYPE_UINT16: \ 389 GAL_TILE_PO_OSET(uint16_t,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 390 break; \ 391 case GAL_TYPE_INT16: \ 392 GAL_TILE_PO_OSET(int16_t, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 393 break; \ 394 case GAL_TYPE_UINT32: \ 395 GAL_TILE_PO_OSET(uint32_t,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 396 break; \ 397 case GAL_TYPE_INT32: \ 398 GAL_TILE_PO_OSET(int32_t, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 399 break; \ 400 case GAL_TYPE_UINT64: \ 401 GAL_TILE_PO_OSET(uint64_t,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 402 break; \ 403 case GAL_TYPE_INT64: \ 404 GAL_TILE_PO_OSET(int64_t, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 405 break; \ 406 case GAL_TYPE_FLOAT32: \ 407 GAL_TILE_PO_OSET(float, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 408 break; \ 409 case GAL_TYPE_FLOAT64: \ 410 GAL_TILE_PO_OSET(double, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 411 break; \ 412 default: \ 413 { \ 414 fprintf(stderr, "type code %d not recognized in " \ 415 "'GAL_TILE_PARSE_OPERATE'", tpo_oblock->type); \ 416 exit(EXIT_FAILURE); \ 417 } \ 418 } \ 419 else \ 420 /* When 'OTHER==NULL', its type is irrelevant, we'll just use */ \ 421 /*'int' as a place holder. */ \ 422 GAL_TILE_PO_OSET(int, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \ 423 } 424 425 426 427 __END_C_DECLS /* From C++ preparations */ 428 429 #endif 430