[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 75ba49e: Tile and Arithmetic libraries documen
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 75ba49e: Tile and Arithmetic libraries documented in manual |
|
Date: |
Wed, 3 May 2017 00:57:48 -0400 (EDT) |
branch: master
commit 75ba49ed10b30262cde6d3e132deb39df9908d4d
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Tile and Arithmetic libraries documented in manual
The `tile.h' and `arithmetic.h' headers are now fully documented in the
manual. In the process, the following changes were also made:
- The operator argument to `GAL_TILE_PARSE_OPERATE' was moved to its end,
similar to the `GAL_DIMENSION_NEIGHBOR_OP' function-like macro. It is
much more convenient (and readable) if the expression-block is at the
end so it can take its own line and with the finishing of the block the
full macro actually finishes.
- The internally necessary `gal_arithmetic_*' functions in `arithmetic.h'
were moved to a new `arithmetic-internal.h' header.
- BuildProgram was added to the README and also the generation of `Man'
pages. Its description was also updated in its `ui.c'.
- MakeCatalog will now complain with an error if the object label image
doesn't have any non-zero pixels.
- If a clump image was provided, but contained no clumps, MakeCatalog will
now behave as if no clumps image was provided.
---
README | 6 +
bin/buildprog/ui.c | 6 +-
bin/mkcatalog/ui.c | 38 +-
bin/noisechisel/clumps.c | 13 +-
bin/noisechisel/detection.c | 2 +-
bin/noisechisel/sky.c | 21 +-
bin/noisechisel/threshold.c | 8 +-
bin/noisechisel/ui.c | 15 +-
doc/Makefile.am | 17 +-
doc/gnuastro.texi | 963 +++++++++++++++++++++++++---
lib/Makefile.am | 12 +-
lib/arithmetic-binary.c | 1 +
lib/arithmetic-onlyint.c | 1 +
lib/arithmetic.c | 4 +-
lib/binary.c | 12 +-
lib/blank.c | 2 +-
lib/gnuastro-internal/arithmetic-internal.h | 74 +++
lib/gnuastro/arithmetic.h | 22 +-
lib/gnuastro/tile.h | 260 +++-----
lib/statistics.c | 21 +-
lib/tile.c | 161 ++---
21 files changed, 1209 insertions(+), 450 deletions(-)
diff --git a/README b/README
index 5f0e85e..7f8d3b7 100644
--- a/README
+++ b/README
@@ -37,6 +37,12 @@ context under categories/chapters.
and growing set of arithmetic, mathematical, and even statistical
operators (for example +, -, *, /, sqrt, log, min, average, median).
+ - BuildProgram (astbuildprog) Compile, link and run programs that depend
+ on the Gnuastro library. BuildProgram will automatically link with the
+ libraries that Gnuastro depends on, so there is no need to explicily
+ mention them every time you are compiling a Gnuastro library dependent
+ program.
+
- ConvertType (astconvertt): Convert astronomical data files (FITS or
IMH) to and from several other standard image and data formats, for
example TXT, JPEG, EPS or PDF.
diff --git a/bin/buildprog/ui.c b/bin/buildprog/ui.c
index b949b27..1814dc0 100644
--- a/bin/buildprog/ui.c
+++ b/bin/buildprog/ui.c
@@ -65,8 +65,10 @@ doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" will compile
and run a "
"depends on. Hence you do not have to worry about explicitly linking "
"with CFITSIO for example if you want to work on a FITS file, or with "
"GSL if you want to use GNU Scientific Library's functions. You can "
- "optionally use the compiler `-I' and `-L' options to use other "
- "include headers or link with other libraries respectively.\n"
+ "optionally use the standard compiler options `-I', `-L', and `-l'. "
+ "You can use these to specify non-standard directories to look for "
+ "headers, directories to look for compiled libraries, or libraries to "
+ "link with, respectively.\n"
GAL_STRINGS_MORE_HELP_INFO
/* After the list of options: */
"\v"
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index ee64156..e6db9ff 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -575,6 +575,14 @@ ui_preparations_read_keywords(struct mkcatalogparams *p)
}
+ /* If there were no objects in the input, then inform the user with an
+ error (no catalog was built). */
+ if(p->numobjects==0)
+ error(EXIT_FAILURE, 0, "no object labels (non-zero pixels) in "
+ "%s (hdu %s). To make a catalog, labeled regions must be defined",
+ objectsfile, p->objectshdu);
+
+
/* Read the keywords from the clumps image if necessary. */
if(p->clumps)
{
@@ -589,6 +597,16 @@ ui_preparations_read_keywords(struct mkcatalogparams *p)
}
+ /* If there were no clumps, then free the clumps array and set it to
+ NULL, so for the rest of the processing, MakeCatalog things that no
+ clumps image was given. */
+ if(p->numclumps==0)
+ {
+ gal_data_free(p->clumps);
+ p->clumps=NULL;
+ }
+
+
/* Clean up. */
keys[0].name=keys[1].name=NULL;
keys[0].array=keys[1].array=NULL;
@@ -705,14 +723,18 @@ ui_preparations_outnames(struct mkcatalogparams *p)
`tableformat' option. */
gal_tableintern_check_fits_format(p->cp.output, p->cp.tableformat);
- /* If a clumps image has been read, we have two outputs. */
- if(p->clumps) ui_preparations_both_names(p, p->cp.output);
- else p->objectsout=p->cp.output;
+ /* If a clumps image has been read (and there are any clumps in the
+ image), then we have two outputs. */
+ if(p->clumps)
+ ui_preparations_both_names(p, p->cp.output);
+ else
+ p->objectsout=p->cp.output;
}
else
{
/* Both clumps and object catalogs are necessary. */
- if(p->clumps) ui_preparations_both_names(p, p->inputname);
+ if(p->clumps)
+ ui_preparations_both_names(p, p->inputname);
/* We only need one objects catalog. */
else
@@ -784,15 +806,15 @@ ui_preparations(struct mkcatalogparams *p)
ui_preparations_read_inputs(p);
- /* Set the output filename(s). */
- ui_preparations_outnames(p);
-
-
/* Read the helper keywords from the inputs and if they aren't present
then calculate the necessary parameters. */
ui_preparations_read_keywords(p);
+ /* Set the output filename(s). */
+ ui_preparations_outnames(p);
+
+
/* Allocate the output columns to fill up with the program. */
columns_define_alloc(p);
diff --git a/bin/noisechisel/clumps.c b/bin/noisechisel/clumps.c
index 6f81a3b..20f9414 100644
--- a/bin/noisechisel/clumps.c
+++ b/bin/noisechisel/clumps.c
@@ -965,8 +965,8 @@ clumps_correct_sky_labels_for_check(struct
clumps_thread_params *cltprm,
/* Go over this tile and correct the values. */
- GAL_TILE_PARSE_OPERATE({if(*i>0) *i=ninds[ *(int32_t *)i ];},
- tile, NULL, 0, 1);
+ GAL_TILE_PARSE_OPERATE( tile, NULL, 0, 1,
+ {if(*i>0) *i=ninds[ *(int32_t *)i ];} );
/* Clean up. */
gal_data_free(newinds);
@@ -1066,7 +1066,7 @@ clumps_find_make_sn_table(void *in_prm)
used. */
c=0;
indarr=cltprm.indexs->array;
- GAL_TILE_PO_OISET(int32_t, int, {
+ GAL_TILE_PO_OISET(int32_t, int, tile, NULL, 0, 1, {
/* This pixel's index over all the image. */
ind = (int32_t *)i - (int32_t *)(p->clabel->array);
gal_dimension_index_to_coord(ind, ndim, dsize, icoord);
@@ -1106,7 +1106,7 @@ clumps_find_make_sn_table(void *in_prm)
}
*/
}
- }, tile, NULL, 0, 1);
+ });
/* Correct the number of indexs. */
cltprm.indexs->size=cltprm.indexs->dsize[0]=c;
@@ -1116,9 +1116,8 @@ clumps_find_make_sn_table(void *in_prm)
/* Set all river pixels to CLUMPS_INIT (to be distinguishable
from the detected regions). */
- GAL_TILE_PO_OISET(int32_t, int,
- {if(*i==CLUMPS_RIVER) *i=CLUMPS_INIT;},
- tile, NULL, 0, 1);
+ GAL_TILE_PO_OISET( int32_t, int, tile, NULL, 0, 1,
+ {if(*i==CLUMPS_RIVER) *i=CLUMPS_INIT;} );
/* For a check, the step variable will be set. */
if(clprm->step==1)
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index f528927..417b158 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -182,7 +182,7 @@ static void
detection_write_in_large(gal_data_t *tile, gal_data_t *copy)
{
uint8_t *c=copy->array;
- GAL_TILE_PARSE_OPERATE({*i=*c++;}, tile, NULL, 0, 0);
+ GAL_TILE_PARSE_OPERATE(tile, NULL, 0, 0, {*i=*c++;});
}
diff --git a/bin/noisechisel/sky.c b/bin/noisechisel/sky.c
index 92ee414..7f369e4 100644
--- a/bin/noisechisel/sky.c
+++ b/bin/noisechisel/sky.c
@@ -89,7 +89,7 @@ sky_mean_std_undetected(void *in_prm)
bintile->size=tile->size;
bintile->dsize=tile->dsize;
bintile->array=gal_tile_block_relative_to_other(tile, p->binary);
- GAL_TILE_PARSE_OPERATE({if(!*o) numsky++;}, tile, bintile, 1, 1);
+ GAL_TILE_PARSE_OPERATE(tile, bintile, 1, 1, {if(!*o) numsky++;});
/* Only continue, if the fraction of Sky values are less than the
requested fraction. */
@@ -97,14 +97,13 @@ sky_mean_std_undetected(void *in_prm)
{
/* Calculate the mean and STD over this tile. */
s=s2=0.0f;
- GAL_TILE_PARSE_OPERATE(
- {
- if(!*o)
- {
- s += *i;
- s2 += *i * *i;
- }
- }, tile, bintile, 1, 1);
+ GAL_TILE_PARSE_OPERATE(tile, bintile, 1, 1, {
+ if(!*o)
+ {
+ s += *i;
+ s2 += *i * *i;
+ }
+ } );
darr[0]=s/numsky;
darr[1]=sqrt( (s2-s*s/numsky)/numsky );
@@ -251,7 +250,7 @@ sky_subtract(struct noisechiselparams *p)
tile=&p->cp.tl.tiles[tid];
/* First subtract the Sky value from the input image. */
- GAL_TILE_PARSE_OPERATE({*i-=sky[tid];}, tile, NULL, 0, 0);
+ GAL_TILE_PARSE_OPERATE(tile, NULL, 0, 0, {*i-=sky[tid];});
/* Change to the convolved image. */
tarray=tile->array;
@@ -268,7 +267,7 @@ sky_subtract(struct noisechiselparams *p)
the threshold have to be checked for a NaN which are by
definition a very small fraction of the total pixels. And if
there are NaN pixels in the image. */
- GAL_TILE_PARSE_OPERATE({*i-=sky[tid];}, tile, NULL, 0, 0);
+ GAL_TILE_PARSE_OPERATE(tile, NULL, 0, 0, {*i-=sky[tid];});
/* Revert back to the original block. */
tile->array=tarray;
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index a555164..55bc66c 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -110,12 +110,12 @@ threshold_apply_on_thread(void *in_prm)
will help in efficiency, because the compiler can move this
check out of the loop and only check for NaN values when we
know the tile has blank pixels. */
- GAL_TILE_PO_OISET(float, uint8_t, {
+ GAL_TILE_PO_OISET(float, uint8_t, tile, p->binary, 1, 0, {
*o = ( *i > value1[tid]
? ( *i > value2[tid] ? THRESHOLD_NO_ERODE_VALUE : 1 )
: ( (tile->flag & GAL_DATA_FLAG_HASBLANK) && !(*i==*i)
? GAL_BLANK_UINT8 : 0 ) );
- }, tile, p->binary, 1, 0);
+ });
/* Revert the tile's pointers back to what they were. */
if(p->conv) { tile->array=tarray; tile->block=tblock; }
@@ -127,12 +127,12 @@ threshold_apply_on_thread(void *in_prm)
/* See the explanation above the same step in the quantile
threshold for an explanation. */
- GAL_TILE_PO_OISET(float, uint8_t, {
+ GAL_TILE_PO_OISET(float, uint8_t, tile, p->binary, 1, 0, {
*o = ( ( *i - value1[tid] > p->dthresh * value2[tid] )
? 1
: ( (tile->flag & GAL_DATA_FLAG_HASBLANK) && !(*i==*i)
? GAL_BLANK_UINT8 : 0 ) );
- }, tile, p->binary, 1, 0);
+ });
break;
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index 2db2aa8..7848a8e 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -474,7 +474,7 @@ ui_prepare_kernel(struct noisechiselparams *p)
static void
ui_prepare_tiles(struct noisechiselparams *p)
{
- gal_data_t *check, *tile;
+ gal_data_t *check;
struct gal_tile_two_layer_params *tl=&p->cp.tl, *ltl=&p->ltl;
@@ -512,18 +512,13 @@ ui_prepare_tiles(struct noisechiselparams *p)
}
- /* If the input has blank elements, then go over all the tiles and set
- their blank flag to 1 if they have blank values within
- them. Afterwards, set the use-zero flag for the tiles in any case.*/
+ /* If the input has blank elements, then set teh appropriate flag for
+ each tile.*/
if( p->input->flag & GAL_DATA_FLAG_HASBLANK )
{
- gal_tile_full_blank_flag(tl->tiles, p->cp.numthreads);
- gal_tile_full_blank_flag(ltl->tiles, p->cp.numthreads);
+ gal_tile_block_blank_flag(tl->tiles, p->cp.numthreads);
+ gal_tile_block_blank_flag(ltl->tiles, p->cp.numthreads);
}
- for(tile=tl->tiles;tile!=NULL;tile=tile->next)
- tile->flag |= GAL_DATA_FLAG_BLANK_CH;
- for(tile=ltl->tiles;tile!=NULL;tile=tile->next)
- tile->flag |= GAL_DATA_FLAG_BLANK_CH;
/* Make the tile check image if requested. */
diff --git a/doc/Makefile.am b/doc/Makefile.am
index 90f54fb..ff81d4f 100644
--- a/doc/Makefile.am
+++ b/doc/Makefile.am
@@ -84,6 +84,9 @@ dist_infognuastro_DATA = $(top_srcdir)/doc/gnuastro-figures/*
if COND_ARITHMETIC
MAYBE_ARITHMETIC_MAN = man/astarithmetic.1
endif
+if COND_BUILDPROG
+ MAYBE_BUILDPROG_MAN = man/astbuildprog.1
+endif
if COND_CONVERTT
MAYBE_CONVERTT_MAN = man/astconvertt.1
endif
@@ -123,11 +126,11 @@ endif
#if COND_TEMPLATE
# MAYBE_TEMPLATE_MAN = man/astTEMPLATE.1
#endif
-dist_man_MANS = $(MAYBE_ARITHMETIC_MAN) $(MAYBE_CONVERTT_MAN) \
- $(MAYBE_CONVOLVE_MAN) $(MAYBE_COSMICCAL_MAN) $(MAYBE_CROP_MAN) \
- $(MAYBE_FITS_MAN) $(MAYBE_WARP_MAN) $(MAYBE_MKCATALOG_MAN) \
- $(MAYBE_MKNOISE_MAN) $(MAYBE_MKPROF_MAN) $(MAYBE_NOISECHISEL_MAN) \
- $(MAYBE_STATISTICS_MAN) $(MAYBE_TABLE_MAN)
+dist_man_MANS = $(MAYBE_ARITHMETIC_MAN) $(MAYBE_BUILDPROG_MAN) \
+ $(MAYBE_CONVERTT_MAN) $(MAYBE_CONVOLVE_MAN) $(MAYBE_COSMICCAL_MAN) \
+ $(MAYBE_CROP_MAN) $(MAYBE_FITS_MAN) $(MAYBE_WARP_MAN) \
+ $(MAYBE_MKCATALOG_MAN) $(MAYBE_MKNOISE_MAN) $(MAYBE_MKPROF_MAN) \
+ $(MAYBE_NOISECHISEL_MAN) $(MAYBE_STATISTICS_MAN) $(MAYBE_TABLE_MAN)
## See if help2man is present or not. When help2man doesn't exist, we don't
@@ -149,6 +152,10 @@ man/astarithmetic.1: $(top_srcdir)/bin/arithmetic/args.h
$(ALLMANSDEP)
$(MAYBE_HELP2MAN) -n "arithmetic operations on images and numbers" \
--libtool $(toputildir)/arithmetic/astarithmetic
+man/astbuildprog.1: $(top_srcdir)/bin/buildprog/args.h $(ALLMANSDEP)
+ $(MAYBE_HELP2MAN) -n "compile, link with Gnuastro library and its
dependencies, and run a C program" \
+ --libtool $(toputildir)/convertt/astconvertt
+
man/astconvertt.1: $(top_srcdir)/bin/convertt/args.h $(ALLMANSDEP)
$(MAYBE_HELP2MAN) -n "convert known data types to each other" \
--libtool $(toputildir)/convertt/astconvertt
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 92fbe12..b34d510 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -518,7 +518,7 @@ CosmicCalculator
Library
* Review of library fundamentals:: Guide on libraries and linking.
-* BuildProgram:: Link and run source files using Gnuastro's
library.
+* BuildProgram:: Link and run source files with this library.
* Gnuastro library:: Description of all library functions.
* Library demo programs:: Demonstration for using libraries.
@@ -535,7 +535,7 @@ BuildProgram
Gnuastro library
* Configuration information:: General information about library config.
-* Library data types:: Definitions and functions for types in
Gnuastro.
+* Library data types:: Definitions and functions for types.
* Library blank values:: Blank values and functions to deal with them.
* Library data container:: General data container in Gnuastro.
* Dimensions:: Dealing with coordinates and dimensions.
@@ -574,13 +574,18 @@ Linked lists (@file{list.h})
FITS files (@file{fits.h})
-* FITS macros errors filenames:: General macros, errors and checking FITS
names
+* FITS macros errors filenames:: General macros, errors and checking names.
* CFITSIO and Gnuastro types:: Conversion between FITS and Gnuastro types.
* FITS HDUs:: Opening and getting information about HDUs.
* FITS header keywords:: Reading and writing FITS header keywords.
* FITS arrays:: Reading and writing FITS images/arrays.
* FITS tables:: Reading and writing FITS tables.
+Tessellation library (@file{tile.h})
+
+* Independent tiles:: Work on or check independent tiles.
+* Tile grid:: Cover a full dataset with non-overlapping
tiles.
+
Multithreaded programming (@file{threads.h})
* Implementation of pthread_barrier:: Some systems don't have pthread_barrier
@@ -8644,21 +8649,22 @@ condition. Three operands are required for
@command{where}. The input
format is demonstrated in this simplified example:
@example
-$ astarithmetic modify.fits cond.fits if-true.fits where
+$ astarithmetic modify.fits binary.fits if-true.fits where
@end example
-The value of any pixel in @file{out.fits} that corresponds to a non-zero
-pixel of @file{cond.fits} will be changed to the value of the same pixel in
address@hidden (or a fixed number). The third and second popped
-operands (@file{modify.fits} and @file{cond.fits} above, see @ref{Reverse
-polish notation}) have to have the same
+The value of any pixel in @file{modify.fits} that corresponds to a non-zero
+pixel of @file{binary.fits} will be changed to the value of the same pixel
+in @file{if-true.fits} (this may also be a number). The 3rd and 2nd popped
+operands (@file{modify.fits} and @file{binary.fits} respectively, see
address@hidden polish notation}) have to have the same
dimensions/size. @file{if-true.fits} can be either a number, or have the
same dimension/size as the other two.
-That @file{cond.fits} has to have @code{uint8} (or @code{unsigned char} in
-standard C) type (see @ref{Numeric data types} and the type conversion
-operators below). However, ommonly this is no problem because in most cases
-you won't be dealing with an actual FITS file of a conditional image. You
+The 2nd popped operand (@file{binary.fits}) has to have @code{uint8} (or
address@hidden char} in standard C) type (see @ref{Numeric data types}). It
+is treated as a binary dataset (with only two values: zero and non-zero,
+hence the name @code{binary.fits} in this example). However, commonly you
+won't be dealing with an actual FITS file of a condition/binary image. You
will probably define the condition in the same run based on some other
reference image and use the conditional and logical operators above to make
a true/false (or one/zero) image for you internally. For example the case
@@ -8670,15 +8676,27 @@ $ astarithmetic in.fits reference.fits 100 gt new.fits
where
In the example above, any of the @file{in.fits} pixels that has a value in
@file{reference.fits} greater than @command{100}, will be replaced with the
-corresponding pixel in @file{new.fits}. Finally, the input operands are
-read and used independently, so you can use the same file more than once as
-any of the operands.
+corresponding pixel in @file{new.fits}. Effectively the
address@hidden 100 gt} part created the condition/binary image which
+was added to the stack (in memory) and later used by @code{where}. The
+command above is thus equivalent to these two commands:
+
address@hidden
+$ astarithmetic reference.fits 100 gt --output=binary.fits
+$ astarithmetic in.fits binary.fits new.fits where
address@hidden example
+
+Finally, the input operands are read and used independently, so you can use
+the same file more than once as any of the operands.
-When the first popped operand (@file{if-true.fits}) is a single number, it
-can be a NaN value (or any blank value, depending on its type) like the
-example below (see @ref{Blank pixels}). In this example, all the pixels in
address@hidden that have a value greater than 100, will become blank
-in the natural data type of @file{in.fits}.
+When the 1st popped operand to @code{where} (@file{if-true.fits}) is a
+single number, it may be a NaN value (or any blank value, depending on its
+type) like the example below (see @ref{Blank pixels}). When the number is
+blank, it will be converted to the blank value of the type of the 3rd
+popped operand (@code{in.fits}). Hence, in the example below, all the
+pixels in @file{reference.fits} that have a value greater than 100, will
+become blank in the natural data type of @file{in.fits} (even though NaN
+values are only defined for floating point types).
@example
$ astarithmetic in.fits reference.fits 100 gt nan where
@@ -15885,7 +15903,7 @@ program} to easily define your own programs(s).
@menu
* Review of library fundamentals:: Guide on libraries and linking.
-* BuildProgram:: Link and run source files using Gnuastro's
library.
+* BuildProgram:: Link and run source files with this library.
* Gnuastro library:: Description of all library functions.
* Library demo programs:: Demonstration for using libraries.
@end menu
@@ -16609,7 +16627,7 @@ documentation will correspond to your installed version.
@menu
* Configuration information:: General information about library config.
-* Library data types:: Definitions and functions for types in
Gnuastro.
+* Library data types:: Definitions and functions for types.
* Library blank values:: Blank values and functions to deal with them.
* Library data container:: General data container in Gnuastro.
* Dimensions:: Dealing with coordinates and dimensions.
@@ -17428,70 +17446,9 @@ for more on these wonderful high-level constructs.
Pointer to the start of the complete allocated block of memory. When this
pointer is not @code{NULL}, the dataset is not treated as a contiguous
patch of memory. Rather, it is seen as covering only a portion of the
-larger patch of memory that @code{block} points to.
-
-In many contexts, it is desirable to slice the dataset into subsets or
-tiles (not necessarily overlapping). In such a way that you can work on
-each tile independently. One method would be to copy that region to a
-separate allocated space, but in many contexts this isn't necessary and
-infact can be a big burden on CPU/Memory usage. The @code{block} pointer in
address@hidden is defined for such situations: where allocation is not
-necessary. You just want to read the data or write to it independently (or
-in coordination with) other regions of the dataset. Added with parallel
-processing, this can greatly improve the time/memory consumption.
-
-See the figure below for example: assume the @code{larger} dataset is a
-contiguous block of memory that you are interpretting as a 2D array. But
-you only want to work on the smaller @code{tile} region.
-
address@hidden
- larger
- ---------------------------------
- | |
- | tile |
- | ---------- |
- | | | |
- | |_ | |
- | |*| | |
- | ---------- |
- | tile->block = larger |
- |_ |
- |*| |
- ---------------------------------
address@hidden example
-
-To use @code{gal_data_t}'s @code{block} concept, you allocate a
address@hidden *tile} which is initialized with the pointer to the first
-element in the sub-array (as its @code{array} argument). Note that this is
-not necessarily the first element in the larger array. You can set the size
-of the tile along with the initialization as you please. Recall that, when
-given a address@hidden pointer as @code{array}, @code{gal_data_initialize}
-(and thus @code{gal_data_alloc}) do not allocate any space and just uses
-the given pointer for the new @code{array} element of the
address@hidden So your @code{tile} data structure will not be pointing
-to a separately allocated space.
-
-After the allocation is done, you just point @code{tile->block} to the
address@hidden dataset which hosts the full block of memory. Gnuastro's
-library will check the @code{block} pointer of their input dataset to see
-how to deal with dimensions and increments so they can always remain within
-the tile, see @ref{Tessellation library}.
-
-Since the block structure is defined as a pointer, arbitrary levels of
-tesselation/griding are possible (@code{tile->block} may itself be a tile
-in an even larger allocated space). Therefore, just like a linked-list (see
address@hidden lists}), it is important to have the @code{block} pointer of
-the largest dataset set to @code{NULL}. Normally, you won't have to worry
-about this, because @code{gal_data_initialize} (and thus
address@hidden) will set the @code{block} element to @code{NULL} by
-default (just remember not to change it). You can then only change the
address@hidden element for the tiles you define over the allocated space.
-
-The functions in @ref{Tessellation library} are defined to help you most
-effectively use/check/define a tessellation or tiles. This approach to
-dealing with parts of a larger block was inspired from the way the GNU
-Scientific Library (GSL) does it in its ``Vectors and Matrices'' chapter.
-
+larger patch of memory that @code{block} points to. See @ref{Tessellation
+library} for a more thorough explanation and functions to help work with
+tiles that are created from this pointer.
@end table
@@ -18638,7 +18595,7 @@ functions in your code along with those discussed here.
@menu
-* FITS macros errors filenames:: General macros, errors and checking FITS
names
+* FITS macros errors filenames:: General macros, errors and checking names.
* CFITSIO and Gnuastro types:: Conversion between FITS and Gnuastro types.
* FITS HDUs:: Opening and getting information about HDUs.
* FITS header keywords:: Reading and writing FITS header keywords.
@@ -19416,9 +19373,839 @@ saying that the @code{filename} has been created.
@node Arithmetic on datasets, Tessellation library, Table input output,
Gnuastro library
@subsection Arithmetic on datasets (@file{arithmetic.h})
+When the dataset's type and other information are already known, any
+programming language (including C) provides some very good tools for
+various operations (including arithmetic operations like addition) on the
+dataset with a simple loop. However, as an author of a program, making
+assumptions about the type of data, its dimensionality and other basic
+characteristics will come with a large processing burdern.
+
+For example if you always read your data as double precision floating
+points for a simple operation like addition with an integer constant, you
+will be wasting a lot of CPU and memory when the input dataset is
address@hidden type for example (see @ref{Numeric data types}). This overhead
+may be small for small images, but as you scale your process up and work
+with hundred/thousands of files that can be very large, this overhead will
+take a significant portion of the processing power. The functions and
+macros in this section are designed precisely for this purpose: to allow
+you to do any of the defined operations on any dataset with no overhead (in
+the native type of the dataset).
+
+Gnuastro's Arithmetic program uses the functions and macros of this
+section, so please also have a look at the @ref{Arithmetic} program and in
+particular @ref{Arithmetic operators} for a better description of the
+operators discussed here.
+
+The main function of this library is @code{gal_arithmetic} that is
+described below. It can take an arbitrary number of arguments as operands
+(depending on the operator, similar to @code{printf}). Its first two
+arguments are integers specifying the flags and operator. So first we will
+review the constants for the recognized flags and operators and discuss
+them, then introduce the actual function.
+
address@hidden Macro GAL_ARITHMETIC_INPLACE
address@hidden Macro GAL_ARITHMETIC_FREE
address@hidden Macro GAL_ARITHMETIC_NUMOK
address@hidden Macro GAL_ARITHMETIC_FLAGS_ALL
address@hidden Bitwise Or
+Bit-wise flags to pass onto @code{gal_arithmetic} (see below). To pass
+multiple flags, use the bitwise-or operator, for example
address@hidden |
+GAL_ARITHMETIC_FREE}. @code{GAL_ARITHMETIC_FLAGS_ALL} is a combination of
+all flags to shorten your code if you want all flags activated. Each flag
+is described below:
address@hidden @code
address@hidden GAL_ARITHMETIC_INPLACE
+Do the operation in-place (in the input dataset, thus modifying it) to
+improve CPU and memory usage. If this flag is used, after
address@hidden finishes, the input dataset will be modified. It is
+thus useful if you have no more need for the input after the operation.
+
address@hidden GAL_ARITHMETIC_FREE
+Free (all the) input dataset(s) after the operation is done. Hence the
+inputs are no longer usable after @code{gal_arithmetic}.
+
address@hidden GAL_ARITHMETIC_NUMOK
+It is acceptable to use a number and an array together. For example if you
+want to add all the pixels in an image with a single number you can pass
+this flag to avoid having to allocate a constant array the size of the
+image (with all the pixels having the same number).
address@hidden table
address@hidden deffn
+
address@hidden Macro GAL_ARITHMETIC_OP_PLUS
address@hidden Macro GAL_ARITHMETIC_OP_MINUS
address@hidden Macro GAL_ARITHMETIC_OP_MULTIPLY
address@hidden Macro GAL_ARITHMETIC_OP_DIVIDE
address@hidden Macro GAL_ARITHMETIC_OP_LT
address@hidden Macro GAL_ARITHMETIC_OP_LE
address@hidden Macro GAL_ARITHMETIC_OP_GT
address@hidden Macro GAL_ARITHMETIC_OP_GE
address@hidden Macro GAL_ARITHMETIC_OP_EQ
address@hidden Macro GAL_ARITHMETIC_OP_NE
address@hidden Macro GAL_ARITHMETIC_OP_AND
address@hidden Macro GAL_ARITHMETIC_OP_OR
+Binary operators (requiring two operands) that accept datasets of any
+recognized type (see @ref{Numeric data types}). When @code{gal_arithmetic}
+is called with any of these operators, it expects two datasets as
+arguments. For a full description of these operators with the same name,
+see @ref{Arithmetic operators}. The first dataset/operand will be put on
+the left of the operator and the second will be put on the right. The
+output type of the first four is determined from the input types (largest
+type of the inputs). The rest (which are all conditional operators) will
+output a binary @code{uint8_t} (or @code{unsigned char}) dataset with
+values of either @code{0} (zero) or @code{1} (one).
address@hidden deffn
+
address@hidden Macro GAL_ARITHMETIC_OP_NOT
+The logical NOT operator. When @code{gal_arithmetic} is called with this
+operator, it only expects one operand (dataset), since this is a unary
+operator. The output is @code{uint8_t} (or @code{unsigned char}) dataset of
+the same size as the input. Any non-zero element in the input will be
address@hidden (zero) in the output and any @code{0} (zero) will have a value of
address@hidden (one).
address@hidden deffn
+
address@hidden Macro GAL_ARITHMETIC_OP_ISBLANK
+A unary operator with output that is @code{1} for any element in the input
+that is blank, and @code{0} for any non-blank element. When
address@hidden is called with this operator, it will only expect one
+input dataset. The output dataset will have @code{uint8_t} (or
address@hidden char}) type.
+
address@hidden with this operator is just a wrapper for the
address@hidden function of @ref{Library blank values} and this
+operator is just included for completeness in arithmetic operations. So in
+your program, it might be easier to just call @code{gal_blank_flag}.
address@hidden deffn
+
address@hidden Macro GAL_ARITHMETIC_OP_WHERE
+The three-operand @emph{where} operator thoroughly discussed in
address@hidden operators}. When @code{gal_arithmetic} is called with this
+operator, it will only expect three input datasets: the first (which is the
+same as the returned dataset) is the array that will be modified. The
+second is the condition dataset (that must have a @code{uint8_t} or
address@hidden char} type), and the third is the value to be used if
+condition is non-zero.
+
+As a result, note that the order of operands when calling
address@hidden with @code{GAL_ARITHMETIC_OP_WHERE} is the opposite
+of running Gnuastro's Arithmetic program with the @code{where} operator
+(see @ref{Arithmetic}). This is because the latter uses the reverse-Polish
+notation which isn't necessary when calling a function (see @ref{Reverse
+polish notation}).
address@hidden deffn
+
address@hidden Macro GAL_ARITHMETIC_OP_SQRT
address@hidden Macro GAL_ARITHMETIC_OP_LOG
address@hidden Macro GAL_ARITHMETIC_OP_LOG10
+Unary operator functions for calculating the square root
+(@mymath{\sqrt{i}}), @mymath{ln(i)} and @mymath{log(i)} mathematic
+operators on each element of the input dataset. The output will have the
+same type as the input, so if your inputs are integer types be careful.
+
+If you want your output to be floating point but your input is an integer
+type, you can convert the input to a floating point type with
address@hidden or
address@hidden(see @ref{Copying datasets}).
address@hidden deffn
+
address@hidden Macro GAL_ARITHMETIC_OP_MINVAL
address@hidden Macro GAL_ARITHMETIC_OP_MAXVAL
address@hidden Macro GAL_ARITHMETIC_OP_NUMVAL
address@hidden Macro GAL_ARITHMETIC_OP_SUMVAL
address@hidden Macro GAL_ARITHMETIC_OP_MEANVAL
address@hidden Macro GAL_ARITHMETIC_OP_STDVAL
address@hidden Macro GAL_ARITHMETIC_OP_MEDIANVAL
+Unary operand statistical operators that will return a single value for
+datasets of any size. These are just wrappers around similar functions in
address@hidden operations} and are included in @code{gal_arithmetic} only
+for completeness (to use easily in @ref{Arithmetic}). In your programs, it
+will probably be easier if you use those @code{gal_statistics_} functions
+directly.
address@hidden deffn
+
address@hidden Macro GAL_ARITHMETIC_OP_ABS
+Unary operand absolute-value operator.
address@hidden deffn
+
address@hidden Macro GAL_ARITHMETIC_OP_MIN
address@hidden Macro GAL_ARITHMETIC_OP_MAX
address@hidden Macro GAL_ARITHMETIC_OP_NUM
address@hidden Macro GAL_ARITHMETIC_OP_SUM
address@hidden Macro GAL_ARITHMETIC_OP_MEAN
address@hidden Macro GAL_ARITHMETIC_OP_STD
address@hidden Macro GAL_ARITHMETIC_OP_MEDIAN
+Multi-operand statistical operations. When @code{gal_arithmetic} is called
+with any of these operators, it will expect only a single operand that will
+be interpretted as a list of datasets (see @ref{List of gal_data_t}. The
+output will be a single dataset with each of its elements replaced by the
+respective statistical operation on the whole list. See the discussion
+under the @code{min} operator in @ref{Arithmetic operators}.
address@hidden deffn
+
address@hidden Macro GAL_ARITHMETIC_OP_POW
+Binary operator to-power operator. When @code{gal_arithmetic} is called
+with any of these operators, it will expect two operands: rasing the first
+by the second. This operator only accepts floating point inputs and the
+output is also floating point.
address@hidden deffn
+
address@hidden Macro GAL_ARITHMETIC_OP_BITAND
address@hidden Macro GAL_ARITHMETIC_OP_BITOR
address@hidden Macro GAL_ARITHMETIC_OP_BITXOR
address@hidden Macro GAL_ARITHMETIC_OP_BITLSH
address@hidden Macro GAL_ARITHMETIC_OP_BITRSH
address@hidden Macro GAL_ARITHMETIC_OP_MODULO
+Binary integer-only operand operators. These operators are only defined on
+integer data types. When @code{gal_arithmetic} is called with any of these
+operators, it will expect two operands: the first is put on the left of the
+operator and the second on the right. The ones starting with @code{BIT} are
+the respective bit-wise operators in C and @code{MODULO} is the
+modulo/remainder operator. For a discussion on these operators, please see
address@hidden operators}.
+
+The output type is determined from the input types and C's internal
+conversions: it is strongly recommended that both inputs have the same type
+(any integer type), otherwise the bit-wise behavior will be determined by
+your compiler.
address@hidden deffn
+
address@hidden Macro GAL_ARITHMETIC_OP_BITNOT
+The unary bit-wise NOT operator. When @code{gal_arithmetic} is called with
+any of these operators, it will expect one operand of an integer type and
+preform the bitwise-NOT operation on it. The output will have the same type
+as the input.
address@hidden deffn
+
+
address@hidden Macro GAL_ARITHMETIC_OP_TO_UINT8
address@hidden Macro GAL_ARITHMETIC_OP_TO_INT8
address@hidden Macro GAL_ARITHMETIC_OP_TO_UINT16
address@hidden Macro GAL_ARITHMETIC_OP_TO_INT16
address@hidden Macro GAL_ARITHMETIC_OP_TO_UINT32
address@hidden Macro GAL_ARITHMETIC_OP_TO_INT32
address@hidden Macro GAL_ARITHMETIC_OP_TO_UINT64
address@hidden Macro GAL_ARITHMETIC_OP_TO_INT64
address@hidden Macro GAL_ARITHMETIC_OP_TO_FLOAT32
address@hidden Macro GAL_ARITHMETIC_OP_TO_FLOAT64
+Unary type-conversion operators. When @code{gal_arithmetic} is called with
+any of these operators, it will expect one operand and convert it to the
+requested type. Note that with these operators, @code{gal_arithmetic} is
+just a wrapper over the @code{gal_data_copy_to_new_type} or
address@hidden that are discussed in @code{Copying
+datasets}. It accepts these operators only for completeness and easy usage
+in @ref{Arithmetic}. So in your programs, it might be preferrable to
+directly use those functions.
address@hidden deffn
+
+
address@hidden {gal_data_t *} gal_arithmetic (int operator, unsigned char
flags, ...)
+Do the arithmetic operation of @code{operator} on the given operands (the
+third argument and any further argument). Certain special conditions can
+also be specified with the @code{flag} operator. The acceptable values for
address@hidden are defined in the macros above.
+
address@hidden is a multi-argument function (like C's
address@hidden). In other words, the number of necessary arguments is not
+fixed and depends on the value to @code{operator}. Here are a few examples
+showing this variability:
+
address@hidden
+out_1=gal_arithmetic(GAL_ARITHMETIC_OP_LOG, 0, in_1);
+out_2=gal_arithmetic(GAL_ARITHMETIC_OP_PLUS, 0, in_1, in_2);
+out_3=gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, 0, in_1, in_2, in_3);
address@hidden example
+
+The number of necessary operands for each operator (and thus the number of
+necessary arguments to @code{gal_arithmetic}) are described above under
+each operator.
address@hidden deftypefun
+
@node Tessellation library, Bounding box, Arithmetic on datasets, Gnuastro
library
@subsection Tessellation library (@file{tile.h})
+In many contexts, it is desirable to slice the dataset into subsets or
+tiles (overlapping or not). In such a way that you can work on each tile
+independently. One method would be to copy that region to a separate
+allocated space, but in many contexts this isn't necessary and infact can
+be a big burden on CPU/Memory usage. The @code{block} pointer in Gnuastro's
address@hidden data container} is defined for such situations: where
+allocation is not necessary. You just want to read the data or write to it
+independently (or in coordination with) other regions of the dataset. Added
+with parallel processing, this can greatly improve the time/memory
+consumption.
+
+See the figure below for example: assume the @code{larger} dataset is a
+contiguous block of memory that you are interpretting as a 2D array. But
+you only want to work on the smaller @code{tile} region.
+
address@hidden
+ larger
+ ---------------------------------
+ | |
+ | tile |
+ | ---------- |
+ | | | |
+ | |_ | |
+ | |*| | |
+ | ---------- |
+ | tile->block = larger |
+ |_ |
+ |*| |
+ ---------------------------------
address@hidden example
+
+To use @code{gal_data_t}'s @code{block} concept, you allocate a
address@hidden *tile} which is initialized with the pointer to the first
+element in the sub-array (as its @code{array} argument). Note that this is
+not necessarily the first element in the larger array. You can set the size
+of the tile along with the initialization as you please. Recall that, when
+given a address@hidden pointer as @code{array}, @code{gal_data_initialize}
+(and thus @code{gal_data_alloc}) do not allocate any space and just uses
+the given pointer for the new @code{array} element of the
address@hidden So your @code{tile} data structure will not be pointing
+to a separately allocated space.
+
+After the allocation is done, you just point @code{tile->block} to the
address@hidden dataset which hosts the full block of memory. Where relevant,
+Gnuastro's library functions will check the @code{block} pointer of their
+input dataset to see how to deal with dimensions and increments so they can
+always remain within the tile. The tools introduced in this section are
+designed to help in defining and working with tiles that are created in
+this manner.
+
+Since the block structure is defined as a pointer, arbitrary levels of
+tesselation/griding are possible (@code{tile->block} may itself be a tile
+in an even larger allocated space). Therefore, just like a linked-list (see
address@hidden lists}), it is important to have the @code{block} pointer of
+the largest (allocated) dataset set to @code{NULL}. Normally, you won't
+have to worry about this, because @code{gal_data_initialize} (and thus
address@hidden) will set the @code{block} element to @code{NULL} by
+default, just remember not to change it. You can then only change the
address@hidden element for the tiles you define over the allocated space.
+
+Below, we will first review constructs for @ref{Independent tiles} and then
+define the current approach to fully tessellating a dataset (or covering
+every pixel/data-element with a non-overlapping tile grid in @ref{Tile
+grid}. This approach to dealing with parts of a larger block was inspired
+from a similarly named concept in the GNU Scientific Library (GSL), see its
+``Vectors and Matrices'' chapter for their implementation.
+
+
+
address@hidden
+* Independent tiles:: Work on or check independent tiles.
+* Tile grid:: Cover a full dataset with non-overlapping
tiles.
address@hidden menu
+
address@hidden Independent tiles, Tile grid, Tessellation library, Tessellation
library
address@hidden Independent tiles
+
+The most general application of tiles is to treat each independently, for
+example they may overlap, or they may not cover the full image. This
+section provides functions to help in checking/inspecting such tiles. In
address@hidden grid} we will discuss functions that define/work-with a tile grid
+(where the tiles don't overlap and fully cover the input
+dataset). Therefore, the functions in this section are general and can be
+used for the tiles produced by that section also.
+
address@hidden void gal_tile_start_coord (gal_data_t @code{*tile}, size_t
@code{*start_coord})
+Calculate the starting coordinates of a tile in its allocated block of
+memory and write them in the memory that @code{start_coord} points to
+(which must have @code{tile->ndim} elements).
address@hidden deftypefun
+
address@hidden void gal_tile_start_end_coord (gal_data_t @code{*tile}, size_t
@code{*start_end}, int @code{rel_block})
+Put the starting and ending (end point is not inclusive) coordinates of
address@hidden into the @code{start_end} array. It is assumed that a space of
address@hidden>ndim} has been already allocated (static or dynamic) for
address@hidden before this function is called.
+
address@hidden (or relative-to-block) is only relevant when @code{tile}
+has an intermediate tile between it and the allocated space (like a
+channel, see @code{gal_tile_full_two_layers}). If it doesn't
+(@code{tile->block} points the allocated dataset), then the value to
address@hidden is irrelevant.
+
+When @code{tile->block} is its self a larger block and @code{rel_block} is
+set to 0, then the starting and ending positions will be based on the
+position within @code{tile->block}, not the allocated space.
address@hidden deftypefun
+
address@hidden {void *} gal_tile_start_end_ind_inclusive (gal_data_t
@code{*tile}, gal_data_t @code{*work}, size_t @code{*start_end_inc})
+Put the indexs of the first/start and last/end pixels (inclusive) in a tile
+into the @code{start_end} array (that must have two elements). NOTE: this
+function stores the index of each point, not its coordinates. It will then
+return the pointer to the start of the tile in the @code{work} data
+structure (which doesn't have to be equal to @code{tile->block}.
+
+The outputs of this function are defined to make it easy to parse over an
+n-dimensional tile. For example, this function is one of the most important
+parts of the internal processing of in @code{GAL_TILE_PARSE_OPERATE}
+function-like macro that is described below.
address@hidden deftypefun
+
address@hidden {gal_data_t *} gal_tile_series_from_minmax (gal_data_t
@code{*block}, size_t @code{*minmax}, size_t @code{number})
+Construct a list of tile(s) given coordinates of the minimum and maximum of
+each tile. The minimum and maximums are assumed to be inclusive. The
+returned pointer is an allocated @code{gal_data_t} array that can later be
+freed with @code{gal_data_array_free} (see @ref{Arrays of
+datasets}). Internally, each element of the output array points to the next
+element, so the output may also be treated as a list of datasets (see
address@hidden of gal_data_t}) and passed onto the other functions described in
+this section.
+
+The array keeping the minmium and maximum coordinates for each tile must
+have the following format. So in total @code{minmax} must have
address@hidden elements.
+
address@hidden
+| min0_d0 | min0_d1 | max0_d0 | max0_d1 | ...
+ ... | minN_d0 | minN_d1 | maxN_d0 | maxN_d1 |
address@hidden example
address@hidden deftypefun
+
address@hidden {gal_data_t *} gal_tile_block (gal_data_t @code{*tile})
+Return the dataset that contains @code{tile}'s allocated block of
+memory. If tile is immediately defined as part of the allocated block, then
+this is equivalent to @code{tile->block}. However, it is possible to have
+multiple layers of tiles (where @code{tile->block} is itself a tile). So
+this function is the most generic way to get to the actual allocated
+dataset.
address@hidden deftypefun
+
address@hidden size_t gal_tile_block_increment (gal_data_t @code{*block},
size_t @code{*tsize}, size_t @code{num_increment}, size_t @code{*coord})
+Return the increment necessary to start at the next contiguous patch memory
+associated with a tile. @code{block} is the allocated block of memory and
address@hidden is the size of the tile along every dimension. If @code{coord}
+is @code{NULL}, it is ignored. Otherwise, it will contain the coordinate of
+the start of the next contiguous patch of memory.
+
+This function is intended to be used in a loop and @code{num_increment} is
+the main variable to this function. For the first time you call this
+function, it should be @code{1}. In subsequent calls (while you are parsing
+a tile), it should be increased by one.
address@hidden deftypefun
+
address@hidden {gal_data_t *} gal_tile_block_write_const_value (gal_data_t
@code{*tilevalues}, gal_data_t @code{*tilesll}, int @code{initialize})
+Write a constant value for each tile over the area it covers in an
+allocated dataset that is the size of @code{tile}'s allocated block of
+memory (found through @code{gal_tile_block} described above). The arguments
+to this function are:
+
address@hidden @code
address@hidden tilevalues
+This must be an array that has the same number of elements as the nodes in
+in @code{tilesll} and in the same order that `tilesll' elements are parsed
+(from top to bottom, see @ref{Linked lists}). As a result the
+dimensionality of this array is irrelevant, it will be parsed contiguously.
+
address@hidden tilesll
+The list of input tiles (see @ref{List of gal_data_t}). Internally, it
+might be stored as an array (for example the output of
address@hidden described above), but this function
+doesn't care, it will parse the @code{next} elements to go to the next
+tile. This function will not pop-from or free the @code{tilesll}, it will
+only parse it from start to end.
+
address@hidden initialize
+Initialize the allocated space with blank values before writing in the
+constant values. This can be useful when the tiles don't cover the full
+allocated block.
address@hidden table
address@hidden deftypefun
+
address@hidden {gal_data_t *} gal_tile_block_check_tiles (gal_data_t
@code{*tilesll})
+Make a copy of the memory block and fill it with the index of each tile in
address@hidden (counting from 0). The non-filled areas will have blank
+values. The output dataset will have a type of @code{GAL_TYPE_INT32} (see
address@hidden data types}).
+
+This function can be used when you want to check the coverage of each tile
+over the allocated block of memory. It is just a wrapper over the
address@hidden
address@hidden deftypefun
+
address@hidden {void *} gal_tile_block_relative_to_other (gal_data_t
@code{*tile}, gal_data_t @code{*other})
+Return the pointer corresponding to the start of the region covered by
address@hidden over the @code{other} dataset. See the examples in
address@hidden for some example applications of this
+function.
address@hidden deftypefun
+
+
address@hidden void gal_tile_block_blank_flag (gal_data_t @code{*tilell},
size_t @code{numthreads})
+Check if each tile in the list has blank values and update its @code{flag}
+to mark this check and its result (see @ref{Generic data container}). The
+operation will be done on @code{numthreads} threads.
address@hidden deftypefun
+
+
address@hidden {Function-like macro} GAL_TILE_PARSE_OPERATE (@code{IN},
@code{OTHER}, @code{PARSE_OTHER}, @code{CHECK_BLANK}, @code{OP})
+Parse over @code{IN} (which can be a tile or a fully allocated block of
+memory) and do the @code{OP} operation on it (can be any combination of C
+expresssions). If @code{OTHER!=NULL}, this macro will allow access to its
+element(s) and it can optionally be parsed while parsing over
address@hidden You can use this to parse the same region over two arrays. The
+input arguments are (the expected types of each argument are also written
+here):
+
address@hidden @code
address@hidden IN (gal_data_t)
+Input dataset, this can be a tile or an allocated block of memory.
+
address@hidden OTHER (gal_data_t)
+Dataset (@code{gal_data_t}) to parse along with @code{IN}. It can be
address@hidden In that case, @code{o} (see description of @code{OP} below)
+will be @code{NULL} and should not be used. If @code{PARSE_OTHER} is zero,
+the size of this dataset is irrelevant. Otherwise, it has to have the same
+size as the allocated block of @code{IN}.
+
address@hidden PARSE_OTHER (int)
+Parse the other dataset along with the input. When this is non-zero and
address@hidden, then the @code{o} pointer will be incremented to cover
+the @code{OTHER} tile at the same rate as @code{i}, see description of
address@hidden for @code{i} and @code{o}.
+
address@hidden CHECK_BLANK (int)
+If it is non-zero, then the input will be checked for blank values and
address@hidden will only be called when we are not on a blank element.
+
address@hidden OP
+Operator: this can be any number of C experssions. This macro is going to
+define a @code{itype *i} variable which will increment over each element of
+the input array/tile. @code{itype} will be replaced with the C type that
+corresponds to the type of @code{INPUT}. As an example, if @code{INPUT}'s
+type is @code{GAL_DATA_UINT16} or @code{GAL_DATA_FLOAT32}, @code{i} will be
+defined as @code{uint16} or @code{float} respectively.
+
+This function-like macro will also define an @code{otype *o} which you can
+use to access an element of the @code{OTHER} dataset (if
address@hidden). @code{o} will correspond to the type of @code{OTHER}
+(similar to @code{itype} and @code{INPUT} discussed above). If
address@hidden is non-zero, then @code{o} will also be incremented to
+the same index element but in the other array. You can use these along with
+any other variable you define before this macro to process the input and
+store it in the other.
+
+All variables within this function-like macro begin with @code{tpo_} except
+for the three variables listed above. So as long as you don't start the
+names of your variables with this prefix everything will be fine. Note that
address@hidden (and possibly @code{o}) will be incremented once by this
+function-like macro, so don't increment them within @code{OP}. As a
+summary, the three variables you have access to within @code{OP} are:
+
address@hidden @code
address@hidden i
+Pointer to the element of @code{INPUT} that is being parsed with the proper
+type.
+
address@hidden o
+Pointer to the element of @code{OTHER} that is being parsed with the proper
+type. @code{o} can only be used if @code{OTHER!=NULL} and it will be
+parsed/incremented if @code{PARSE_OTHER} is non-zero.
+
address@hidden b
+Blank value in the type of @code{INPUT}.
address@hidden table
address@hidden table
+
+You can use a given tile (@code{tile} on a dataset that it was not
+initialized with (but has the same size, let's call it @code{new}) with the
+following steps:
+
address@hidden
+void *tarray;
+gal_data_t *tblock;
+
+/* `tile->block' must corrected AFTER `tile->array' */
+tarray = tile->array;
+tblock = tile->block;
+tile->array = gal_tile_block_relative_to_other(tile, new);
+tile->block = new;
+
+/* Parse and operate over this region of the `new' dataset. */
+GAL_TILE_PARSE_OPERATE(tile, NULL, 0, 0, @{
+ YOUR_PROCESSING;
+ @});
+
+/* Reset `tile->block' and `tile->array'. */
+tile->array=tarray;
+tile->block=tblock;
address@hidden example
+
+You can work on the same region of another block in one run of this
+function-like macro. To do that, you can make a fake tile and pass that as
+the @code{OTHER} argument. Below is a demonstration, @code{tile} is the
+actual tile that you start with and @code{new} is the other block of
+allocated memory.
+
address@hidden
+size_t zero=0;
+gal_data_t *faketile;
+
+/* Allocate the fake tile, these can be done outside a loop
+ * (over many tiles). */
+faketile=gal_data_alloc(NULL, new->type, 1, &zero,
+ NULL, 0, -1, NULL, NULL, NULL);
+free(faketile->array); /* To keep things clean. */
+free(faketile->dsize); /* To keep things clean. */
+faketile->block = new;
+faketile->ndim = new->ndim;
+
+/* These can be done in a loop (over many tiles). */
+faketile->size = tile->size;
+faketile->dsize = tile->dsize;
+faketile->array = gal_tile_block_relative_to_other(tile, new);
+
+// Do your processing....
+GAL_TILE_PARSE_OPERATE(tile, faketile, 1, 1, @{
+ YOUR_PROCESSING_EXPRESSIONS;
+ @});
+
+// Clean up.
+bintile->array=NULL;
+bintile->dsize=NULL;
+gal_data_free(bintile);
address@hidden example
address@hidden deffn
+
+
+
address@hidden Tile grid, , Independent tiles, Tessellation library
address@hidden Tile grid
+
+One very useful application of tiles is to completely cover an input
+dataset with tiles. Such that you know every pixel/data-element of the
+input image is covered by only one tile. The constructs in this section
+allow easy definition of such a tile structure. They will create lists of
+tiles that are also usable by the general tools discussed in
address@hidden tiles}.
+
+As discussed in @ref{Tessellation}, (mainly raw) astronomical images will
+mostly require two layers of tessellation, one for amplifier channels which
+all have the same size and another (smaller tile-size) tessellation over
+each channel. Hence, in this section we define a general structure to keep
+the main parameters of this two-layer tessellation and help in benefiting
+from it.
+
address@hidden {Type (C @code{struct})} gal_tile_two_layer_params
+This is the general structure to keep all the necessary parameters for a
+two-layer tessellation.
+
address@hidden
+struct gal_tile_two_layer_params
address@hidden
+ /* Inputs */
+ size_t *tilesize; /*******************************/
+ size_t *numchannels; /* These parameters have to be */
+ float remainderfrac; /* filled manually before */
+ uint8_t workoverch; /* calling the functions in */
+ uint8_t checktiles; /* this section. */
+ uint8_t oneelempertile; /*******************************/
+
+ /* Internal parameters. */
+ size_t ndim;
+ size_t tottiles;
+ size_t tottilesinch;
+ size_t totchannels;
+ size_t *channelsize;
+ size_t *numtiles;
+ size_t *numtilesinch;
+ char *tilecheckname;
+ size_t *permutation;
+ size_t *firsttsize;
+
+ /* Tile and channel arrays (which are also lists). */
+ gal_data_t *tiles;
+ gal_data_t *channels;
address@hidden;
address@hidden example
address@hidden deftp
+
address@hidden {size_t *} gal_tile_full (gal_data_t @code{*input}, size_t
@code{*regular}, float @code{remainderfrac}, gal_data_t @code{**out}, size_t
@code{multiple}, size_t @code{**firsttsize})
+Cover the full dataset with (mostly) identical tiles and return the number
+of tiles created along each dimension. The regular tile size (along each
+dimension) is determined from the @code{regular} array. If @code{input}'s
+size is not an exact multiple of @code{regular} for each dimension, then
+the tiles touching the edges in that dimension will have a different size
+to fully cover every element of the input (depending on
address@hidden).
+
+The output is an array with the same dimensions as @code{input} which
+contains the number of tiles along each dimension. See @ref{Tessellation}
+for a description of its application in Gnuastro's programs and
address@hidden, just note that this function defines only onelayer of
+tiles.
+
+This is a low-level function (independent of the
address@hidden strcuture defined above). If you want a
+two-layer tessellation, directly call @code{gal_tile_full_two_layers} that
+is described below. The input arguments to this function are:
+
address@hidden @code
address@hidden input
+The main dataset (allocated block) which you want to create a tessellation
+over (only used for its sizes). So @code{input} may be a tile also.
+
address@hidden regular
+The the size of the regular tiles along each of the input's dimensions. So
+it must have the same number of elements as the dimensions of @code{input}
+(or @code{input->ndim}).
+
address@hidden remainderfrac
+The significant fraction of the remainder space to see if it should be
+split into two and put on both sides of a dimension or not. This is thus
+only relevant @code{input} length along a dimension isn't an exact multiple
+of the regular tile size along that dimension. See @ref{Tessellation} for a
+more thorough discussion.
+
address@hidden out
+Pointer to the array of data structures that will keep all the tiles (see
address@hidden of datasets}). If @code{*out==NULL}, then the necessary space
+to keep all the tiles will be allocated. If not, then all the tile
+information will be filled from the dataset that @code{*out} points to, see
address@hidden for more.
+
address@hidden multiple
+When @code{*out==NULL} (and thus will be allocated by this function),
+allocate space for @code{multiple} times the number of tiles needed. This
+can be very useful when you have several more identically sized
address@hidden, and you want all their tiles to be allocated (and thus
+indexed) together, even though they have different @code{block} datasets
+(that then link to one allocated space). See the definition of channels in
address@hidden and @code{gal_tile_full_two_layers} below.
+
address@hidden firsttsize
+The size of the first tile along every dimension. This is only different
+from the regular tile size when @code{regular} is not an exact multiple of
address@hidden's length along every dimension. This array is allocated
+internally by this function.
address@hidden table
address@hidden deftypefun
+
+
address@hidden void gal_tile_full_sanity_check (char @code{*filename}, char
@code{*hdu}, gal_data_t @code{*input}, struct gal_tile_two_layer_params
@code{*tl})
+Make sure that the input parameters (in @code{tl}, short for two-layer)
+correspond to the input dataset. @code{filename} and @code{hdu} are only
+required for error messages. Also, allocate and fill the
address@hidden>channelsize} array.
address@hidden deftypefun
+
address@hidden void gal_tile_full_two_layers (gal_data_t @code{*input}, struct
gal_tile_two_layer_params @code{*tl})
+Create the two layered tessellation in @code{tl}. The general set of steps
+you need to take to define the two-layered tessellation over an image can
+be seen in the example code below.
+
address@hidden
+gal_data_t *input;
+struct gal_tile_two_layer_params tl;
+char *filename="input.fits", *hdu="1";
+
+/* Set all the inputs shown in the structure definition. */
+...
+
+/* Read the input dataset. */
+input=gal_fits_img_read(filename, hdu, -1);
+
+/* Do a sanity check and preparations. */
+gal_tile_full_sanity_check(filename, hdu, input, &tl);
+
+/* Build the two-layer tessellation*/
+gal_tile_full_two_layers(input, &tl);
+
+/* `tl.tiles' and `tl.channels' are now a lists of tiles.*/
address@hidden example
address@hidden deftypefun
+
+
address@hidden void gal_tile_full_permutation (struct gal_tile_two_layer_params
@code{*tl})
+Make a permutation to allow the conversion of tile location in memory to
+its location in the full input dataset and put it in
address@hidden>permutation}. If a permutation has already been defined for the
+tessellation, this function will not do anything. If permutation won't be
+necessary (there is only one channel or one dimension), then this function
+will not do anything (@code{tl->permutation} must have been initialized to
address@hidden).
+
+When there is only one channel OR one dimension, the tiles are allocated in
+memory in the same order that they represent the input data. However, to
+make channel-independent processing possible in a generic way, the tiles of
+each channel are allocated contiguously. So, when there is more than one
+channel AND more than one dimension, the index of the tile does not
+correspond to its position in the grid covering the input dataset.
+
+The example below may help clarify: assume you have a 6x6 tessellation with
+two channels in the horizontal and one in the vertical. On the left you can
+see how the tile IDs correspond to the input dataset. NOTE how `03' is on
+the second row, not on the first after `02'. On the right, you can see how
+the tiles are stored in memory (and shown if you simply write the array
+into a FITS file for example).
+
address@hidden
+ Corresponding to input In memory
+ ---------------------- --------------
+ 15 16 17 33 34 35 30 31 32 33 34 35
+ 12 13 14 30 31 32 24 25 26 27 28 29
+ 09 10 11 27 28 29 18 19 20 21 22 23
+ 06 07 08 24 25 26 <-- 12 13 14 15 16 17
+ 03 04 05 21 22 23 06 07 08 09 10 11
+ 00 01 02 18 19 20 00 01 02 03 04 05
address@hidden example
+
+As a result, if your values are stored in same order as the tiles, and you
+want them in over-all memory (for example to save as a FITS file), you need
+to permute the values:
+
address@hidden
+gal_permutation_apply(values, tl->permutation);
address@hidden example
+
+If you have values over-all and you want them in tile-order, you can apply
+the inverse permutation:
+
address@hidden
+gal_permutation_apply_inverse(values, tl->permutation);
address@hidden example
+
+Recall that this is the definition of permutation in this context:
+
address@hidden
+permute: IN_ALL[ i ] = IN_MEMORY[ perm[i] ]
+inverse: IN_ALL[ perm[i] ] = IN_MEMORY[ i ]
address@hidden example
address@hidden deftypefun
+
address@hidden void gal_tile_full_values_write (gal_data_t @code{*tilevalues},
struct gal_tile_two_layer_params @code{*tl}, char @code{*filename},
gal_fits_list_key_t @code{*keys}, char @code{*program_string})
+Write one value for each tile into a file. It is important to note that the
+values in @code{tilevalues} must be ordered in the same manner as the
+tiles, so @code{tilevalues->array[i]} is the value that should be given to
address@hidden>tiles[i]}. The @code{tl->permutation} array must have been
+initialized before calling this function with
address@hidden
address@hidden deftypefun
+
address@hidden {gal_data_t *} gal_tile_full_values_smooth (gal_data_t
@code{*tilevalues}, struct gal_tile_two_layer_params @code{*tl}, size_t
@code{width}, size_t @code{numthreads})
+Smooth the given values with a flat kernel of the given @code{width}. This
+cannot be done manually because if @code{tl->workoverch==0}, tiles in
+different channels must not be mixed/smoothed. Also the tiles are
+contiguous within the channel, not within the image, see the description
+under @code{gal_tile_full_permutation}.
address@hidden deftypefun
+
address@hidden size_t gal_tile_full_id_from_coord (struct
gal_tile_two_layer_params @code{*tl}, size_t @code{*coord})
+Return the ID of the tile that corresponds to the coordinates
address@hidden Having this ID, you can use the @code{tl->tiles} array to get
+to the proper tile or read/write a value into an array that has one value
+per tile.
address@hidden deftypefun
+
address@hidden void gal_tile_full_free_contents (struct
gal_tile_two_layer_params @code{*tl})
+Free all the allocated arrays within @code{tl}.
address@hidden deftypefun
+
+
@node Bounding box, Git wrappers, Tessellation library, Gnuastro library
@subsection Bounding box (@file{box.h})
diff --git a/lib/Makefile.am b/lib/Makefile.am
index ef281ea..d68448a 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -85,12 +85,12 @@ pkginclude_HEADERS = gnuastro/config.h
$(headersdir)/arithmetic.h \
# will not distributed, so we need to explicitly tell Automake to
# distribute them here.
internaldir=$(top_srcdir)/lib/gnuastro-internal
-EXTRA_DIST = gnuastro.pc.in $(headersdir)/README $(internaldir)/README \
- $(internaldir)/arithmetic-binary.h $(internaldir)/arithmetic-onlyint.h \
- $(internaldir)/checkset.h $(internaldir)/commonopts.h \
- $(internaldir)/config.h.in $(internaldir)/fixedstringmacros.h \
- $(internaldir)/options.h $(internaldir)/tableintern.h \
- $(internaldir)/timing.h
+EXTRA_DIST = gnuastro.pc.in $(headersdir)/README $(internaldir)/README \
+ $(internaldir)/arithmetic-binary.h $(internaldir)/arithmetic-internal.h \
+ $(internaldir)/arithmetic-onlyint.h $(internaldir)/checkset.h \
+ $(internaldir)/commonopts.h $(internaldir)/config.h.in \
+ $(internaldir)/fixedstringmacros.h $(internaldir)/options.h \
+ $(internaldir)/tableintern.h $(internaldir)/timing.h
diff --git a/lib/arithmetic-binary.c b/lib/arithmetic-binary.c
index 886537e..651c023 100644
--- a/lib/arithmetic-binary.c
+++ b/lib/arithmetic-binary.c
@@ -29,6 +29,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/arithmetic.h>
+#include <gnuastro-internal/arithmetic-internal.h>
diff --git a/lib/arithmetic-onlyint.c b/lib/arithmetic-onlyint.c
index 53db601..564c849 100644
--- a/lib/arithmetic-onlyint.c
+++ b/lib/arithmetic-onlyint.c
@@ -30,6 +30,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/arithmetic.h>
#include <gnuastro-internal/arithmetic-onlyint.h>
+#include <gnuastro-internal/arithmetic-internal.h>
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 969a9b0..d73e189 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -35,7 +35,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro-internal/arithmetic-binary.h>
#include <gnuastro-internal/arithmetic-onlyint.h>
-
+#include <gnuastro-internal/arithmetic-internal.h>
@@ -1429,7 +1429,7 @@ gal_arithmetic(int operator, unsigned char flags, ...)
break;
case GAL_ARITHMETIC_OP_WHERE:
- d1 = va_arg(va, gal_data_t *); /* Output value/array. */
+ d1 = va_arg(va, gal_data_t *); /* To modify value/array. */
d2 = va_arg(va, gal_data_t *); /* Condition (unsigned char). */
d3 = va_arg(va, gal_data_t *); /* If true value/array. */
arithmetic_where(flags, d1, d2, d3);
diff --git a/lib/binary.c b/lib/binary.c
index 8ac04af..ee3ced4 100644
--- a/lib/binary.c
+++ b/lib/binary.c
@@ -626,8 +626,8 @@ binary_make_padded_inverse(gal_data_t *input, gal_data_t
**outtile)
/* Fill the central regions. */
in=input->array;
- GAL_TILE_PARSE_OPERATE({*i = *in==GAL_BLANK_UINT8 ? *in : !*in; ++in;},
- tile, NULL, 0, 0);
+ GAL_TILE_PARSE_OPERATE( tile, NULL, 0, 0,
+ {*i = *in==GAL_BLANK_UINT8 ? *in : !*in; ++in;} );
/* Clean up and return. */
@@ -694,12 +694,12 @@ gal_binary_fill_holes(gal_data_t *input)
/* The type of the tile is already known (it is `int32_t') and we have no
output, so we'll just put `int' as a place-holder. In this way we can
- avoid the switch statement of GAL_TILE_PARSE_OPERATE, and directly
- use the workhorse macro `GAL_TILE_PO_OISET'. */
- GAL_TILE_PO_OISET(int32_t, int, {
+ avoid the switch statement of GAL_TILE_PARSE_OPERATE, and directly use
+ the workhorse macro `GAL_TILE_PO_OISET'. */
+ GAL_TILE_PO_OISET(int32_t, int, tile, NULL, 0, 0, {
*in = *i>1 && *i!=GAL_BLANK_INT32 ? 1 : *in;
++in;
- }, tile, NULL, 0, 0);
+ });
/* Clean up and return. */
diff --git a/lib/blank.c b/lib/blank.c
index 80ed2e2..3a21971 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -105,7 +105,7 @@ gal_blank_alloc_write(uint8_t type)
void
gal_blank_initialize(gal_data_t *input)
{
- GAL_TILE_PARSE_OPERATE({*i=b;}, input, NULL, 0, 0);
+ GAL_TILE_PARSE_OPERATE(input, NULL, 0, 0, {*i=b;});
}
diff --git a/lib/gnuastro-internal/arithmetic-internal.h
b/lib/gnuastro-internal/arithmetic-internal.h
new file mode 100644
index 0000000..34578a9
--- /dev/null
+++ b/lib/gnuastro-internal/arithmetic-internal.h
@@ -0,0 +1,74 @@
+/*********************************************************************
+Arithmetic -- Preform arithmetic operations on datasets.
+This is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+ Mohammad Akhlaghi <address@hidden>
+Contributing author(s):
+Copyright (C) 2017, Free Software Foundation, Inc.
+
+Gnuastro is free software: you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation, either version 3 of the License, or (at your
+option) any later version.
+
+Gnuastro is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
+**********************************************************************/
+#ifndef __GAL_ARITHMETIC_INTERNAL_H__
+#define __GAL_ARITHMETIC_INTERNAL_H__
+
+/* Include other headers if necessary here. Note that other header files
+ must be included before the C++ preparations below */
+#include <gnuastro/data.h>
+
+
+/* When we are within Gnuastro's building process, `IN_GNUASTRO_BUILD' is
+ defined. In the build process, installation information (in particular
+ `GAL_CONFIG_ARITH_CHAR' and the rest of the types that we needed in the
+ arithmetic function) is kept in `config.h'. When building a user's
+ programs, this information is kept in `gnuastro/config.h'. Note that all
+ `.c' files must start with the inclusion of `config.h' and that
+ `gnuastro/config.h' is only created at installation time (not present
+ during the building of Gnuastro).*/
+#ifndef IN_GNUASTRO_BUILD
+#include <gnuastro/config.h>
+#endif
+
+
+/* C++ Preparations */
+#undef __BEGIN_C_DECLS
+#undef __END_C_DECLS
+#ifdef __cplusplus
+# define __BEGIN_C_DECLS extern "C" {
+# define __END_C_DECLS }
+#else
+# define __BEGIN_C_DECLS /* empty */
+# define __END_C_DECLS /* empty */
+#endif
+/* End of C++ preparations */
+
+
+
+/* Actual header contants (the above were for the Pre-processor). */
+__BEGIN_C_DECLS /* From C++ preparations */
+
+int
+gal_arithmetic_binary_out_type(int operator, gal_data_t *l, gal_data_t *r);
+
+char *
+gal_arithmetic_operator_string(int operator);
+
+gal_data_t *
+gal_arithmetic_convert_to_compiled_type(gal_data_t *in, unsigned char flags);
+
+
+
+__END_C_DECLS /* From C++ preparations */
+
+#endif /* __GAL_ARITHMETIC_H__ */
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index 23b513f..5d50518 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -1,11 +1,11 @@
/*********************************************************************
-data -- Structure and functions to represent/work with data
+Arithmetic -- Preform arithmetic operations on datasets.
This is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi <address@hidden>
Contributing author(s):
-Copyright (C) 2015, Free Software Foundation, Inc.
+Copyright (C) 2017, Free Software Foundation, Inc.
Gnuastro is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
@@ -61,9 +61,9 @@ __BEGIN_C_DECLS /* From C++ preparations */
/* Arithmetic flags. */
-#define GAL_ARITHMETIC_INPLACE 1
-#define GAL_ARITHMETIC_FREE 2
-#define GAL_ARITHMETIC_NUMOK 4
+#define GAL_ARITHMETIC_INPLACE 0x1
+#define GAL_ARITHMETIC_FREE 0x2
+#define GAL_ARITHMETIC_NUMOK 0x4
#define GAL_ARITHMETIC_FLAGS_ALL ( GAL_ARITHMETIC_INPLACE \
| GAL_ARITHMETIC_FREE \
@@ -136,23 +136,11 @@ enum gal_arithmetic_operators
-
-int
-gal_arithmetic_binary_out_type(int operator, gal_data_t *l, gal_data_t *r);
-
-char *
-gal_arithmetic_operator_string(int operator);
-
-gal_data_t *
-gal_arithmetic_convert_to_compiled_type(gal_data_t *in, unsigned char flags);
-
gal_data_t *
gal_arithmetic(int operator, unsigned char flags, ...);
-
-
__END_C_DECLS /* From C++ preparations */
#endif /* __GAL_ARITHMETIC_H__ */
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index ef4d882..0ae7af4 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -71,13 +71,11 @@ gal_tile_series_from_minmax(gal_data_t *block, size_t
*minmax, size_t number);
-
-
/***********************************************************************/
/************** Allocated block **********************/
/***********************************************************************/
gal_data_t *
-gal_tile_block(gal_data_t *input);
+gal_tile_block(gal_data_t *tile);
size_t
gal_tile_block_increment(gal_data_t *block, size_t *tsize,
@@ -93,7 +91,8 @@ gal_tile_block_check_tiles(gal_data_t *tiles);
void *
gal_tile_block_relative_to_other(gal_data_t *tile, gal_data_t *other);
-
+void
+gal_tile_block_blank_flag(gal_data_t *tile_ll, size_t numthreads);
@@ -173,45 +172,45 @@ gal_tile_full_free_contents(struct
gal_tile_two_layer_params *tl);
/***********************************************************************/
/************** Function-like macros ******************/
/***********************************************************************/
-
-/* Useful when the input and output types are already known. We want this
- to be self-sufficient (and be possible to call it independent of
- `GAL_TILE_PARSE_OPERATE'), so some variables (basic definitions) are
- re-defined here are in `GAL_TILE_PARSE_OPERATE'.*/
-#define GAL_TILE_PO_OISET(IT, OT, OP, IN, OUT, PARSE_OUT, CHECK_BLANK) { \
- int parse_out=(OUT && PARSE_OUT); \
- size_t increment=0, num_increment=1; \
- gal_data_t *out_w=OUT; /* `OUT' may be NULL. */ \
- gal_data_t *iblock = gal_tile_block(IN); \
- IT b=0, *st=NULL, *i=IN->array, *f=i+IN->size; \
- OT *ost=NULL, *o = out_w ? out_w->array : NULL; \
- gal_data_t *oblock = OUT ? gal_tile_block(OUT) : NULL; \
- int hasblank = CHECK_BLANK ? gal_blank_present(IN, 0) : 0; \
- size_t s_e_ind[2]={0,iblock->size-1}; /* -1: this is INCLUSIVE */ \
+/* Useful when the input and other types are already known. We want this to
+ be self-sufficient (and be possible to call it independent of
+ `GAL_TILE_PARSE_OPERATE'), so some variables (basic definitions) that
+ are already defined in `GAL_TILE_PARSE_OPERATE' re-defined here. */
+#define GAL_TILE_PO_OISET(IT, OT, IN, OTHER, PARSE_OTHER, CHECK_BLANK, OP) { \
+ int tpo_parse_other=(OTHER && PARSE_OTHER); \
+ size_t tpo_increment=0, tpo_num_increment=1; \
+ gal_data_t *tpo_other_w=OTHER; /* `OTHER' may be NULL. */ \
+ gal_data_t *tpo_iblock = gal_tile_block(IN); \
+ IT b, *tpo_st=NULL, *i=IN->array, *tpo_f=i+IN->size; \
+ OT *tpo_ost=NULL, *o = tpo_other_w ? tpo_other_w->array : NULL; \
+ gal_data_t *tpo_oblock = OTHER ? gal_tile_block(OTHER) : NULL; \
+ int tpo_hasblank = CHECK_BLANK ? gal_blank_present(IN, 0) : 0; \
+ size_t tpo_s_e_i[2]={0,tpo_iblock->size-1}; /* -1: this is INCLUSIVE */ \
\
- /* Write the blank value for the input type into `b'). Note that */ \
- /* a tile doesn't necessarily have to have a type. */ \
- if(hasblank) gal_blank_write(&b, iblock->type); \
+ /* Write the blank value for the input type into `b'. */ \
+ gal_blank_write(&b, tpo_iblock->type); \
\
/* If this is a tile, not a full block. */ \
- if(IN!=iblock) \
+ if(IN!=tpo_iblock) \
{ \
- st=gal_tile_start_end_ind_inclusive(IN, iblock, s_e_ind); \
- if(parse_out) \
- ost = (OT *)(oblock->array) + ( st - (IT *)(iblock->array) ); \
+ tpo_st = gal_tile_start_end_ind_inclusive(IN, tpo_iblock, \
+ tpo_s_e_i); \
+ if( tpo_parse_other ) \
+ tpo_ost = ( (OT *)(tpo_oblock->array) \
+ + ( tpo_st - (IT *)(tpo_iblock->array) ) ); \
} \
\
/* Go over contiguous patches of memory. */ \
- while( s_e_ind[0] + increment <= s_e_ind[1] ) \
+ while( tpo_s_e_i[0] + tpo_increment <= tpo_s_e_i[1] ) \
{ \
\
/* If we are on a tile, reset `i' and `f'. Also, if there */ \
- /* is more than one element in `OUT', then set that. Note */ \
+ /* is more than one element in `OTHER', then set that. Note */ \
/* that it is upto the caller to increment `o' in `OP'. */ \
- if(IN!=iblock) \
+ if(IN!=tpo_iblock) \
{ \
- f = ( i = st + increment ) + IN->dsize[IN->ndim-1]; \
- if(parse_out) o = ost + increment; \
+ tpo_f = ( i = tpo_st + tpo_increment ) + IN->dsize[IN->ndim-1]; \
+ if(tpo_parse_other) o = tpo_ost + tpo_increment; \
} \
\
/* Do the operation depending the nature of the blank value. */ \
@@ -223,23 +222,26 @@ gal_tile_full_free_contents(struct
gal_tile_two_layer_params *tl);
/* case, the only way to check if a data element is blank or */ \
/* not is to see if the value of each element is equal to */ \
/* itself or not. */ \
- if(hasblank) \
+ if(tpo_hasblank) \
{ \
- if(b==b) do{if(*i!=b) {OP;}if(parse_out) ++o;} while(++i<f); \
- else do{if(*i==*i){OP;}if(parse_out) ++o;} while(++i<f); \
+ if(b==b) \
+ do{if(*i!=b) {OP;} if(tpo_parse_other) ++o;} while(++i<tpo_f);\
+ else \
+ do{if(*i==*i) {OP;} if(tpo_parse_other) ++o;} while(++i<tpo_f);\
} \
- else do {{OP;}if(parse_out) ++o;} while(++i<f); \
+ else \
+ do { {OP;} if(tpo_parse_other) ++o;} while(++i<tpo_f);\
\
\
/* Set the incrementation. On a fully allocated iblock (when */ \
- /* `IN==iblock'), we have already gone through the whole */ \
+ /* `IN==tpo_iblock'), we have already gone through the whole */ \
/* array, so we'll set the incrementation to the size of the */ \
/* while block which will stop the `while' loop above. On a */ \
/* tile, we need to increment to the next contiguous patch */ \
/* of memory to continue parsing this tile. */ \
- increment += ( IN==iblock ? iblock->size \
- : gal_tile_block_increment(iblock, IN->dsize, \
- num_increment++, \
+ tpo_increment += ( IN==tpo_iblock ? tpo_iblock->size \
+ : gal_tile_block_increment(tpo_iblock, IN->dsize, \
+ tpo_num_increment++, \
NULL) ); \
} \
\
@@ -249,43 +251,43 @@ gal_tile_full_free_contents(struct
gal_tile_two_layer_params *tl);
}
-#define GAL_TILE_PO_OSET(OT, OP, IN, OUT, PARSE_OUT, CHECK_BLANK) { \
- switch(iblock->type) \
+#define GAL_TILE_PO_OSET(OT, IN, OTHER, PARSE_OTHER, CHECK_BLANK, OP) { \
+ switch(tpo_iblock->type) \
{ \
case GAL_TYPE_UINT8: \
- GAL_TILE_PO_OISET(uint8_t, OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OISET(uint8_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_INT8: \
- GAL_TILE_PO_OISET(int8_t, OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OISET(int8_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_UINT16: \
- GAL_TILE_PO_OISET(uint16_t, OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OISET(uint16_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_INT16: \
- GAL_TILE_PO_OISET(int16_t, OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OISET(int16_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_UINT32: \
- GAL_TILE_PO_OISET(uint32_t, OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OISET(uint32_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_INT32: \
- GAL_TILE_PO_OISET(int32_t, OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OISET(int32_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_UINT64: \
- GAL_TILE_PO_OISET(uint64_t, OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OISET(uint64_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_INT64: \
- GAL_TILE_PO_OISET(int64_t, OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OISET(int64_t, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_FLOAT32: \
- GAL_TILE_PO_OISET(float, OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OISET(float, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_FLOAT64: \
- GAL_TILE_PO_OISET(double, OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OISET(double, OT,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
default: \
- { \
- fprintf(stderr, "type code %d not recognized in " \
- "`GAL_TILE_PO_OSET'", iblock->type); \
+ { /* `error' function might not be available for the user. */ \
+ fprintf(stderr, "GAL_TILE_PO_OSET: type code %d not recognized",\
+ tpo_iblock->type); \
exit(EXIT_FAILURE); \
} \
} \
@@ -293,159 +295,75 @@ gal_tile_full_free_contents(struct
gal_tile_two_layer_params *tl);
/* Parse over a region of memory (can be an n-dimensional tile or a fully
- allocated block of memory) and do a certain operation. If `OUT' is not
- NULL, it can also save the output of the operation into the `gal_data_t'
- that it points to. Note that OUT must either have only one element (for
- the whole input) or have exactly the same number of elements as the
- input (one value for one pixel/element of the input). The input
- arguments are:
-
- `OP': Operator: this can be any number of C experssions. This macro
- is going to define an `IT *i' variable which will increment
- over each element of the input array/tile. It will also define
- an `OT *o' which you can use to access the output. If
- `PARSE_OUT' is non-zero, then `o' will also be incremented to
- the same index element but in the output array. You can use
- these along with any other variable you define before this
- macro to process the input and store it in the output. Note
- that `i' and `o' will be incremented once in here, so don't
- increment them in `OP'.
-
- `IN': Input `gal_data_t', can be a tile or an allocated block of
- memory.
-
- `OUT': Output `gal_data_t'. It can be NULL. In that case, `o' will
- be NULL and should not be used. If given,
-
- `PARSE_OUT': Parse the output along with the input. When this is
- non-zero, then the `o' pointer (described in `OP') will be
- incremented to cover the same region of the input in the
- output.
-
- `CHECK_BLANK': If it is non-zero, then the input will be checked for
- blank values.
-
- For `OP' you have access to some other variables besides `i' and `o':
-
- `i': Pointer to this element of input to parse.
-
- `o': Pointer to corresponding element in output.
-
- `b': blank value in input't type.
-
- See `lib/statistics.c' for some example applications of this function.
-
- -----------
- TIPS/TRICKS
- -----------
-
- You can use a given tile on a dataset that it was not initialized with
- (but has the same size). You can do that like this:
-
- void *tarray;
- gal_data_t *tblock;
-
- tarray=tile->array;
- tblock=tile->block;
- tile->array=gal_tile_block_relative_to_other(tile, THE_OTHER_DATASET);
- tile->block=THE_OTHER_DATASET; <<-- `tile->block' must corrected
- AFTER `tile->array'.
- GAL_TILE_PARSE_OPERATE(...)
-
- tile->array=tarray;
- tile->block=tblock;
-
- You can work on one other array while working on the one that tile
- actually points to. To do that, you can make a fake tile and pass that as
- the `OUT' argument to this macro. While the name is `OUT' it can be used
- in any manner you like. To do this, you can take the following steps. In
- this example, `tile' is the actual tile that you have.
-
- gal_data_t *faketile;
-
- // These can be done outside a loop.
- faketile=gal_data_alloc(NULL, OTHER_DATASET->type, 1, &dsize,
- NULL, 0, -1, NULL, NULL, NULL);
- free(faketile->array);
- free(faketile->dsize);
- faketile->block=OTHER_DATASET;
- faketile->ndim=OTHER_DATASET->ndim;
-
- // These can be done in the loop over tiles.
- faketile->size=tile->size;
- faketile->dsize=tile->dsize;
- faketile->array=gal_tile_block_relative_to_other(tile, OTHER_DATASET);
-
- // Do your processing....
- GAL_TILE_PARSE_OPERATE(PROCESSING, tile, faketile, 1, 1);
-
- // Clean up.
- bintile->array=NULL;
- bintile->dsize=NULL;
- gal_data_free(bintile);
-*/
-#define GAL_TILE_PARSE_OPERATE(OP, IN, OUT, PARSE_OUT, CHECK_BLANK) { \
- int parse_out=(OUT && PARSE_OUT); \
- gal_data_t *iblock = gal_tile_block(IN); \
- gal_data_t *oblock = OUT ? gal_tile_block(OUT) : NULL; \
+ allocated block of memory) and do a certain operation. If `OTHER' is not
+ NULL, this macro will also parse it at the same time . Note that OTHER
+ must either have only one element (for the whole input) or have exactly
+ the same number of elements as the input (one value for one
+ pixel/element of the input). See the documentation for more on this
+ macro and some examples. */
+#define GAL_TILE_PARSE_OPERATE(IN, OTHER, PARSE_OTHER, CHECK_BLANK, OP) { \
+ int tpo_parse_other=(OTHER && PARSE_OTHER); \
+ gal_data_t *tpo_iblock = gal_tile_block(IN); \
+ gal_data_t *tpo_oblock = OTHER ? gal_tile_block(OTHER) : NULL; \
\
/* A small sanity check. */ \
- if( parse_out && gal_data_dsize_is_different(iblock, oblock) ) \
+ if( tpo_parse_other \
+ && gal_data_dsize_is_different(tpo_iblock, tpo_oblock) ) \
{ \
/* The `error' function, is a GNU extension. */ \
- fprintf(stderr, "when `PARSE_OUT' is non-zero, the " \
- "allocated block size of the input and output of " \
- "`GAL_TILE_PARSE_OPERATE' must be equal, but they are " \
- "not: %zu and %zu elements respectively)", \
- iblock->size, oblock->size); \
+ fprintf(stderr, "GAL_TILE_PARSE_OPERATE: when `PARSE_OTHER' is "\
+ "non-zero, the allocated block size of `IN' and " \
+ "`OTHER' must be equal, but they are not: %zu and %zu " \
+ "elements respectively)", tpo_iblock->size, \
+ tpo_oblock->size); \
exit(EXIT_FAILURE); \
} \
\
- /* First set the OUTPUT type. */ \
- if(OUT) \
- switch(oblock->type) \
+ /* First set the OTHER type. */ \
+ if(OTHER) \
+ switch(tpo_oblock->type) \
{ \
case GAL_TYPE_UINT8: \
- GAL_TILE_PO_OSET(uint8_t, OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OSET(uint8_t, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_INT8: \
- GAL_TILE_PO_OSET(int8_t, OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OSET(int8_t, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_UINT16: \
- GAL_TILE_PO_OSET(uint16_t, OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OSET(uint16_t,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_INT16: \
- GAL_TILE_PO_OSET(int16_t, OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OSET(int16_t, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_UINT32: \
- GAL_TILE_PO_OSET(uint32_t, OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OSET(uint32_t,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_INT32: \
- GAL_TILE_PO_OSET(int32_t, OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OSET(int32_t, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_UINT64: \
- GAL_TILE_PO_OSET(uint64_t, OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OSET(uint64_t,IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_INT64: \
- GAL_TILE_PO_OSET(int64_t, OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OSET(int64_t, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_FLOAT32: \
- GAL_TILE_PO_OSET(float, OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OSET(float, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
case GAL_TYPE_FLOAT64: \
- GAL_TILE_PO_OSET(double, OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OSET(double, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
break; \
default: \
{ \
fprintf(stderr, "type code %d not recognized in " \
- "`GAL_TILE_PARSE_OPERATE'", oblock->type); \
+ "`GAL_TILE_PARSE_OPERATE'", tpo_oblock->type); \
exit(EXIT_FAILURE); \
} \
} \
else \
- /* When `OUT==NULL', its type is irrelevant, we'll just use */ \
+ /* When `OTHER==NULL', its type is irrelevant, we'll just use */ \
/*`int' as a place holder. */ \
- GAL_TILE_PO_OSET(int, OP,IN,OUT,PARSE_OUT,CHECK_BLANK); \
+ GAL_TILE_PO_OSET(int, IN,OTHER,PARSE_OTHER,CHECK_BLANK,OP); \
}
diff --git a/lib/statistics.c b/lib/statistics.c
index 8922fe4..8db333e 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -58,17 +58,20 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
gal_data_t *
gal_statistics_number(gal_data_t *input)
{
+ uint64_t counter=0;
size_t dsize=1;
gal_data_t *out=gal_data_alloc(NULL, GAL_TYPE_UINT64, 1, &dsize,
NULL, 1, -1, NULL, NULL, NULL);
/* If there is no blank values in the input, then the total number is
just the size. */
- if(gal_blank_present(input, 0))
- { GAL_TILE_PARSE_OPERATE({++(*o);}, input, out, 0, 1) }
+ if(gal_blank_present(input, 0)) /* `{}' necessary for `else'. */
+ { GAL_TILE_PARSE_OPERATE(input, NULL, 0, 1, {++counter;}); }
else
- *((uint64_t *)(out->array)) = input->size;
+ counter = input->size;
+ /* Write the value into memory. */
+ *((uint64_t *)(out->array)) = counter;
return out;
}
@@ -91,7 +94,7 @@ gal_statistics_minimum(gal_data_t *input)
gal_type_max(out->type, out->array);
/* Parse the full input. */
- GAL_TILE_PARSE_OPERATE({*o = *i < *o ? *i : *o; ++n;}, input, out, 0, 0);
+ GAL_TILE_PARSE_OPERATE(input, out, 0, 0, {*o = *i < *o ? *i : *o; ++n;});
/* If there were no usable elements, set the output to blank, then
return. */
@@ -116,7 +119,7 @@ gal_statistics_maximum(gal_data_t *input)
gal_type_min(out->type, out->array);
/* Parse the full input. */
- GAL_TILE_PARSE_OPERATE({*o = *i > *o ? *i : *o; ++n;}, input, out, 0, 0);
+ GAL_TILE_PARSE_OPERATE(input, out, 0, 0, {*o = *i > *o ? *i : *o; ++n;});
/* If there were no usable elements, set the output to blank, then
return. */
@@ -139,7 +142,7 @@ gal_statistics_sum(gal_data_t *input)
/* Parse the dataset. Note that in `gal_data_alloc' we set the `clear'
flag to 1, so it will be 0.0f. */
- GAL_TILE_PARSE_OPERATE({++n; *o += *i;}, input, out, 0, 1);
+ GAL_TILE_PARSE_OPERATE(input, out, 0, 1, {++n; *o += *i;});
/* If there were no usable elements, set the output to blank, then
return. */
@@ -162,7 +165,7 @@ gal_statistics_mean(gal_data_t *input)
/* Parse the dataset. Note that in `gal_data_alloc' we set the `clear'
flag to 1, so it will be 0.0f. */
- GAL_TILE_PARSE_OPERATE({++n; *o += *i;}, input, out, 0, 1);
+ GAL_TILE_PARSE_OPERATE(input, out, 0, 1, {++n; *o += *i;});
/* Above, we calculated the sum and number, so if there were any elements
in the dataset (`n!=0'), divide the sum by the number, otherwise, put
@@ -187,7 +190,7 @@ gal_statistics_std(gal_data_t *input)
NULL, 1, -1, NULL, NULL, NULL);
/* Parse the input. */
- GAL_TILE_PARSE_OPERATE({++n; s += *i; s2 += *i * *i;}, input, out, 0, 1);
+ GAL_TILE_PARSE_OPERATE(input, out, 0, 1, {++n; s += *i; s2 += *i * *i;});
/* Write the standard deviation into the output. */
*((double *)(out->array)) = n ? sqrt( (s2-s*s/n)/n ) : GAL_BLANK_FLOAT64;
@@ -210,7 +213,7 @@ gal_statistics_mean_std(gal_data_t *input)
NULL, 1, -1, NULL, NULL, NULL);
/* Parse the input. */
- GAL_TILE_PARSE_OPERATE({++n; s += *i; s2 += *i * *i;}, input, out, 0, 1);
+ GAL_TILE_PARSE_OPERATE(input, out, 0, 1, {++n; s += *i; s2 += *i * *i;});
/* Write in the output values and return. */
((double *)(out->array))[0] = n ? s/n : GAL_BLANK_FLOAT64;
diff --git a/lib/tile.c b/lib/tile.c
index b6f061f..24b4029 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -297,10 +297,10 @@ gal_tile_series_from_minmax(gal_data_t *block, size_t
*minmax, size_t number)
the block pointer and return the `dsize' element of lowest-level
structure. */
gal_data_t *
-gal_tile_block(gal_data_t *input)
+gal_tile_block(gal_data_t *tile)
{
- while(input->block!=NULL) input=input->block;
- return input;
+ while(tile->block!=NULL) tile=tile->block;
+ return tile;
}
@@ -308,9 +308,7 @@ gal_tile_block(gal_data_t *input)
/* Return the increment necessary to start at the next series of contiguous
- memory (fastest dimension) associated with a tile. See
- `gal_tile_block_check_tiles' as one example application of this
- function.
+ memory (fastest dimension) associated with a tile.
1D and 2D cases are simple and need no extra explanation, but the case
for higher dimensions can be alittle more complicated, So we will go
@@ -409,7 +407,7 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
dimensionality of this array is irrelevant. Note that unlike
`tiles', `tilevalues' must be an array.
- `tiles': This will be parsed as a linked list (using the `next'
+ `tilesll': This will be parsed as a linked list (using the `next'
element). Internally, it might be stored as an array, but this
function doesn't care! The position of the tile over its block will
be determined according to the `block' element and the pointer of
@@ -457,15 +455,16 @@ gal_tile_block_write_const_value(gal_data_t *tilevalues,
gal_data_t *tilesll,
}
}
- /* Go over the tiles and write the values in. Note Recall that `tofill'
- has the same type as `tilevalues'. So we are using memcopy. */
+ /* Go over the tiles and write the values in. Recall that `tofill' has
+ the same type as `tilevalues'. So we are using memcopy. */
tile_ind=0;
for(tile=tilesll; tile!=NULL; tile=tile->next)
{
/* Set the pointer to use as input. */
in=gal_data_ptr_increment(tilevalues->array, tile_ind++, type);;
- GAL_TILE_PARSE_OPERATE({memcpy(o, in, gal_type_sizeof(type));},
- tile, tofill, 1, 0);
+ GAL_TILE_PARSE_OPERATE( tile, tofill, 1, 0, {
+ memcpy(o, in, gal_type_sizeof(type));
+ } );
}
return tofill;
@@ -475,12 +474,9 @@ gal_tile_block_write_const_value(gal_data_t *tilevalues,
gal_data_t *tilesll,
-/* Make a copy of the memory block in integer type and fill it with the ID
- of each tile, the non-filled areas have blank values. Finally, save the
- final array into a FITS file, specified with `filename'. This is done
- mainly for inspecting the positioning of tiles. We are using a signed
- 32-bit type because this is the standard FITS standard type for
- integers. */
+/* Make a copy of the memory block and fill it with the index of each tile
+ in `tilesll' (counting from 0). The non-filled areas will have blank
+ values. The output dataset will have a type of `GAL_TYPE_INT32'. */
gal_data_t *
gal_tile_block_check_tiles(gal_data_t *tilesll)
{
@@ -523,6 +519,45 @@ gal_tile_block_relative_to_other(gal_data_t *tile,
gal_data_t *other)
+/* To use within `gal_tile_full_blank_flag'. */
+static void *
+tile_block_blank_flag(void *in_prm)
+{
+ struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
+ gal_data_t *tile_ll=(gal_data_t *)(tprm->params);
+
+ size_t i;
+ gal_data_t *tile;
+
+ /* Check all the tiles given to this thread. */
+ for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
+ {
+ tile=&tile_ll[ tprm->indexs[i] ];
+ gal_blank_present(tile, 1);
+ }
+
+ /* Wait for all the other threads to finish. */
+ if(tprm->b) pthread_barrier_wait(tprm->b);
+ return NULL;
+}
+
+
+
+
+
+/* Update the blank flag on the tiles within the list of input tiles. */
+void
+gal_tile_block_blank_flag(gal_data_t *tile_ll, size_t numthreads)
+{
+ /* Go over all the tiles and update their blank flag. */
+ gal_threads_spin_off(tile_block_blank_flag, tile_ll,
+ gal_list_data_number(tile_ll), numthreads);
+}
+
+
+
+
+
@@ -663,6 +698,11 @@ gal_tile_full_regular_first(gal_data_t *parent, size_t
*regular,
though they have different `block' datasets (that then link to one
allocated space). See the `gal_tile_full_two_layers' below.
+ `firsttsize': The size of the first tile along every dimension. This
+ is only different from the regular tile size when `regular' is not
+ an exact multiple of `input''s length along every dimension. This
+ array is allocated internally by this function.
+
Output
------
@@ -682,8 +722,7 @@ gal_tile_full_regular_first(gal_data_t *parent, size_t
*regular,
block in each dimension.
2. parent-size (or `psize') which is the size of the parent in each
- dimension (we don't want to go out of the paren't range).
-*/
+ dimension (we don't want to go out of the paren't range). */
size_t *
gal_tile_full(gal_data_t *input, size_t *regular,
float remainderfrac, gal_data_t **out, size_t multiple,
@@ -970,50 +1009,7 @@ gal_tile_full_two_layers(gal_data_t *input,
`permutation' element. If a permutation has already been defined for the
tessellation, this function will not do anythin. If permutation won't be
necessary, then this function will just return (the permutation must
- have been initialized to NULL).
-
-
- Theory
- ------
-
- When there is only one channel OR one dimension, the tiles are allocated
- in memory in the same order that they represent the input data. However,
- to make channel-independent processing possible in a generic way, the
- tiles of each channel are allocated contiguously. So, when there is more
- than one channel AND more than one dimension, the index of the tile does
- not correspond to its position in the grid covering the input dataset.
-
- The example below may help clarify: assume you have a 6x6 tessellation
- with two channels in the horizontal and one in the vertical. On the left
- you can see how the tile IDs correspond to the input dataset. NOTE how
- `03' is on the second row, not on the first after `02'. On the right,
- you can see how the tiles are stored in memory (and shown if you simply
- write the array into a FITS file for example).
-
- Corresponding to input In memory
- ---------------------- --------------
- 15 16 17 33 34 35 30 31 32 33 34 35
- 12 13 14 30 31 32 24 25 26 27 28 29
- 09 10 11 27 28 29 18 19 20 21 22 23
- 06 07 08 24 25 26 <-- 12 13 14 15 16 17
- 03 04 05 21 22 23 06 07 08 09 10 11
- 00 01 02 18 19 20 00 01 02 03 04 05
-
- As a result, if your values are stored in same order as the tiles, and
- you want them in over-all memory (for example to save as a FITS file),
- you need to permute the values:
-
- gal_permutation_apply(values, tl->permutation);
-
- If you have values over-all and you want them in tile-order, you can
- apply the inverse permutation:
-
- gal_permutation_apply_inverse(values, tl->permutation);
-
- Recall that this is the definition of permutation in this context:
-
- permute: IN_ALL[ i ] = IN_MEMORY[ perm[i] ]
- inverse: IN_ALL[ perm[i] ] = IN_MEMORY[ i ] */
+ have been initialized to NULL). */
void
gal_tile_full_permutation(struct gal_tile_two_layer_params *tl)
{
@@ -1217,45 +1213,6 @@ gal_tile_full_id_from_coord(struct
gal_tile_two_layer_params *tl,
-/* To use within `gal_tile_full_blank_bit'. */
-static void *
-tile_full_blank_flag(void *in_prm)
-{
- struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
- gal_data_t *tile_ll=(gal_data_t *)(tprm->params);
-
- size_t i;
- gal_data_t *tile;
-
- /* Check all the tiles given to this thread. */
- for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
- {
- tile=&tile_ll[ tprm->indexs[i] ];
- if( gal_blank_present(tile, 0) )
- tile->flag |= GAL_DATA_FLAG_HASBLANK;
- }
-
- /* Wait for all the other threads to finish. */
- if(tprm->b) pthread_barrier_wait(tprm->b);
- return NULL;
-}
-
-
-
-
-
-/* Update the blank flag on the tiles within the list of input tiles. */
-void
-gal_tile_full_blank_flag(gal_data_t *tile_ll, size_t numthreads)
-{
- /* Go over all the tiles and update their blank flag. */
- gal_threads_spin_off(tile_full_blank_flag, tile_ll,
- gal_list_data_number(tile_ll), numthreads);
-}
-
-
-
-
/* Clean up the allocated spaces in the parameters. */
void
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 75ba49e: Tile and Arithmetic libraries documented in manual,
Mohammad Akhlaghi <=