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 #include <config.h>
24
25 #include <stdio.h>
26 #include <errno.h>
27 #include <error.h>
28 #include <stdlib.h>
29 #include <string.h>
30
31 #include <gnuastro/fits.h>
32 #include <gnuastro/tile.h>
33 #include <gnuastro/blank.h>
34 #include <gnuastro/threads.h>
35 #include <gnuastro/pointer.h>
36 #include <gnuastro/convolve.h>
37 #include <gnuastro/dimension.h>
38 #include <gnuastro/interpolate.h>
39 #include <gnuastro/permutation.h>
40
41 #include <gnuastro-internal/checkset.h>
42
43
44
45
46
47
48
49
50
51 /***********************************************************************/
52 /************** Single tile ******************/
53 /***********************************************************************/
54 /* Calculate the starting coordinates of a tile in the allocated block of
55 memory. */
56 void
gal_tile_start_coord(gal_data_t * tile,size_t * start_coord)57 gal_tile_start_coord(gal_data_t *tile, size_t *start_coord)
58 {
59 size_t ind, ndim=tile->ndim;
60 gal_data_t *block=gal_tile_block(tile);
61
62 /* If the input tile is actually the same as the block, then the start is
63 at 0 (in all dimensions). */
64 if(block==tile)
65 memset(start_coord, 0, ndim*gal_type_sizeof(GAL_TYPE_SIZE_T));
66 else
67 {
68 /* Calculate the coordinates of the first pixel of the tile. */
69 ind = gal_pointer_num_between(block->array, tile->array, block->type);
70 gal_dimension_index_to_coord(ind, ndim, block->dsize, start_coord);
71 }
72 }
73
74
75
76
77
78 /* Put the starting and ending (end point is not inclusive) coordinates of
79 a tile into the 'start_end' array. It is assumed that a space of
80 '2*tile->ndim' has been already allocated (static or dynamic) before
81 this function is called.
82
83 'rel_block' (or relative-to-block) is only relevant when the tile has an
84 intermediate tile between it and the allocated space (like a channel,
85 see 'gal_tile_full_two_layers'). If it doesn't ('tile->block' points the
86 allocated dataset), then the value to 'rel_block' is irrelevant.
87
88 However, when 'tile->block' is its self a larger block and 'rel_block'
89 is set to 0, then the starting and ending positions will be based on the
90 position within 'tile->block', not the allocated space. */
91 void
gal_tile_start_end_coord(gal_data_t * tile,size_t * start_end,int rel_block)92 gal_tile_start_end_coord(gal_data_t *tile, size_t *start_end, int rel_block)
93 {
94 size_t *s, *sf, *h;
95 gal_data_t *block=gal_tile_block(tile);
96 gal_data_t *host=rel_block ? block : tile->block;
97 size_t *hcoord, start_ind, ndim=tile->ndim, *end=start_end+ndim;
98
99 /* Get the starting index. Note that for the type we need the allocated
100 block dataset and can't rely on the tiles. */
101 start_ind=gal_pointer_num_between(block->array, tile->array, block->type);
102
103 /* Get the coordinates of the starting point relative to the allocated
104 block. */
105 gal_dimension_index_to_coord(start_ind, ndim, block->dsize, start_end);
106
107 /* When the host is different from the block, the tile's starting
108 position needs to be corrected. */
109 if(host!=block)
110 {
111 /* Get the host's starting coordinates. */
112 start_ind=gal_pointer_num_between(block->array, host->array,
113 block->type);
114
115 /* Temporarily put the host's coordinates in the place held for the
116 ending coordinates. */
117 hcoord=end;
118 gal_dimension_index_to_coord(start_ind, ndim, block->dsize, hcoord);
119 sf=(s=start_end)+ndim; h=hcoord; do *s++ -= *h++; while(s<sf);
120 }
121
122 /* Add the dimensions of the tile to the starting coordinate. Note that
123 the ending coordinates are stored immediately after the start.*/
124 gal_dimension_add_coords(start_end, tile->dsize, end, ndim);
125 }
126
127
128
129
130
131 /* Put the indexs of the first/start and last/end pixels (inclusive) in a
132 tile into the 'start_end' array (that has two elements). It will then
133 return the pointer to the start of the tile in the 'work' data
134 structure. */
135 void *
gal_tile_start_end_ind_inclusive(gal_data_t * tile,gal_data_t * work,size_t * start_end_inc)136 gal_tile_start_end_ind_inclusive(gal_data_t *tile, gal_data_t *work,
137 size_t *start_end_inc)
138 {
139 gal_data_t *block=gal_tile_block(tile);
140 size_t ndim=tile->ndim, *s, *e, *l, *sf;
141 size_t *start_coord = gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0,
142 __func__, "start_coord");
143 size_t *end_coord = gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0,
144 __func__, "end_coord");
145
146
147 /* The starting index can be found from the distance of the 'tile->array'
148 pointer and 'block->array' pointer. IMPORTANT: with the type of the
149 block array. */
150 start_end_inc[0]=gal_pointer_num_between(block->array, tile->array,
151 block->type);
152
153
154 /* To find the end index, we need to know the coordinates of the starting
155 point in the allocated block. */
156 gal_dimension_index_to_coord(start_end_inc[0], ndim, block->dsize,
157 start_coord);
158
159
160 /* 'end_coord' is one unit ahead of the last element in the tile in every
161 dimension. To have less potential for bugs, we will remove that extra
162 value, so we get the coordinates of the last pixel in the tile
163 (inclusive). We will finally, increment that value by one to get to
164 the pixel immediately outside of the tile.*/
165 e=end_coord;
166 l=tile->dsize;
167 sf=(s=start_coord)+ndim; do *e++ = *s + *l++ - 1; while(++s<sf);
168
169
170 /* Convert the (inclusive) ending point's coordinates into an index. */
171 start_end_inc[1]=gal_dimension_coord_to_index(ndim, block->dsize,
172 end_coord);
173
174
175 /* For a check:
176 printf("\ntile_dsize: %zu, %zu, %zu\n", tile->dsize[0], tile->dsize[1],
177 tile->dsize[2]);
178 printf("start_coord: %zu, %zu, %zu\n", start_coord[0], start_coord[1],
179 start_coord[2]);
180 printf("end_coord: %zu, %zu, %zu\n", end_coord[0], end_coord[1],
181 end_coord[2]);
182 printf("start_index: %zu\n", start_end_inc[0]);
183 printf("end_index: %zu\n", start_end_inc[1]);
184 exit(1);
185 */
186
187
188 /* Clean up and return the pointer in the work array that the tile starts
189 from. */
190 free(end_coord);
191 free(start_coord);
192 return gal_pointer_increment(work->array, start_end_inc[0], work->type);
193 }
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214 /***********************************************************************/
215 /************** Series of tiles ******************/
216 /***********************************************************************/
217 /* Construct a list of tile(s) given positional minimum(s) and maximum(s).
218 The output is an allocated an allocated array that can later be freed
219 with 'gal_data_array_free'. The minimum and maximums are assumed to be
220 inclusive.
221
222 The array keeping the minmium and maximum coordinates for each label
223 will have the following format:
224
225 | min0_d0 | min0_d1 | max0_d0 | max0_d1 | ...
226
227 ... | minN_d0 | minN_d1 | maxN_d0 | maxN_d1 | */
228 gal_data_t *
gal_tile_series_from_minmax(gal_data_t * block,size_t * minmax,size_t number)229 gal_tile_series_from_minmax(gal_data_t *block, size_t *minmax, size_t number)
230 {
231 size_t ndim=block->ndim;
232
233 size_t *min, *max;
234 size_t i, d, ind, size, width=2*ndim;
235 gal_data_t *tiles=gal_data_array_calloc(number);
236
237 /* Fill the tile information. */
238 for(i=0;i<number;++i)
239 {
240 /* To make things more readable. */
241 min = &minmax[ i * width ];
242 max = &minmax[ i * width + ndim ];
243
244 /* Tile types should be invalid (we shouldn't use tiles directly),
245 also se the other simple values. */
246 tiles[i].flag = 0;
247 tiles[i].block = block;
248 tiles[i].type = GAL_TYPE_INVALID;
249 tiles[i].next = i==number-1 ? NULL : &tiles[i+1];
250
251 /* Set the size related constants. */
252 size = 1;
253 tiles[i].ndim = ndim;
254 tiles[i].dsize = gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0,
255 __func__, "tiles[i].dsize");
256 for(d=0;d<ndim;++d) size *= tiles[i].dsize[d] = max[d] - min[d] + 1;
257 tiles[i].size = size;
258
259 /* Tile's array pointer. */
260 ind=gal_dimension_coord_to_index(ndim, block->dsize, min);
261 tiles[i].array = gal_pointer_increment(block->array, ind, block->type);
262 }
263
264 /* For a check (put all the objects in an extension of a test file).
265 {
266 gal_data_t *copy;
267 for(i=0;i<number;++i)
268 {
269 copy=gal_data_copy(&tiles[i]);
270 gal_fits_img_write(copy, "tiles.fits", NULL, NULL);
271 }
272 }
273 */
274
275 /* Return the final pointer. */
276 return tiles;
277 }
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298 /***********************************************************************/
299 /************** Allocated block of memory ******************/
300 /***********************************************************************/
301 /* When you are working on an array, it important to know the size of the
302 allocated space in each dimension. This simple function will just follow
303 the block pointer and return the 'dsize' element of lowest-level
304 structure. */
305 gal_data_t *
gal_tile_block(gal_data_t * tile)306 gal_tile_block(gal_data_t *tile)
307 {
308 while(tile->block!=NULL) tile=tile->block;
309 return tile;
310 }
311
312
313
314
315
316 /* Return the increment necessary to start at the next series of contiguous
317 memory (fastest dimension) associated with a tile.
318
319 1D and 2D cases are simple and need no extra explanation, but the case
320 for higher dimensions can be alittle more complicated, So we will go
321 over some examples. The notations below are:
322
323 'n' number of dimensions (same in tile and block).
324 't[]' size of the tile in each dimension.
325 'b[]' size of the allocated block in each dimension.
326
327 It is just important to see the output of this function as an increment
328 from the the last patch of contiguous memory associated with the
329 tile. So when the increment number is 't[n-1]' (the first 2D slice of
330 the tile has been parsed), simply incrementing by 'b[n-2] * b[n-1]' will
331 take us to the last row of
332
333 num_increment coord increment
334 ------------- ----- ---------
335 1 (...0,0,0) b[n-1]: fastest dimension of the block.
336 2 (...0,1,0) Similar to previous
337 . . .
338 . . .
339 t[n-2] (...1,0,0) (b[n-2] * b[n-1]) - ( (t[n-2]-1) * b[n-1] )
340 t[n-2] + 1 (...1,1,0) b[n-1]
341 . . .
342 . . .
343 2 * t[n-2] (...2,0,0) b[n-2] * b[n-1]
344 t[n-2]+1 (...2,1,0) b[n-1]
345 . . .
346 . . .
347 t[n-3] * t[n-2] (..1,0,0,0) b[n-3] * b[n-2] * b[n-1]
348
349 */
350 size_t
gal_tile_block_increment(gal_data_t * block,size_t * tsize,size_t num_increment,size_t * coord)351 gal_tile_block_increment(gal_data_t *block, size_t *tsize,
352 size_t num_increment, size_t *coord)
353 {
354 size_t n=block->ndim;
355 size_t *b=block->dsize, *t=tsize;
356 size_t increment=GAL_BLANK_SIZE_T;
357
358 if(n>3)
359 error(EXIT_FAILURE, 0, "%s: currently only implemented for at most 3 "
360 "dimensions", __func__);
361
362 switch(n)
363 {
364 /* A zero-dimensional dataset is not defined. */
365 case 0:
366 error(EXIT_FAILURE, 0, "%s: zero dimensional input is not acceptable",
367 __func__);
368
369 /* 1D: the increment is just the tile size. */
370 case 1:
371 increment=t[0];
372 if(coord) coord[0]+=increment;
373 break;
374
375 /* 2D: the increment is the block's number of fastest axis pixels. */
376 case 2:
377 increment=b[1];
378 if(coord) ++coord[0];
379 break;
380
381 /* 3D: The increment depends on which dimension we are reaching. */
382 case 3:
383 if(num_increment % t[1])
384 {
385 increment = b[2];
386 if(coord) ++coord[1];
387 }
388 else
389 {
390 increment=(b[1] * b[2]) - ( (t[1]-1) * b[2] );
391 if(coord) { ++coord[0]; coord[1] -= t[1]-1; coord[2]=0; }
392 }
393 break;
394 }
395
396 /* Return the final increment value. */
397 return increment;
398 }
399
400
401
402
403
404 /* Write a constant value for each tile into each pixel covered by the
405 input tiles in an array the size of the block and return it.
406
407 Arguments
408 ---------
409
410 'tilevalues': This must be an array that has the same number of
411 elements as that in 'tilesll' and in the same order that 'tilesll'
412 elements are parsed (from first to last). As a result the
413 dimensionality of this array is irrelevant. Note that unlike
414 'tiles', 'tilevalues' must be an array.
415
416 'tilesll': This will be parsed as a linked list (using the 'next'
417 element). Internally, it might be stored as an array, but this
418 function doesn't care! The position of the tile over its block will
419 be determined according to the 'block' element and the pointer of
420 its 'array' as fully described in 'gnuastro/data.h'. This function
421 will not pop/free the list, it will only parse it from start to
422 end.
423
424 'initialize': Initialize the allocated space with blank values before
425 writing in the constant values. This can be useful when the tiles
426 don't cover the full allocated block. */
427 gal_data_t *
gal_tile_block_write_const_value(gal_data_t * tilevalues,gal_data_t * tilesll,int withblank,int initialize)428 gal_tile_block_write_const_value(gal_data_t *tilevalues, gal_data_t *tilesll,
429 int withblank, int initialize)
430 {
431 void *in;
432 int type=tilevalues->type;
433 size_t tile_ind, nt=0, nv=tilevalues->size;
434 gal_data_t *tofill, *tile, *block=gal_tile_block(tilesll);
435
436 /* A small sanity check. */
437 for(tile=tilesll; tile!=NULL; tile=tile->next) ++nt;
438 if(nt!=nv)
439 error(EXIT_FAILURE, 0, "%s: the number of elements in 'tilevalues' (%zu) "
440 "and 'tilesll' (%zu) must be the same", __func__, nv, nt);
441
442 /* Allocate the output array. */
443 tofill=gal_data_alloc(NULL, type, block->ndim, block->dsize, block->wcs,
444 0, block->minmapsize, block->quietmmap,
445 tilevalues->name, tilevalues->unit,
446 tilevalues->comment);
447
448 /* If requested, initialize 'tofill', otherwise it is assumed that the
449 full area of the output is covered by the tiles. */
450 if(withblank || initialize) gal_blank_initialize(tofill);
451 else
452 {
453 /* Copy the flags. */
454 tofill->flag=tilevalues->flag;
455
456 /* If we have more than one dimension, then remove the possibly
457 sorted flags. */
458 if(block->ndim>1)
459 {
460 tofill->flag &= ~GAL_DATA_FLAG_SORTED_I;
461 tofill->flag &= ~GAL_DATA_FLAG_SORTED_D;
462 }
463 }
464
465 /* Go over the tiles and write the values in. Recall that 'tofill' has
466 the same type as 'tilevalues'. So we are using memcopy. */
467 tile_ind=0;
468 for(tile=tilesll; tile!=NULL; tile=tile->next)
469 {
470 /* Set the pointer to use as input. The 'if(o)' statement is set
471 because GCC 7.1.1 complained about the possiblity of the first
472 argument of 'memcpy' being NULL. Recall that 'o' is a pointer. */
473 in=gal_pointer_increment(tilevalues->array, tile_ind++, type);
474 GAL_TILE_PARSE_OPERATE( tile, tofill, 1, withblank, {
475 if(o) memcpy(o, in, gal_type_sizeof(type));
476 } );
477 }
478
479 return tofill;
480 }
481
482
483
484
485
486 /* Make a copy of the memory block and fill it with the index of each tile
487 in 'tilesll' (counting from 0). The non-filled areas will have blank
488 values. The output dataset will have a type of 'GAL_TYPE_INT32'. */
489 gal_data_t *
gal_tile_block_check_tiles(gal_data_t * tilesll)490 gal_tile_block_check_tiles(gal_data_t *tilesll)
491 {
492 int32_t *arr;
493 size_t i, dsize=gal_list_data_number(tilesll);
494 gal_data_t *ids, *out, *block=gal_tile_block(tilesll);
495
496 /* Allocate the array to keep the IDs of each tile. */
497 ids=gal_data_alloc(NULL, GAL_TYPE_INT32, 1, &dsize,
498 NULL, 0, block->minmapsize, block->quietmmap,
499 NULL, NULL, NULL);
500
501 /* Put the IDs into the array. */
502 arr=ids->array; for(i=0;i<dsize;++i) arr[i]=i;
503
504 /* Make the output. */
505 out=gal_tile_block_write_const_value(ids, tilesll, 0, 1);
506
507 /* Clean up and return. */
508 gal_data_free(ids);
509 return out;
510 }
511
512
513
514
515
516 /* Return the pointer corresponding to the tile in another data
517 structure (can have another type). */
518 void *
gal_tile_block_relative_to_other(gal_data_t * tile,gal_data_t * other)519 gal_tile_block_relative_to_other(gal_data_t *tile, gal_data_t *other)
520 {
521 gal_data_t *block=gal_tile_block(tile);
522 return gal_pointer_increment(other->array,
523 gal_pointer_num_between(block->array,
524 tile->array,
525 block->type),
526 other->type);
527 }
528
529
530
531
532
533 /* To use within 'gal_tile_full_blank_flag'. */
534 static void *
tile_block_blank_flag(void * in_prm)535 tile_block_blank_flag(void *in_prm)
536 {
537 struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
538 gal_data_t *tile_ll=(gal_data_t *)(tprm->params);
539
540 size_t i;
541 gal_data_t *tile;
542
543 /* Check all the tiles given to this thread. */
544 for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
545 {
546 tile=&tile_ll[ tprm->indexs[i] ];
547 gal_blank_present(tile, 1);
548 }
549
550 /* Wait for all the other threads to finish. */
551 if(tprm->b) pthread_barrier_wait(tprm->b);
552 return NULL;
553 }
554
555
556
557
558
559 /* Update the blank flag on the tiles within the list of input tiles. */
560 void
gal_tile_block_blank_flag(gal_data_t * tile_ll,size_t numthreads)561 gal_tile_block_blank_flag(gal_data_t *tile_ll, size_t numthreads)
562 {
563 /* Go over all the tiles and update their blank flag. */
564 gal_threads_spin_off(tile_block_blank_flag, tile_ll,
565 gal_list_data_number(tile_ll), numthreads,
566 tile_ll->minmapsize, tile_ll->quietmmap);
567 }
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588 /***********************************************************************/
589 /************** Tile full dataset ********************/
590 /***********************************************************************/
591 /* The user's specified tile size might not be an exact multiple of the
592 parent's size. This function is useful in such cases. It will give the
593 starting tile's size along each dimension.
594
595 The most simplistic way to manage the tiles is to put the regular tiles
596 at the start. The line below can be the length along any dimension, and
597 the tile size along that dimension.
598
599 | tile size | tile size | tile size | tile size | remainder
600 | | | | | |
601 ---------------------------------------------------------
602
603 The remainder of the scenario above will always be smaller than 'tile
604 size' (can be even 1-pixel wide). So, we will merge the first tile size
605 with the remainder. In this way, the size of the first tile will always
606 be between between one and two times the size of the regular tile:
607
608 | first tile | tile size | tile size | tile size |
609 | | | | |
610 ---------------------------------------------------------
611
612 When there is only a small remainder (for example one or two pixels),
613 then this layout is fine. But when the remainder is significant compared
614 to the regular tile size (like the example above), then it will make
615 more sense to cut the first tile into two halfs ('f-half' and 'l-half')
616 and put them at the start and end of the full length:
617
618
619 | f-half | tile size | tile size | tile size | l-half |
620 | | | | | |
621 ---------------------------------------------------------
622
623 So in any case, knowing the size of the first tile, will allow us to
624 parse all the tiles. We just have to make sure we don't go over the full
625 input's length. */
626 static void
gal_tile_full_regular_first(gal_data_t * parent,size_t * regular,float remainderfrac,size_t * first,size_t * last,size_t * tsize)627 gal_tile_full_regular_first(gal_data_t *parent, size_t *regular,
628 float remainderfrac, size_t *first, size_t *last,
629 size_t *tsize)
630 {
631 size_t i, remainder, *dsize=parent->dsize;;
632
633 /* For each dimension, set the size of the first tile. */
634 for(i=0;i<parent->ndim;++i)
635 {
636 /* It might happen that the tile size is bigger than the parent size
637 in a dimension, in that case the analysis in the comments above
638 are useless and only one tile should cover this dimension with the
639 size of the parent. */
640 if( regular[i] >= dsize[i] )
641 {
642 tsize[i]=1;
643 first[i]=last[i]=dsize[i];
644 }
645 else
646 {
647 /* Calculate the remainder in this dimension. */
648 remainder=dsize[i] % regular[i];
649
650 /* Depending on the remainder, set the first tile size and
651 number. */
652 if(remainder)
653 {
654 if( remainder > remainderfrac * regular[i] )
655 {
656 first[i] = ( remainder + regular[i] )/2;
657 tsize[i] = dsize[i]/regular[i] + 1 ;
658
659 /* If we only have one tile along the dimension, then
660 'first[i]==dsize[i]'. In this case, the first and last
661 tiles are the same and must have the same size. */
662 last[i] = ( first[i]==dsize[i]
663 ? first[i]
664 : ( dsize[i]
665 - ( first[i] + regular[i]*(tsize[i]-2) ) ) );
666 }
667 else
668 {
669 first[i] = remainder + regular[i];
670 tsize[i] = dsize[i]/regular[i];
671 last[i] = first[i]==dsize[i] ? first[i] : regular[i];
672 }
673 }
674 else
675 {
676 first[i] = last[i] = regular[i];
677 tsize[i] = dsize[i]/regular[i];
678 }
679 }
680 }
681
682 /* For a check:
683 printf("%s: first: %zu, %zu\n", __func__, first[0], first[1]);
684 printf("%s: last: %zu, %zu\n", __func__, last[0], last[1]);
685 */
686 }
687
688
689
690
691
692 /* Cover the full dataset with (mostly) identical tiles. The regular tile
693 size is determined from the 'size' array. If the input data's size is
694 not an exact multiple of 'size' for each dimension, then the tiles
695 touching the edges in that dimension will have a different size to fully
696 cover every element of the input. For a full description of tiling in
697 'gal_data_t', please see 'data.h'.
698
699 Inputs
700 ------
701
702 'input' is the gal_data_t which you want to tile (only used for its
703 sizes).
704
705 'regular' is the size of the regular tiles along each of the input's
706 dimensions. So it must have the same number of elements as the
707 dimensions of 'input'.
708
709 'remainderfrac' is the significant fraction of the remainder space if
710 the width of the input isn't an exact multiple of the tile size
711 along a dimension, see 'gal_tile_full_regular_first'.
712
713 'out' is the pointer to the array of data structures that is to keep
714 the tile parameters. If '*out==NULL', then the necessary space will
715 be allocated. If it is not NULL, then all the tile information will
716 be filled from the given element, see 'multiple' for more.
717
718 'multiple': When the '*out' array is to be allocated, allocate
719 'multiple' times the necessary space. This can be very useful when
720 you have several more identically sized 'inputs', and you want all
721 their tiles to be allocated (and thus indexed) together, even
722 though they have different 'block' datasets (that then link to one
723 allocated space). See the 'gal_tile_full_two_layers' below.
724
725 'firsttsize': The size of the first tile along every dimension. This
726 is only different from the regular tile size when 'regular' is not
727 an exact multiple of 'input''s length along every dimension. This
728 array is allocated internally by this function.
729
730 Output
731 ------
732
733 The returned output is an array of numbers (the same size as the input
734 data structure's dimensions) keeping the number of tiles along each
735 dimension.
736
737
738 Implementation
739 --------------
740
741 In the most general case, to set the starting pointers for each tile
742 we need the following sizes. If the input array has no parent/block,
743 then both these sizes are equal to it's own size:
744
745 1. block-size (or 'bsize'), which is the size of the allocated
746 block in each dimension.
747
748 2. parent-size (or 'psize') which is the size of the parent in each
749 dimension (we don't want to go out of the paren't range). */
750 size_t *
gal_tile_full(gal_data_t * input,size_t * regular,float remainderfrac,gal_data_t ** out,size_t multiple,size_t ** firsttsize)751 gal_tile_full(gal_data_t *input, size_t *regular,
752 float remainderfrac, gal_data_t **out, size_t multiple,
753 size_t **firsttsize)
754 {
755 size_t i, d, tind, numtiles, *start=NULL;
756 gal_data_t *tiles, *block=gal_tile_block(input);
757 size_t *last = gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim, 0,
758 __func__, "last");
759 size_t *first = gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim, 0,
760 __func__, "first");
761 size_t *coord = gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim, 0,
762 __func__, "coord");
763 size_t *tcoord = gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim, 0,
764 __func__, "tcoord");
765 size_t *tsize = gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim+1, 0,
766 __func__, "tsize");
767
768
769 /* Set the first tile size and total number of tiles along each
770 dimension, then allocate the array of tiles. */
771 gal_tile_full_regular_first(input, regular, remainderfrac,
772 first, last, tsize);
773 numtiles=gal_dimension_total_size(input->ndim, tsize);
774
775
776 /* Allocate the necessary space for all the tiles (if necessary). */
777 if(*out) tiles = *out;
778 else *out = tiles = gal_data_array_calloc(numtiles*multiple);
779
780
781 /* It is possible that the 'input' dataset is its-self a larger tile over
782 a region of the allocated block. In that case, we need to account for
783 the block's dimensions when calculating the position of this block. */
784 if(input->block)
785 {
786 start=gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim, 0, __func__,
787 "start");
788 gal_tile_start_coord(input, start);
789 }
790
791
792 /* Initialize each tile. */
793 for(i=0;i<numtiles;++i)
794 {
795 /* Specify the coordinates of the tile between the other tiles. Note
796 that we are dealing with tiles here, not pixels. */
797 gal_dimension_index_to_coord(i, input->ndim, tsize, tcoord);
798
799 /* The coordinates are currently in units of tiles, not
800 pixels. Convert them to the coordinates of the first pixel in each
801 tile. */
802 for(d=0;d<input->ndim;++d)
803 {
804 /* Convert the tile coordinates to pixel coordinates within
805 'input'. See the comments above 'gal_tile_full_regular_first':
806 The first tile in every dimension can be different from the
807 regular tile size. */
808 coord[d] = tcoord[d] ? first[d] + (tcoord[d]-1)*regular[d] : 0;
809
810 /* When the 'input' data structure (that is to be tiled here) was
811 itself a tile over a larger allocated array, a 'start' array
812 has been allocated to correct the coordinates so they refer to
813 a physical position on the allocated block of memory. */
814 if(start)
815 coord[d] += start[d];
816 }
817
818 /* Convert the coordinates (that are now in element/pixel units on
819 the allocated block of memory) into an index. */
820 tind=gal_dimension_coord_to_index(block->ndim, block->dsize, coord);
821
822 /* Now that we have the index of this tile's starting point compared
823 to the allocated block, put it in to the tile's 'array'
824 pointer. */
825 tiles[i].array=gal_pointer_increment(block->array, tind, block->type);
826
827 /* Set the sizes of the tile. */
828 tiles[i].size=1; /* Just an initializer, will be changed. */
829 tiles[i].ndim=input->ndim;
830 tiles[i].minmapsize=input->minmapsize;
831 tiles[i].dsize=gal_pointer_allocate(GAL_TYPE_SIZE_T,input->ndim, 0,
832 __func__, "tiles[i].dsize");
833 for(d=0;d<input->ndim;++d)
834 {
835 /* The size of the first and last tiles can be different from the
836 majority of the 'regular' tiles that have the same size. When
837 a tile is on the edge in one of the dimensions, then its
838 'tcoord[d]' will be either 0 or the last. */
839 if( first[d] != regular[d]
840 && ( tcoord[d]==0 || tcoord[d]==tsize[d]-1 ) )
841 {
842 if( tcoord[d] == 0 ) tiles[i].dsize[d] = first[d];
843 if( tcoord[d] == tsize[d]-1 ) tiles[i].dsize[d] = last[d];
844 }
845 else
846 tiles[i].dsize[d]=regular[d];
847
848 /* Set the size value. */
849 tiles[i].size *= tiles[i].dsize[d];
850 }
851
852 /* Set the block structure for this tile to the 'input', and set the
853 next pointer as the next tile. Note that only when we are dealing
854 with the last tile should the 'next' pointer be set to NULL.*/
855 tiles[i].flag = 0;
856 tiles[i].block = input;
857 tiles[i].next = i==numtiles-1 ? NULL : &tiles[i+1];
858
859 /* For a check:
860 printf("%zu:\n\tStart index: %zu\n\tsize: %zu x %zu\n", i, tind,
861 tiles[i].dsize[1], tiles[i].dsize[0]);
862 exit(0);
863 */
864 }
865
866
867 /* Clean up and return. */
868 free(last);
869 free(coord);
870 free(tcoord);
871 *firsttsize=first;
872 if(start) free(start);
873 tsize[input->ndim]=-1; /* 'tsize' had ndim+1 values, we will mark the */
874 return tsize; /* extra space with the largest possible value: */
875 } /* -1, see 'gal_tile_full_sanity_check'. */
876
877
878
879
880
881 /* Make sure that the input parameters (in 'tl', short for two-layer) fit
882 with the input dataset. The filename and HDU are only required for error
883 messages. Also, allocate and fill the 'channelsize' array. */
884 void
gal_tile_full_sanity_check(char * filename,char * hdu,gal_data_t * input,struct gal_tile_two_layer_params * tl)885 gal_tile_full_sanity_check(char *filename, char *hdu, gal_data_t *input,
886 struct gal_tile_two_layer_params *tl)
887 {
888 double d;
889 size_t i, ndim=input->ndim;
890
891 /* Check the tile's dimensions. */
892 for(i=0;tl->tilesize[i]!=-1;++i)
893 {
894 /* Not equal to zero. */
895 if(tl->tilesize[i]==0)
896 error(EXIT_FAILURE, 0, "'--tilesize' must be larger than zero, "
897 "the given value for dimension %zu was zero", ndim-i);
898
899 /* If the tile size is larger than the dataset size in this
900 dimension, then quietly change the tile size to the dataset size
901 along that dimension. */
902 if( tl->tilesize[i] > input->dsize[i] )
903 tl->tilesize[i] = input->dsize[i];
904 }
905
906
907 /* Make sure the number of tile sizes (tile dimensions) are the same as
908 the dataset's dimensions). */
909 if(i!=ndim)
910 error(EXIT_FAILURE, 0, "%s (hdu: %s): has %zu dimensions, but only %zu "
911 "value(s) given for the tile size ('--tilesize' option).",
912 filename, hdu, ndim, i);
913
914
915 /* Check the channel's dimensions. */
916 for(i=0; tl->numchannels[i]!=-1; ++i)
917 if(tl->numchannels[i]==0)
918 error(EXIT_FAILURE, 0, "the number of channels in all dimensions must "
919 "be larger than zero. The number for dimension %zu was zero",
920 i+1);
921 if(i!=ndim)
922 error(EXIT_FAILURE, 0, "%s (hdu: %s): has %zu dimensions, but only %zu "
923 "value(s) given for the number of channels", filename, hdu, ndim,
924 i);
925
926
927 /* Allocate space for the channel sizes. */
928 tl->channelsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
929 "tl->channelsize");
930
931
932 /* Check if the channels are exactly divisible by the input's size along
933 each dimension and set the correct size. */
934 for(i=0;i<ndim;++i)
935 {
936 /* Check if the number of channels is not more than the size of the
937 image. Note that the reported dimension must be in FITS format.*/
938 if( input->dsize[i] < tl->numchannels[i] )
939 error(EXIT_FAILURE, 0, "the number of channels in dimension %zu "
940 "(%zu) is more than the size of the '%s' (hdu: %s) in that "
941 "dimension", ndim-i, tl->numchannels[i], filename, hdu);
942
943 /* Also check the tile size. */
944 if( input->dsize[i] < tl->tilesize[i] )
945 error(EXIT_FAILURE, 0, "the tile size in dimension %zu (%zu) is "
946 "more than the size of the '%s' (hdu: %su) in that dimension",
947 ndim-i, tl->tilesize[i], filename, hdu);
948
949 /* First check. */
950 d=(double)input->dsize[i]/(double)(tl->numchannels[i]);
951 if(ceil(d)!=d)
952 error(EXIT_FAILURE, 0, "%zu (number of channels along dimension "
953 "%zu) is not exactly divisible by %zu (the length of '%s' "
954 "(hdu: %s) that dimension). The channels cover the input "
955 "dataset, hence, they must be identical", tl->numchannels[i],
956 ndim-i, input->dsize[i], filename, hdu);
957
958 /* Put the channel size into the output. */
959 tl->channelsize[i]=d;
960 }
961 }
962
963
964
965
966
967 /* A dataset can be tiled with two layers that are related:
968
969 Channels: A tesselation of larger tile sizes that all have the same
970 size ('channel_size' must be an exact multiple of 'input's size
971 along every dimension. In astronomy images, this can be seen as
972 CCD amplifiers, that cover large parts of the image. If
973 '*channels!=NULL' then it is assumed to be already present and
974 will not be allocated.
975
976 Tiles: A combined tesselation of each channel with smaller
977 tiles. These tiles can be used to calculate things like
978 gradients over each channel and thus over the whole image. */
979 void
gal_tile_full_two_layers(gal_data_t * input,struct gal_tile_two_layer_params * tl)980 gal_tile_full_two_layers(gal_data_t *input,
981 struct gal_tile_two_layer_params *tl)
982 {
983 gal_data_t *t;
984 size_t i, *junk, *junk2, ndim=tl->ndim=input->ndim;
985
986 /* Initialize. */
987 tl->channels=tl->tiles=NULL;
988
989 /* Initialize necessary values and do the channels tessellation. */
990 junk = gal_tile_full(input, tl->channelsize, tl->remainderfrac,
991 &tl->channels, 1, &junk2);
992 tl->totchannels = gal_dimension_total_size(ndim, tl->numchannels);
993 for(i=0;i<ndim;++i)
994 if(junk[i]!=tl->numchannels[i])
995 error(EXIT_FAILURE, 0, "%s: the input and output number of channels "
996 "don't match in dimension %zu: %zu and %zu respectively.",
997 __func__, ndim-i, tl->numchannels[i], junk[i]);
998 free(junk);
999 free(junk2);
1000
1001 /* Tile each channel. While tiling the first channel, we are also going
1002 to allocate the space for the other channels. Then pass those pointers
1003 when we want to fill in each tile of the other channels. */
1004 tl->numtilesinch = gal_tile_full(tl->channels, tl->tilesize,
1005 tl->remainderfrac, &tl->tiles,
1006 tl->totchannels, &tl->firsttsize);
1007 tl->tottilesinch = gal_dimension_total_size(ndim, tl->numtilesinch);
1008 for(i=1; i<tl->totchannels; ++i)
1009 {
1010 /* Set the first tile in this channel. Then use it it fill the 'next'
1011 pointer of the previous channel's tiles. Note that 'gal_tile_full'
1012 set this 'next' element to NULL. */
1013 t = tl->tiles + i * tl->tottilesinch;
1014 tl->tiles[ i * tl->tottilesinch - 1 ].next = t;
1015
1016 /* Fill in the information for all the tiles in this channel. Note
1017 that we already have the returned value, so it isn't important.*/
1018 junk=gal_tile_full(&tl->channels[i], tl->tilesize, tl->remainderfrac,
1019 &t, 1, &junk2);
1020 free(junk);
1021 free(junk2);
1022 }
1023
1024 /* Multiply the number of tiles along each dimension OF ONE CHANNEL by
1025 the number of channels in each dimension to get the dimensionality of
1026 the full tile structure. */
1027 tl->numtiles = gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
1028 "tl->numtiles");
1029 for(i=0;i<ndim;++i)
1030 tl->numtiles[i] = tl->numtilesinch[i] * tl->numchannels[i];
1031 tl->tottiles = gal_dimension_total_size(ndim, tl->numtiles);
1032 }
1033
1034
1035
1036
1037
1038 /* Usage
1039 -----
1040
1041 Make a permutation to allow the conversion of tile location in memory to
1042 its location in the full input dataset and put it in the input's
1043 'permutation' element. If a permutation has already been defined for the
1044 tessellation, this function will not do anythin. If permutation won't be
1045 necessary, then this function will just return (the permutation must
1046 have been initialized to NULL). */
1047 void
gal_tile_full_permutation(struct gal_tile_two_layer_params * tl)1048 gal_tile_full_permutation(struct gal_tile_two_layer_params *tl)
1049 {
1050 size_t *ch_coord, *tinch_coord;
1051 size_t i, p=0, t, ch, ind_in_all, ndim=tl->ndim;
1052
1053 /* If the permutation has already been defined for this tessellation,
1054 then there is no need to do it again here. */
1055 if(tl->permutation) return;
1056
1057 /* If there is only one channel or one dimension, return NULL. The
1058 permutation functions know that the input and output indexs are the
1059 same when the permutation is NULL. */
1060 if( ndim==1 || tl->totchannels==1) return;
1061
1062 /* Allocate the space for the permutation and coordinates. */
1063 ch_coord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
1064 "ch_coord");
1065 tinch_coord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
1066 "tinch_coord");
1067 tl->permutation=gal_pointer_allocate(GAL_TYPE_SIZE_T, tl->tottiles, 0,
1068 __func__, "tl->permutation");
1069
1070 /* Fill in the permutation, we use the fact that the tiles are filled
1071 from the first channel to the last. */
1072 for(ch=0;ch<tl->totchannels;++ch)
1073 {
1074 /* Get the coordinates of this channel's first tile. */
1075 gal_dimension_index_to_coord(ch, ndim, tl->numchannels, ch_coord);
1076 for(i=0;i<ndim;++i) ch_coord[i] *= tl->numtilesinch[i];
1077
1078 /* Go over all the tiles in this channel. */
1079 for(t=0;t<tl->tottilesinch;++t)
1080 {
1081 /* Convert its index to coordinates and add them to the channel's
1082 starting coordinates. */
1083 gal_dimension_index_to_coord(t, ndim, tl->numtilesinch,
1084 tinch_coord);
1085 for(i=0;i<ndim;++i) tinch_coord[i] += ch_coord[i];
1086
1087 /* Convert the coordinates into an index. */
1088 ind_in_all = gal_dimension_coord_to_index(ndim, tl->numtiles,
1089 tinch_coord);
1090 tl->permutation[ind_in_all] = p++;
1091 }
1092 }
1093
1094 /* Clean up and return. */
1095 free(tinch_coord);
1096 free(ch_coord);
1097 }
1098
1099
1100
1101
1102
1103 /* Write one value for each tile into a file.
1104
1105 IMPORTANT: it is assumed that the values are in the same order as the
1106 tiles.
1107
1108 tile[i] --> tilevalues[i] */
1109 void
gal_tile_full_values_write(gal_data_t * tilevalues,struct gal_tile_two_layer_params * tl,int withblank,char * filename,gal_fits_list_key_t * keys,char * program_string)1110 gal_tile_full_values_write(gal_data_t *tilevalues,
1111 struct gal_tile_two_layer_params *tl,
1112 int withblank, char *filename,
1113 gal_fits_list_key_t *keys, char *program_string)
1114 {
1115 gal_data_t *disp;
1116
1117 /* Make the dataset to be displayed. */
1118 if(tl->oneelempertile)
1119 {
1120 if(tl->ndim>1 && tl->totchannels>1)
1121 {
1122 /* A small sanity check. */
1123 if(tl->permutation==NULL)
1124 error(EXIT_FAILURE, 0, "%s: no permutation defined for the input "
1125 "tessellation", __func__);
1126
1127 /* Writing tile values to disk is not done for checking, not for
1128 efficiency. So to be safe (allow the caller to work on
1129 multiple threads), we will copy the tile values, then permute
1130 those. */
1131 disp = gal_data_copy(tilevalues);
1132 gal_permutation_apply(disp, tl->permutation);
1133 }
1134 else disp = tilevalues;
1135 }
1136 else
1137 disp=gal_tile_block_write_const_value(tilevalues, tl->tiles,
1138 withblank, 0);
1139
1140 /* Write the array as a file and then clean up (if necessary). */
1141 gal_fits_img_write(disp, filename, keys, program_string);
1142 if(disp!=tilevalues) gal_data_free(disp);
1143 }
1144
1145
1146
1147
1148
1149 /* Smooth the given values with a flat kernel of the given width. */
1150 gal_data_t *
gal_tile_full_values_smooth(gal_data_t * tilevalues,struct gal_tile_two_layer_params * tl,size_t width,size_t numthreads)1151 gal_tile_full_values_smooth(gal_data_t *tilevalues,
1152 struct gal_tile_two_layer_params *tl,
1153 size_t width, size_t numthreads)
1154 {
1155 size_t *kdsize, knum, i;
1156 gal_data_t *kernel, *smoothed;
1157 struct gal_tile_two_layer_params ttl={0};
1158 int permute=tl->ndim>1 && tl->totchannels>1;
1159
1160
1161 /* Check if the width is odd. */
1162 if(width%2==0)
1163 error(EXIT_FAILURE, 0, "%s: %zu not acceptable as width. It has to be "
1164 "an odd number", __func__, width);
1165
1166
1167 /* Prepare the kernel size along every dimension. */
1168 kdsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, tl->ndim, 0, __func__,
1169 "kdsize");
1170 for(i=0;i<tl->ndim;++i) kdsize[i]=width;
1171
1172
1173 /* Make the kernel. */
1174 kernel=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, tilevalues->ndim,
1175 kdsize, NULL, 0, -1, 1, NULL, NULL, NULL);
1176 knum=gal_dimension_total_size(tl->ndim, kernel->dsize);
1177 for(i=0;i<knum;++i) ((float *)(kernel->array))[i]=1/((double)knum);
1178
1179 /* Permute (if necessary). */
1180 if(permute)
1181 {
1182 gal_tile_full_permutation(tl);
1183 gal_permutation_apply(tilevalues, tl->permutation);
1184 }
1185
1186 /* Do the smoothing. */
1187 if(tl->workoverch)
1188 smoothed=gal_convolve_spatial(tilevalues, kernel, numthreads, 1, 1);
1189 else
1190 {
1191 /* Create the tile structure. */
1192 ttl.tilesize=tl->numtilesinch;
1193 ttl.numchannels=tl->numchannels;
1194 gal_tile_full_sanity_check("IMPOSSIBLE", "IMP_HDU", tilevalues, &ttl);
1195 gal_tile_full_two_layers(tilevalues, &ttl);
1196
1197 /* Do the convolution separately on each channel. */
1198 smoothed=gal_convolve_spatial(ttl.tiles, kernel, numthreads, 1, 0);
1199
1200 /* Clean up. */
1201 ttl.tilesize=ttl.numchannels=NULL;
1202 gal_tile_full_free_contents(&ttl);
1203 }
1204
1205 /* Reverse the permutation. */
1206 if(permute) gal_permutation_apply_inverse(smoothed, tl->permutation);
1207
1208 /* Clean up and return; */
1209 free(kdsize);
1210 gal_data_free(kernel);
1211 return smoothed;
1212 }
1213
1214
1215
1216
1217
1218 size_t
gal_tile_full_id_from_coord(struct gal_tile_two_layer_params * tl,size_t * coord)1219 gal_tile_full_id_from_coord(struct gal_tile_two_layer_params *tl,
1220 size_t *coord)
1221 {
1222 /* This function only works for 10 dimensions. */
1223 size_t i, tr, chid, tile[10];
1224
1225
1226 /* Host channel's ID. */
1227 for(i=0;i<tl->ndim;++i)
1228 tile[i] = tl->totchannels == 1 ? 0 : coord[i] / tl->channelsize[i];
1229 chid=gal_dimension_coord_to_index(tl->ndim, tl->numchannels, tile);
1230
1231
1232 /* Find the tile within the channel. */
1233 for(i=0;i<tl->ndim;++i)
1234 {
1235 tr=coord[i] % tl->channelsize[i];
1236 if( tl->firsttsize[i] != tl->tilesize[i] )
1237 tile[i] = ( tr <= tl->firsttsize[i]
1238 ? 0
1239 : 1 + (tr - tl->firsttsize[i]) / tl->tilesize[i] );
1240 else
1241 tile[i] = tr / tl->tilesize[i];
1242 }
1243
1244
1245 /* Return the tile ID. */
1246 return ( chid * tl->tottilesinch
1247 + gal_dimension_coord_to_index(tl->ndim, tl->numtilesinch, tile) );
1248 }
1249
1250
1251
1252
1253
1254
1255 /* Clean up the allocated spaces in the parameters. */
1256 void
gal_tile_full_free_contents(struct gal_tile_two_layer_params * tl)1257 gal_tile_full_free_contents(struct gal_tile_two_layer_params *tl)
1258 {
1259 /* Free the simply allocated spaces. */
1260 if(tl->tilesize) free(tl->tilesize);
1261 if(tl->numchannels) free(tl->numchannels);
1262 if(tl->channelsize) free(tl->channelsize);
1263 if(tl->numtiles) free(tl->numtiles);
1264 if(tl->numtilesinch) free(tl->numtilesinch);
1265 if(tl->tilecheckname) free(tl->tilecheckname);
1266 if(tl->permutation) free(tl->permutation);
1267 if(tl->firsttsize) free(tl->firsttsize);
1268
1269 /* Free the arrays of 'gal_data_t' for each tile and channel. */
1270 if(tl->tiles) gal_data_array_free(tl->tiles, tl->tottiles, 0);
1271 if(tl->channels) gal_data_array_free(tl->channels, tl->totchannels, 0);
1272 }
1273