[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 31b3b76 099/125: Spatial convolution now using
From: |
Mohammad Akhlaghi |
Subject: |
[gnuastro-commits] master 31b3b76 099/125: Spatial convolution now using new tessellation |
Date: |
Sun, 23 Apr 2017 22:36:47 -0400 (EDT) |
branch: master
commit 31b3b7647d8629f062d19f6b68939d9ce9025a7c
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Spatial convolution now using new tessellation
Spatial convolution (on multiple threads) is one of the most basic and
common applications of the tessellation management that was designed in the
previous commits. With this commit convolution on central tiles of the
image has been applied. Only the edge tiles remain. The new system has
significantly simplified the code (and probably the efficiency) of spatial
convolution in Gnuastro.
Also, reading of the input image has now been changed to float32
type. Previously it was float64, while the kernel was initially read as
float64. For frequency domain convolution this makes no difference, because
the input is put into a complex float64 type array before the processing is
done. For spatial domain convolution it also makes no difference because
the sum is stored as float64.
Finally, as before with the `mesh' library, a new
`gal_tile_function_on_threads' function has been defined to make it easy to
spin off threads working on separate tiles.
---
bin/convolve/convolve.c | 17 ++-
bin/convolve/ui.c | 2 +-
lib/convolve.c | 258 +++++++++++++++++++++++++++++++++++++++++-
lib/gnuastro/convolve.h | 5 +-
lib/gnuastro/tile.h | 41 +++++--
lib/tile.c | 291 ++++++++++++++++++++++++++++++++++++++----------
6 files changed, 535 insertions(+), 79 deletions(-)
diff --git a/bin/convolve/convolve.c b/bin/convolve/convolve.c
index 3c7d1c0..ff523d1 100644
--- a/bin/convolve/convolve.c
+++ b/bin/convolve/convolve.c
@@ -190,10 +190,10 @@ void
makepaddedcomplex(struct convolveparams *p)
{
size_t i, ps0, ps1;
- float *f, *ff, *kernel=p->kernel->array;
+ double *o, *op, *pimg, *pker;
size_t is0=p->input->dsize[0], is1=p->input->dsize[1];
size_t ks0=p->kernel->dsize[0], ks1=p->kernel->dsize[1];
- double *o, *op, *d, *df, *pimg, *pker, *input=p->input->array;
+ float *f, *ff, *input=p->input->array, *kernel=p->kernel->array;
/* Find the sizes of the padded image, note that since the kernel
@@ -220,8 +220,8 @@ makepaddedcomplex(struct convolveparams *p)
op=(o=pimg+i*2*ps1)+2*ps1; /* pimg is complex. */
if(i<is0)
{
- df=(d=input+i*is1)+is1;
- do {*o++=*d; *o++=0.0f;} while(++d<df);
+ ff=(f=input+i*is1)+is1;
+ do {*o++=*f; *o++=0.0f;} while(++f<ff);
}
do *o++=0.0f; while(o<op);
}
@@ -769,6 +769,7 @@ frequencyconvolve(struct convolveparams *p)
void
convolve(struct convolveparams *p)
{
+ gal_data_t *out;
gal_data_t *tiles, *channels=NULL; /* `channels' has to be initialized. */
/* Do the convolution. */
@@ -782,12 +783,18 @@ convolve(struct convolveparams *p)
if(p->tilesname)
gal_tile_block_check_tiles(tiles, p->tilesname, PROGRAM_NAME);
+
/* Do the spatial convolution. */
- gal_convolve_spatial(tiles, p->kernel);
+ out=gal_convolve_spatial(tiles, p->kernel, p->cp.numthreads,
+ p->convoverch);
/* Clean up */
+ gal_data_free(p->input);
gal_data_array_free(tiles, gal_data_num_in_ll(tiles), 0);
gal_data_array_free(channels, gal_data_num_in_ll(channels), 0);
+
+ /* Put the convolved dataset into the `input' to save as output. */
+ p->input=out;
}
else
frequencyconvolve(p);
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index 8ccf534..513c3ee 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -339,7 +339,7 @@ ui_preparations(struct convolveparams *p)
/* Read the input image as a float64 array and its WCS info. */
p->input=gal_fits_img_read_to_type(p->filename, p->cp.hdu,
- GAL_DATA_TYPE_FLOAT64, p->cp.minmapsize);
+ GAL_DATA_TYPE_FLOAT32, p->cp.minmapsize);
gal_wcs_read(p->filename, p->cp.hdu, 0, 0, &p->input->nwcs, &p->input->wcs);
diff --git a/lib/convolve.c b/lib/convolve.c
index 802c50d..47cc530 100644
--- a/lib/convolve.c
+++ b/lib/convolve.c
@@ -25,16 +25,238 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdio.h>
#include <errno.h>
#include <error.h>
+#include <string.h>
#include <stdlib.h>
#include <gnuastro/tile.h>
+#include <gnuastro/threads.h>
+#include <gnuastro/multidim.h>
+#include <gnuastro/convolve.h>
+
+
+
+
+
+
+
+
+/*********************************************************************/
+/******************** Utilities ********************/
+/*********************************************************************/
+/* See if the tile is on the edge of the hosted region or not. It doesn't
+ matter if the host is the allocated block of memory or a region in it (a
+ channel). */
+static int
+convolve_tile_is_on_edge(size_t *h, size_t *start_end_coord, size_t *k,
+ size_t ndim)
+{
+ size_t *s=start_end_coord, *e=start_end_coord+ndim, *kf=k+ndim;
+
+ /* If the starting point of the tile is smaller than the half-kernel
+ length along that dimension, then the tile is on the edge.
+
+ If the end of the tile in this dimension added with the half-kernel
+ length along that dimension is larger than the host's size, then the
+ tile is on the edge. */
+ do if( (*s++ < *k/2) || (*e++ + *k/2 > *h++) ) return 1; while(++k<kf);
+
+ /* If control reaches here, then this is a central tile. */
+ return 0;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*********************************************************************/
+/******************** Spatial convolution ********************/
+/*********************************************************************/
+struct spatial_params
+{
+ gal_data_t *out;
+ gal_data_t *tiles;
+ gal_data_t *kernel;
+ int convoverch;
+};
+
+
+
+
+
+/* Convolve over one tile that is not touching the edge. */
+static void
+convolve_spatial_tile(struct spatial_params *cprm, size_t tile_ind,
+ size_t *parse_coords, gal_data_t *overlap)
+{
+ gal_data_t *tile=&cprm->tiles[tile_ind];
+
+ int on_edge;
+ double sum, ksum;
+ size_t i, ndim=tile->ndim, csize=tile->dsize[ndim-1];
+ gal_data_t *block=gal_tile_block(tile), *kernel=cprm->kernel;
+
+ size_t i_inc, i_ninc, i_st_en[2], o_inc, o_ninc, o_st_en[2];
+ size_t *p, *pf, *k, *o, *e, o_st_inc, k_start, start_fastdim;
+
+ /* These variables depend on the type of the input. */
+ float *kv, *iv, *ivf, *i_start, *o_start;
+ float *in_v, *in=block->array, *out=cprm->out->array;
+
+ /* The `parse_coords' array was allocated once before this function for
+ all the tiles that are given to a thread. It has the space for all the
+ necessary coordinates. At first, `pix' will be the starting element of
+ the tile, but then it will be incremented as we carpet the tile. So we
+ aren't calling it `start'. */
+ size_t *pix = parse_coords;
+ size_t *end = parse_coords + ndim;
+ size_t *overlap_start = parse_coords + ndim * 2;
+ size_t *o_c = parse_coords + ndim * 3;
+
+
+ /* Set the starting and ending coordinates of this tile (recall that the
+ start and end are the first two allocated spaces in
+ parse_coords). When `convoverch' is set, we want to convolve over the
+ whole allocated block, not just one channel. So in effect, it is the
+ same as `rel_block' in `gal_tile_start_end_coord'. */
+ gal_tile_start_end_coord(tile, parse_coords, cprm->convoverch);
+ start_fastdim = pix[ndim-1];
+
+
+ /* See if this tile is on the edge or not. */
+ on_edge=convolve_tile_is_on_edge( ( cprm->convoverch
+ ? block->dsize
+ : tile->block->dsize ),
+ parse_coords, kernel->dsize, ndim);
+
+
+ /* Go over the tile, first find its limits, then loop over the fastest
+ dimension. */
+ i_inc=0; i_ninc=1;
+ i_start=gal_tile_start_end_ind_inclusive(tile, block, i_st_en);
+ while( i_st_en[0] + i_inc <= i_st_en[1] )
+ {
+ /* Initialize the value along the fastest dimension (it is not
+ incremented during `gal_tile_block_increment'). */
+ pix[ndim-1]=start_fastdim;
+
+ /* Go over each pixel to convolve. */
+ for(i=0;i<csize;++i)
+ {
+ /* Pointer to the pixel under consideration. */
+ in_v = i_start + i_inc + i;
+
+ /* If the input on this pixel is a NaN, then just set the output
+ to NaN too and go onto the next pixel. `in_v' is the pointer
+ on this pixel. */
+ if( isnan(*in_v) )
+ out[ in_v - in ]=NAN;
+ else
+ {
+ /* Coordinate to start convolution for this pixel. */
+ pf=(p=pix)+ndim; e=end;
+ o=overlap_start; k=kernel->dsize;
+ do
+ {
+ if(on_edge)
+ {
+ return;
+ }
+ else { *o++ = *p - *k++/2; k_start=0; }
+ }
+ while(++p<pf);
+ o_st_inc=gal_multidim_coord_to_index(ndim, block->dsize,
+ overlap_start);
+
+ /* Set the overlap parameters, then set the region. */
+ overlap->dsize=kernel->dsize;
+ overlap->array=gal_data_ptr_increment(block->array, o_st_inc,
+ block->type);
+ o_start=gal_tile_start_end_ind_inclusive(overlap, block,
+ o_st_en);
+
+ /* Go over the kernel-overlap region. */
+ kv = ( k_start
+ ? gal_data_ptr_increment(kernel->array, k_start,
+ kernel->type)
+ : kernel->array );
+ ksum=sum=0.0f; o_inc=0; o_ninc=1;
+ while( o_st_en[0] + o_inc <= o_st_en[1] )
+ {
+ /* Go over the contiguous region. */
+ ivf = ( iv = o_start + o_inc ) + overlap->dsize[ndim-1];
+ do
+ {
+ if( !isnan(*iv) )
+ { ksum += *kv; sum += *iv * *kv; }
+ ++kv;
+ }
+ while(++iv<ivf);
+
+ /* Update the incrementation to the next contiguous
+ region of memory over this tile. */
+ o_inc += gal_tile_block_increment(block, overlap->dsize,
+ o_ninc++, o_c);
+ }
+
+ /* Set the output value. */
+ out[ in_v - in ] = ksum==0.0f ? NAN : sum/ksum;
+ }
+
+ /* Increment the last coordinate. */
+ pix[ndim-1]++;
+ }
+
+ /* Increase the increment from the start of the tile for the next
+ contiguous patch. */
+ i_inc += gal_tile_block_increment(block, tile->dsize, i_ninc++, pix);
+ }
+}
+
+
/* Do spatial convolution on each mesh. */
static void *
-gal_convolve_spatial_on_thread(void *inparm)
+convolve_spatial_on_thread(void *inparam)
{
+ struct gal_tile_thread_param *tprm=(struct gal_tile_thread_param *)inparam;
+ struct spatial_params *cprm=(struct spatial_params *)(tprm->params);
+
+ size_t i;
+ gal_data_t overlap, *block=gal_tile_block(cprm->tiles);
+ size_t *parse_coords=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T,
+ 4*block->ndim);
+
+ /* Initialize the necessary overlap parameters. */
+ overlap.block=block;
+ overlap.ndim=block->ndim;
+ overlap.dsize=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, block->ndim);
+
+
+ /* Go over all the tiles given to this thread. */
+ for(i=0; tprm->indexs[i] != GAL_THREADS_NON_THRD_INDEX; ++i)
+ convolve_spatial_tile(cprm, tprm->indexs[i], parse_coords, &overlap);
+
+
+ /* Clean up, wait until all other threads finish, then return. In a
+ single thread situation, `tprm->b==NULL'. */
+ if(tprm->b) pthread_barrier_wait(tprm->b);
return NULL;
}
@@ -44,8 +266,38 @@ gal_convolve_spatial_on_thread(void *inparm)
/* Convolve a dataset with a given kernel in the spatial domain. Since
spatial convolution can be very series of tiles arranged as an array. */
-void
-gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel)
+gal_data_t *
+gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel,
+ size_t numthreads, int convoverch)
{
+ struct spatial_params params;
+ gal_data_t *out, *block=gal_tile_block(tiles);
+
+ /* Small sanity checks. */
+ if(tiles->ndim!=kernel->ndim)
+ error(EXIT_FAILURE, 0, "The number of dimensions between the kernel and "
+ "input should be the same in `gal_convolve_spatial'");
+ if( block->type!=GAL_DATA_TYPE_FLOAT32
+ || kernel->type!=GAL_DATA_TYPE_FLOAT32 )
+ error(EXIT_FAILURE, 0, "`gal_convolve_spatial' currently only works on "
+ "`float32' type input and kernel");
+
+ /* Allocate the output convolved dataset. */
+ out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, block->ndim, block->dsize,
+ block->wcs, 0, block->minmapsize, "CONVOLVED",
+ block->unit, NULL);
+ { float *f=out->array, *ff=f+out->size; do *f=NAN; while(++f<ff); }
+
+ /* Set the pointers in the parameters structure. */
+ params.out=out;
+ params.tiles=tiles;
+ params.kernel=kernel;
+ params.convoverch=convoverch;
+
+ /* Do the spatial convolution on threads. */
+ gal_tile_function_on_threads(tiles, convolve_spatial_on_thread,
+ numthreads, ¶ms);
+ /* Clean up and return the output array. */
+ return out;
}
diff --git a/lib/gnuastro/convolve.h b/lib/gnuastro/convolve.h
index 9ed035f..5751168 100644
--- a/lib/gnuastro/convolve.h
+++ b/lib/gnuastro/convolve.h
@@ -44,8 +44,9 @@ __BEGIN_C_DECLS /* From C++ preparations */
-void
-gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel);
+gal_data_t *
+gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel,
+ size_t numthreads, int convoverch);
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index 6891d4d..45d6b63 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -48,23 +48,30 @@ __BEGIN_C_DECLS /* From C++ preparations */
-
/***********************************************************************/
-/************** About block **********************/
+/************** Single tile ******************/
/***********************************************************************/
-gal_data_t *
-gal_tile_block(gal_data_t *input);
+void
+gal_tile_start_coord(gal_data_t *tile, size_t *start_coord);
void
-gal_tile_block_start_coord(gal_data_t *tile, size_t *start_coord);
+gal_tile_start_end_coord(gal_data_t *tile, size_t *start_end, int rel_block);
void *
-gal_tile_block_start_end(gal_data_t *tile, gal_data_t *work,
- size_t *start_end);
+gal_tile_start_end_ind_inclusive(gal_data_t *tile, gal_data_t *work,
+ size_t *start_end_inc);
+
+
+
+/***********************************************************************/
+/************** Allocated block **********************/
+/***********************************************************************/
+gal_data_t *
+gal_tile_block(gal_data_t *input);
size_t
gal_tile_block_increment(gal_data_t *block, size_t *tsize,
- size_t num_increment);
+ size_t num_increment, size_t *coord);
void
gal_tile_block_check_tiles(gal_data_t *tiles, char *filename,
@@ -92,6 +99,24 @@ gal_tile_all_position_two_layers(gal_data_t *input, size_t
*channel_size,
+/*********************************************************************/
+/******************** On threads ********************/
+/*********************************************************************/
+struct gal_tile_thread_param
+{
+ size_t id; /* Id of this thread. */
+ void *params; /* Input structure for higher-level settings. */
+ size_t *indexs; /* Indexes of actions to be done in this thread. */
+ pthread_barrier_t *b; /* Pointer the barrier for all threads. */
+};
+
+void
+gal_tile_function_on_threads(gal_data_t *tiles, void *(*meshfunc)(void *),
+ size_t numthreads, void *caller_params);
+
+
+
+
__END_C_DECLS /* From C++ preparations */
#endif
diff --git a/lib/tile.c b/lib/tile.c
index f1f742d..b2cad46 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/fits.h>
#include <gnuastro/tile.h>
#include <gnuastro/blank.h>
+#include <gnuastro/threads.h>
#include <gnuastro/convolve.h>
#include <gnuastro/multidim.h>
@@ -39,28 +40,18 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
-/***********************************************************************/
-/************** Allocated block of memory ******************/
-/***********************************************************************/
-/* When you are working on an array, it important to know the size of the
- allocated space in each dimension. This simple function will just follow
- the block pointer and return the `dsize' element of lowest-level
- structure. */
-gal_data_t *
-gal_tile_block(gal_data_t *input)
-{
- while(input->block!=NULL) input=input->block;
- return input;
-}
+/***********************************************************************/
+/************** Single tile ******************/
+/***********************************************************************/
/* Calculate the starting coordinates of a tile in the allocated block of
memory. */
void
-gal_tile_block_start_coord(gal_data_t *tile, size_t *start_coord)
+gal_tile_start_coord(gal_data_t *tile, size_t *start_coord)
{
size_t ind, ndim=tile->ndim;
gal_data_t *block=gal_tile_block(tile);
@@ -80,13 +71,55 @@ gal_tile_block_start_coord(gal_data_t *tile, size_t
*start_coord)
-/* Given a tile */
+
+/* Put the starting and ending (end point is not inclusive) coordinates of
+ a tile into the `start_end' array. It is assumed that a space of
+ `2*tile->ndim' has been already allocated (static or dynamic) before
+ this function is called.
+
+ `rel_block' (or relative-to-block) is only relevant when the tile has an
+ intermediate tile between it and the allocated space (like a channel,
+ see `gal_tile_all_position_two_layers'). If it doesn't (`tile->block'
+ points the allocated dataset), then the value to `rel_block' is
+ irrelevant.
+
+ However, when `tile->block' is its self a larger block and `rel_block'
+ is set to 0, then the starting and ending positions will be based on the
+ position within `tile->block', not the allocated space. */
+void
+gal_tile_start_end_coord(gal_data_t *tile, size_t *start_end, int rel_block)
+{
+ size_t start_ind, ndim=tile->ndim;
+ gal_data_t *block=gal_tile_block(tile);
+ gal_data_t *host=rel_block ? block : tile->block;
+
+ /* Get the starting index. Note that for the type we need the allocated
+ block dataset and can't rely on the tiles. */
+ start_ind=gal_data_ptr_dist(host->array, tile->array, block->type);
+
+ /* Get the coordinates of the starting point. */
+ gal_multidim_index_to_coord(start_ind, tile->ndim, host->dsize,
+ start_end);
+
+ /* Add the dimensions of the tile to the starting coordinate. Note that
+ the ending coordinates are stored immediately after the start.*/
+ gal_multidim_add_coords(start_end, tile->dsize, start_end+ndim, ndim);
+}
+
+
+
+
+
+/* Put the indexs of the first/start and last/end pixels (inclusive) in a
+ tile into the `start_end' array (that has two elements). It will then
+ return the pointer to the start of the tile in the `work' data
+ structure. */
void *
-gal_tile_block_start_end(gal_data_t *tile, gal_data_t *work,
- size_t *start_end)
+gal_tile_start_end_ind_inclusive(gal_data_t *tile, gal_data_t *work,
+ size_t *start_end_inc)
{
- size_t ndim=tile->ndim, *s, *sf;
gal_data_t *block=gal_tile_block(tile);
+ size_t ndim=tile->ndim, *s, *e, *l, *sf;
size_t *start_coord = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
size_t *end_coord = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
@@ -94,17 +127,13 @@ gal_tile_block_start_end(gal_data_t *tile, gal_data_t
*work,
/* The starting index can be found from the distance of the `tile->array'
pointer and `block->array' pointer. IMPORTANT: with the type of the
block array. */
- start_end[0]=gal_data_ptr_dist(block->array, tile->array, block->type);
+ start_end_inc[0]=gal_data_ptr_dist(block->array, tile->array, block->type);
/* To find the end index, we need to know the coordinates of the starting
point in the allocated block. */
- gal_multidim_index_to_coord(start_end[0], ndim, block->dsize, start_coord);
-
-
- /* Adding the coordinates of the starting point with the tile's
- dimensions, we will get the ending point's coordinates. */
- gal_multidim_add_coords(start_coord, tile->dsize, end_coord, ndim);
+ gal_multidim_index_to_coord(start_end_inc[0], ndim, block->dsize,
+ start_coord);
/* `end_coord' is one unit ahead of the last element in the tile in every
@@ -112,16 +141,13 @@ gal_tile_block_start_end(gal_data_t *tile, gal_data_t
*work,
value, so we get the coordinates of the last pixel in the tile
(inclusive). We will finally, increment that value by one to get to
the pixel immediately outside of the tile.*/
- sf=(s=end_coord)+ndim; do *s-=1; while(++s<sf);
+ e=end_coord;
+ l=tile->dsize;
+ sf=(s=start_coord)+ndim; do *e++ = *s + *l++ - 1; while(++s<sf);
/* Convert the (inclusive) ending point's coordinates into an index. */
- start_end[1]=gal_multidim_coord_to_index(ndim, block->dsize, end_coord);
-
-
- /* Increment the (inclusive) ending point's index by one to be
- immediately outside the tile. */
- ++start_end[1];
+ start_end_inc[1]=gal_multidim_coord_to_index(ndim, block->dsize, end_coord);
/* For a check:
@@ -131,8 +157,8 @@ gal_tile_block_start_end(gal_data_t *tile, gal_data_t *work,
start_coord[2]);
printf("end_coord: %zu, %zu, %zu\n", end_coord[0], end_coord[1],
end_coord[2]);
- printf("start_index: %zu\n", start_end[0]);
- printf("end_index: %zu\n", start_end[1]);
+ printf("start_index: %zu\n", start_end_inc[0]);
+ printf("end_index: %zu\n", start_end_inc[1]);
exit(1);
*/
@@ -141,7 +167,40 @@ gal_tile_block_start_end(gal_data_t *tile, gal_data_t
*work,
from. */
free(end_coord);
free(start_coord);
- return gal_data_ptr_increment(work->array, start_end[0], work->type);
+ return gal_data_ptr_increment(work->array, start_end_inc[0], work->type);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/***********************************************************************/
+/************** Allocated block of memory ******************/
+/***********************************************************************/
+/* When you are working on an array, it important to know the size of the
+ allocated space in each dimension. This simple function will just follow
+ the block pointer and return the `dsize' element of lowest-level
+ structure. */
+gal_data_t *
+gal_tile_block(gal_data_t *input)
+{
+ while(input->block!=NULL) input=input->block;
+ return input;
}
@@ -167,27 +226,28 @@ gal_tile_block_start_end(gal_data_t *tile, gal_data_t
*work,
the tile has been parsed), simply incrementing by `b[n-2] * b[n-1]' will
take us to the last row of
- num_increment increment
- ------------- ---------
- 0 b[n-1]: fastest dimension of the block.
- 1 Similar to previous
- . .
- . .
- t[n-2] (b[n-2] * b[n-1]) - ( (t[n-2]-1) * b[n-1] )
- t[n-2] + 1 b[n-1]
- . .
- . .
- 2 * t[n-2] b[n-2] * b[n-1]
- t[n-2]+1 b[n-1]
- . .
- . .
- t[n-3] * t[n-2] b[n-3] * b[n-2] * b[n-1]
+ num_increment coord increment
+ ------------- ----- ---------
+ 0 (...0,0,0) b[n-1]: fastest dimension of the block.
+ 1 (...0,1,0) Similar to previous
+ . . .
+ . . .
+ t[n-2] (...1,0,0) (b[n-2] * b[n-1]) - ( (t[n-2]-1) * b[n-1] )
+ t[n-2] + 1 (...1,1,0) b[n-1]
+ . . .
+ . . .
+ 2 * t[n-2] (...2,0,0) b[n-2] * b[n-1]
+ t[n-2]+1 (...2,1,0) b[n-1]
+ . . .
+ . . .
+ t[n-3] * t[n-2] (..1,0,0,0) b[n-3] * b[n-2] * b[n-1]
*/
size_t
gal_tile_block_increment(gal_data_t *block, size_t *tsize,
- size_t num_increment)
+ size_t num_increment, size_t *coord)
{
+ size_t increment;
size_t n=block->ndim;
size_t *b=block->dsize, *t=tsize;
@@ -203,20 +263,34 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
/* 1D: the increment is just the tile size. */
case 1:
- return t[0];
+ increment=t[0];
+ coord[0]+=increment;
break;
/* 2D: the increment is the block's number of fastest axis pixels. */
case 2:
- return b[1];
+ increment=b[1];
+ ++coord[0];
break;
/* Higher dimensions. */
default:
- if(num_increment % t[n-2]) return b[n-1];
- else return (b[n-2] * b[n-1]) - ( (t[n-2]-1) * b[n-1] );
+ if(num_increment % t[n-2])
+ {
+ increment=b[n-1];
+ ++coord[n-2];
+ }
+ else
+ {
+ increment=(b[n-2] * b[n-1]) - ( (t[n-2]-1) * b[n-1] );
+ ++coord[n-3];
+ coord[n-2]=coord[n-1]=0;
+ }
break;
}
+
+ /* Return the final increment value. */
+ return increment;
}
@@ -237,7 +311,8 @@ gal_tile_block_check_tiles(gal_data_t *tiles, char
*filename,
gal_data_t *tofill, *tile;
gal_data_t *block=gal_tile_block(tiles);
int32_t *p, *pf, tile_index=0, *start=NULL;
- size_t ndim=tiles->ndim, increment, start_end[2];
+ size_t ndim=tiles->ndim, increment, start_end_inc[2];
+ size_t *coord=gal_data_calloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
/***************************************************************/
/************* For a check ****************
@@ -266,7 +341,7 @@ gal_tile_block_check_tiles(gal_data_t *tiles, char
*filename,
{
/* Set the starting and ending indexs of this tile over the allocated
block. */
- start=gal_tile_block_start_end(tile, tofill, start_end);
+ start=gal_tile_start_end_ind_inclusive(tile, tofill, start_end_inc);
/* Go over the full area of this tile. The loop will stop as soon as
the incrementation will go over the last index of the tile. Note
@@ -274,7 +349,7 @@ gal_tile_block_check_tiles(gal_data_t *tiles, char
*filename,
of zero is meaningful in the calculation of the increment. */
increment=0;
num_increment=1;
- while( start_end[0] + increment < start_end[1] )
+ while( start_end_inc[0] + increment <= start_end_inc[1] )
{
/* Parse the elements in the fastest-dimension (the contiguous
patch of memory associated with this tile). */
@@ -284,7 +359,7 @@ gal_tile_block_check_tiles(gal_data_t *tiles, char
*filename,
/* Increase the increment from the start of the tile for the next
contiguous patch. */
increment += gal_tile_block_increment(block, tile->dsize,
- num_increment++);
+ num_increment++, coord);
}
/* Increment the index for the next tile. */
@@ -295,6 +370,7 @@ gal_tile_block_check_tiles(gal_data_t *tiles, char
*filename,
gal_fits_img_write(tofill, filename, NULL, program_name);
/* Clean up. */
+ free(coord);
gal_data_free(tofill);
}
@@ -559,7 +635,7 @@ gal_tile_all_position(gal_data_t *input, size_t *regular,
if(input->block)
{
start=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim);
- gal_tile_block_start_coord(input, start);
+ gal_tile_start_coord(input, start);
}
@@ -699,3 +775,98 @@ gal_tile_all_position_two_layers(gal_data_t *input, size_t
*channel_size,
gal_tile_all_position(&ch[i], tile_size, remainderfrac, &t, 1);
}
}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*********************************************************************/
+/******************** On threads ********************/
+/*********************************************************************/
+/* Run a given function on the given tiles. */
+void
+gal_tile_function_on_threads(gal_data_t *tiles, void *(*function)(void *),
+ size_t numthreads, void *caller_params)
+{
+ int err;
+ pthread_t t; /* All thread ids saved in this, not used. */
+ pthread_attr_t attr;
+ pthread_barrier_t b;
+ struct gal_tile_thread_param *prm;
+ size_t i, *indexs, thrdcols, numbarriers;
+ size_t numtiles=gal_data_num_in_ll(tiles);
+
+ /* Allocate the array of parameters structure structures. */
+ prm=malloc(numthreads*sizeof *prm);
+ if(prm==NULL)
+ {
+ fprintf(stderr, "%zu bytes could not be allocated for prm.",
+ numthreads*sizeof *prm);
+ exit(EXIT_FAILURE);
+ }
+
+ /* Distribute the actions into the threads: */
+ gal_threads_dist_in_threads(numtiles, numthreads, &indexs, &thrdcols);
+
+ /* Do the job: when only one thread is necessary, there is no need to
+ spin off one thread, just call the function directly (spinning off
+ threads is expensive). This is for the generic thread spinner
+ function, not this simple function where `numthreads' is a
+ constant. */
+ if(numthreads==1)
+ {
+ prm[0].id=0;
+ prm[0].b=NULL;
+ prm[0].indexs=indexs;
+ prm[0].params=caller_params;
+ function(&prm[0]);
+ }
+ else
+ {
+ /* Initialize the attributes. Note that this running thread
+ (that spinns off the nt threads) is also a thread, so the
+ number the barriers should be one more than the number of
+ threads spinned off. */
+ numbarriers = (numtiles<numthreads ? numtiles : numthreads) + 1;
+ gal_threads_attr_barrier_init(&attr, &b, numbarriers);
+
+ /* Spin off the threads: */
+ for(i=0;i<numthreads;++i)
+ if(indexs[i*thrdcols]!=GAL_THREADS_NON_THRD_INDEX)
+ {
+ prm[i].id=i;
+ prm[i].b=&b;
+ prm[i].params=caller_params;
+ prm[i].indexs=&indexs[i*thrdcols];
+ err=pthread_create(&t, &attr, function, &prm[i]);
+ if(err)
+ {
+ fprintf(stderr, "can't create thread %zu", i);
+ exit(EXIT_FAILURE);
+ }
+ }
+
+ /* Wait for all threads to finish and free the spaces. */
+ pthread_barrier_wait(&b);
+ pthread_attr_destroy(&attr);
+ pthread_barrier_destroy(&b);
+ }
+
+ /* Clean up. */
+ free(indexs);
+}
- [gnuastro-commits] master 18c9823 094/125: Image created on second HDU, common options to gnuastro.conf, (continued)
- [gnuastro-commits] master 18c9823 094/125: Image created on second HDU, common options to gnuastro.conf, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master db9f983 065/125: Arithmetic uses new option management, no more mask images, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master c1dba50 070/125: Column I/O options are now common options, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 5eb817a 052/125: Option values stored in argp_option, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 7e888ce 073/125: First implementation of MakeProfiles with gal_data_t, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e8d9751 087/125: MakeNoise now uses the new gal_data_t infrastructure, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e978d44 068/125: Numeric option sanity checks, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e74eb25 093/125: Changed name of Header program to Fits, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e353932 109/125: NaN magnitude will only blank the profile's area, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 96ab804 098/125: Working on tiles implemented in tile checking, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 31b3b76 099/125: Spatial convolution now using new tessellation,
Mohammad Akhlaghi <=
- [gnuastro-commits] master c0dd6db 102/125: The noedgecorrection option in spatial domain convolution, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0dfdc34 069/125: Option values directly set in library, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 95ba655 100/125: Spatial convolution with new tessellation, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0dec9ac 110/125: Allow for floating point errors in Crop comparisons, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master db4c16d 076/125: Pulled in all the recent work in master, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 64f8db5 080/125: ConvertType working except for text inputs and outputs, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 4d853b4 112/125: Old neighbors.h header removed, used new in MakeProfiles, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 9959e45 004/125: Arithmetic using new data structure, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0a2fd0e 089/125: New implementations for histogram and CFP, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 213d1c3 106/125: gal_data_copy_* functions can now work on tiles, Mohammad Akhlaghi, 2017/04/23