[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master df716bc 097/125: Tesselation with the new data
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master df716bc 097/125: Tesselation with the new data structure: work started |
|
Date: |
Sun, 23 Apr 2017 22:36:47 -0400 (EDT) |
branch: master
commit df716bcc8d1535b4e1d64054fcc017a7187fe7cb
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Tesselation with the new data structure: work started
Work has started on implementing the tessellation processing with the new
`gal_data_t'. Thanks to its starter pointer and sizes, it was just
necessary to define a new `block' element. From now on, when `block' is
non-NULL, it means that the array in this dataset is not actually allocated
but is only one part of it.
The old `mesh' and `spatialconvolve' library source files have been removed
and we now have two new source files called `tile' and `convolve'. The old
routines in the former source files are being transferred to the new source
files.
An intial implementation of defining a full, non-overlapping tessellation
has been done in `gal_tile_all_position'. Using this function, we have
defined `gal_tile_all_position_two_layers' which will define the channels
and tiles tessellations. Previously this was done in a (very complicated)
mesh structure (`gal_mesh_params') that was not very trivial to navigate
and understand. The old structure was mainly focused on efficiency and only
defined for regular non-overlapping grids/tessellations, but was hard to
follow and less general.
The new system will allow arbitrary tiles defined anywhere (even
overlapping) over an allocated array. It is also extremely easy to
understand, since there is no new structure, it is the same old
`gal_data_t' that a developer is already fully familiar with and is
commonly used in many libraries/programs.
This was mainly done in preparation of porting NoiseChisel and SubtractSky
to `gal_data_t', but currently it is being tested on Convolve's
spatial-domain convolution. Once it has been fully tested and implemented
here, it will be really easy to apply it in those programs. Infact,
SubtractSky will probably be removed, because we can define tessellation
feature in Statistics which can effectively estimate the Sky value (and
many other statistics) on a grid. Then it is possible to use Arithmetic to
actually subtract the derived image from the input image.
The tile-related options (currently in Convolve, but later also in
NoiseChisel and Statistics) have also been updated to be more easier. Also,
in the tessellation of the full dataset, the remainder area would go into
the last tile in each dimension. But now it is the first tile that is
affected by any remainder area along each dimension and when the remainder
is significant, it will be divided between the first and last tiles.
Note that Convolve's spatial convolution is still not fully implemented yet
and its test will fail with this commit. The work done in Gnuastro has been
significant until this stage (defining the tessellation). So, this commit
is being made before starting work on each tile.
Some other modifications with this commit:
- The functions that work with arrays of data structures have been
separated as `gal_data_array_*', currently there is only two: one for
clear-allocation and the other freeing the arrays and their contents.
- Since the `--type' option is now a common option to all programs, a new
`gal_fits_img_write_to_type' was defined to make the job easy in each
program.
- In many cases, we want all the flags for `gal_arithmetic', so a new
`GAL_ARITHMETIC_FLAGS_ALL' was defined to make the codes more clear.
- Since the programs are increasingly becoming dimension-independent, a
`multidim' source library file has been defined to help in making the
functions more clear to read.
---
bin/convolve/args.h | 95 +-
bin/convolve/astconvolve.conf | 10 +-
bin/convolve/convolve.c | 123 +-
bin/convolve/main.c | 2 +-
bin/convolve/main.h | 55 +-
bin/convolve/ui.c | 224 +--
bin/convolve/ui.h | 19 +-
bin/fits/fits.c | 2 +-
bin/statistics/statistics.c | 43 +-
bin/table/ui.c | 2 +-
doc/gnuastro.texi | 55 +-
lib/Makefile.am | 22 +-
bin/convolve/main.c => lib/convolve.c | 38 +-
lib/data.c | 127 +-
lib/fits.c | 24 +-
lib/gnuastro/arithmetic.h | 4 +-
bin/convolve/main.c => lib/gnuastro/convolve.h | 51 +-
lib/gnuastro/data.h | 155 +-
lib/gnuastro/fits.h | 5 +
lib/gnuastro/mesh.h | 201 ---
lib/gnuastro/multidim.h | 68 +
lib/gnuastro/spatialconvolve.h | 96 --
lib/gnuastro/tile.h | 81 +
lib/mesh.c | 1868 ------------------------
lib/multidim.c | 170 +++
lib/options.c | 88 +-
lib/options.h | 5 +
lib/spatialconvolve.c | 237 ---
lib/tile.c | 471 ++++++
lib/txt.c | 2 +-
30 files changed, 1464 insertions(+), 2879 deletions(-)
diff --git a/bin/convolve/args.h b/bin/convolve/args.h
index 5a1796d..cbda046 100644
--- a/bin/convolve/args.h
+++ b/bin/convolve/args.h
@@ -119,98 +119,61 @@ struct argp_option program_options[] =
- /* Mesh grid options. */
+ /* Tile grid options. */
{
0, 0, 0, 0,
- "Mesh grid (only for spatial domain):",
+ "Tile grid (only for spatial domain):",
ARGS_GROUP_MESH_GRID
},
{
- "meshsize",
- ARGS_OPTION_KEY_MESHSIZE,
- "INT",
- 0,
- "Size of each mesh (tile) in the grid.",
- ARGS_GROUP_MESH_GRID,
- &p->mp.meshsize,
- GAL_DATA_TYPE_SIZE_T,
- GAL_OPTIONS_RANGE_GT_0,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
- "nch1",
- ARGS_OPTION_KEY_NCH1,
- "INT",
+ "tile",
+ ARGS_OPTION_KEY_TILE,
+ "INT[,INT]",
0,
- "Number of channels along first FITS axis.",
+ "Size of tiles along each dim. (FITS order).",
ARGS_GROUP_MESH_GRID,
- &p->mp.nch1,
+ &p->tile,
GAL_DATA_TYPE_SIZE_T,
GAL_OPTIONS_RANGE_GT_0,
GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
+ GAL_OPTIONS_NOT_SET,
+ gal_options_parse_sizes_reverse
},
{
- "nch2",
- ARGS_OPTION_KEY_NCH2,
- "INT",
+ "numchannels",
+ ARGS_OPTION_KEY_NUMCHANNELS,
+ "INT[,..]",
0,
- "Number of channels along second FITS axis.",
+ "No. of channels along each dim. (FITS order).",
ARGS_GROUP_MESH_GRID,
- &p->mp.nch2,
- GAL_DATA_TYPE_SIZE_T,
- GAL_OPTIONS_RANGE_GT_0,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
- "lastmeshfrac",
- ARGS_OPTION_KEY_LASTMESHFRAC,
- "FLT",
- 0,
- "Fraction of last mesh area to add new.",
- ARGS_GROUP_MESH_GRID,
- &p->mp.lastmeshfrac,
- GAL_DATA_TYPE_FLOAT32,
- GAL_OPTIONS_RANGE_GE_0_LE_1,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
- "lastmeshfrac",
- ARGS_OPTION_KEY_LASTMESHFRAC,
- "FLT",
- 0,
- "Fraction of last mesh area to add new.",
- ARGS_GROUP_MESH_GRID,
- &p->mp.lastmeshfrac,
- GAL_DATA_TYPE_FLOAT32,
- GAL_OPTIONS_RANGE_GE_0_LE_1,
+ &p->numchannels,
+ GAL_DATA_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
+ GAL_OPTIONS_NOT_SET,
+ gal_options_parse_sizes_reverse
},
{
- "fullconvolution",
- ARGS_OPTION_KEY_FULLCONVOLUTION,
+ "convoverch",
+ ARGS_OPTION_KEY_CONVOVERCH,
0,
0,
- "Ignore channels in spatial convolution.",
+ "Convolve over channel borders.",
ARGS_GROUP_MESH_GRID,
- &p->mp.fullconvolution,
+ &p->convoverch,
GAL_OPTIONS_NO_ARG_TYPE,
GAL_OPTIONS_RANGE_0_OR_1,
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
},
{
- "checkmesh",
- ARGS_OPTION_KEY_CHECKMESH,
+ "checktiles",
+ ARGS_OPTION_KEY_CHECKTILES,
0,
0,
- "File created to view mesh structure",
+ "Tile IDs in an image, the size of input.",
ARGS_GROUP_MESH_GRID,
- &p->checkmesh,
+ &p->checktiles,
GAL_OPTIONS_NO_ARG_TYPE,
GAL_OPTIONS_RANGE_0_OR_1,
GAL_OPTIONS_NOT_MANDATORY,
@@ -236,13 +199,13 @@ struct argp_option program_options[] =
{
"makekernel",
ARGS_OPTION_KEY_MAKEKERNEL,
- 0,
+ "INT",
0,
"Make 2*INT kernel to create input image.",
GAL_OPTIONS_GROUP_OPERATING_MODE,
&p->makekernel,
- GAL_OPTIONS_NO_ARG_TYPE,
- GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_DATA_TYPE_SIZE_T,
+ GAL_OPTIONS_RANGE_GT_0,
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
},
diff --git a/bin/convolve/astconvolve.conf b/bin/convolve/astconvolve.conf
index 34f746f..908f523 100644
--- a/bin/convolve/astconvolve.conf
+++ b/bin/convolve/astconvolve.conf
@@ -22,15 +22,13 @@
khdu 1
# Output:
+ type float32
# Mesh grid:
- meshsize 32
- nch1 1
- nch2 1
- lastmeshfrac 0.6
- fullconvolution 0
+ tile 30,30
+ numchannels 1,1
+ convoverch 0
# Operating mode:
domain frequency
- makekernel 0
minsharpspec 0.005
\ No newline at end of file
diff --git a/bin/convolve/convolve.c b/bin/convolve/convolve.c
index 27bf71f..ba51617 100644
--- a/bin/convolve/convolve.c
+++ b/bin/convolve/convolve.c
@@ -30,6 +30,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gsl/gsl_errno.h>
#include <gnuastro/wcs.h>
+#include <gnuastro/tile.h>
#include <gnuastro/fits.h>
#include <gnuastro/threads.h>
#include <gnuastro/spatialconvolve.h>
@@ -188,21 +189,25 @@ void
makepaddedcomplex(struct convolveparams *p)
{
size_t i, ps0, ps1;
- double *o, *op, *pimg, *pker;
- float *f, *fp, *input=p->input, *kernel=p->kernel;
- size_t is0=p->is0, is1=p->is1, ks0=p->ks0, ks1=p->ks1;
+ float *f, *ff, *kernel=p->kernel->array;
+ size_t is0=p->input->dsize[0], is1=p->input->dsize[1];
+ size_t ks0=p->kernel->dsize[0], ks1=p->kernel->dsize[1];
+ double *o, *op, *d, *df, *pimg, *pker, *input=p->input->array;
+
/* Find the sizes of the padded image, note that since the kernel
sizes are always odd, the extra padding on the input image is
always going to be an even number (clearly divisable). */
- ps0=p->ps0 = p->makekernel ? p->is0 : p->is0+p->ks0-1;
- ps1=p->ps1 = p->makekernel ? p->is1 : p->is1+p->ks1-1;
+ ps0=p->ps0 = p->makekernel ? is0 : is0 + ks0 - 1;
+ ps1=p->ps1 = p->makekernel ? is1 : is1 + ks1 - 1;
+
/* The Discrete Fourier transforms operate faster on even-sized
arrays. So if the padded sides are not even, make them so: */
if(ps0%2) ps0=p->ps0=ps0+1;
if(ps1%2) ps1=p->ps1=ps1+1;
+
/* Allocate the space for the padded input image and fill it. */
errno=0;
pimg=p->pimg=malloc(2*ps0*ps1*sizeof *pimg);
@@ -214,8 +219,8 @@ makepaddedcomplex(struct convolveparams *p)
op=(o=pimg+i*2*ps1)+2*ps1; /* pimg is complex. */
if(i<is0)
{
- fp=(f=input+i*is1)+is1;
- do {*o++=*f; *o++=0.0f;} while(++f<fp);
+ df=(d=input+i*is1)+is1;
+ do {*o++=*d; *o++=0.0f;} while(++d<df);
}
do *o++=0.0f; while(o<op);
}
@@ -232,8 +237,8 @@ makepaddedcomplex(struct convolveparams *p)
op=(o=pker+i*2*ps1)+2*ps1; /* pker is complex. */
if(i<ks0)
{
- fp=(f=kernel+i*ks1)+ks1;
- do {*o++=*f; *o++=0.0f;} while(++f<fp);
+ ff=(f=kernel+i*ks1)+ks1;
+ do {*o++=*f; *o++=0.0f;} while(++f<ff);
}
do *o++=0.0f; while(o<op);
}
@@ -243,19 +248,19 @@ makepaddedcomplex(struct convolveparams *p)
-/* Remove the padding from the final convolved image and also correct
- for roundoff errors. Notice that we are putting the pixels in the
- input image of float type.
+/* Remove the padding from the final convolved image and also correct for
+ roundoff errors.
- NOTE: The padding to the input image (on the first axis for
- example) was p->ks0-1. Since p->ks0 is always odd, the
- padding will alwys be even. */
+ NOTE: The padding to the input image (on the first axis for example)
+ was `p->kernel->dsize[0]-1'. Since `p->kernel->dsize[0]' is
+ always odd, the padding will always be even. */
void
removepaddingcorrectroundoff(struct convolveparams *p)
{
size_t ps1=p->ps1;
- float *o, *input=p->input;
- size_t i, hi0, hi1, is0, is1;
+ size_t i, hi0, hi1;
+ size_t *isize=p->input->dsize;
+ double *o, *input=p->input->array;
double *d, *df, *start, *pimg=p->pimg;
/* Set all the necessary parameters to crop the desired region. hi0
@@ -265,26 +270,24 @@ removepaddingcorrectroundoff(struct convolveparams *p)
region that contains non-zero rows and columns.*/
if(p->makekernel)
{
- hi0 = 2*p->makekernel-1<p->is0 ? p->ps0/2-p->makekernel : 0;
- hi1 = 2*p->makekernel-1<p->is1 ? p->ps1/2-p->makekernel : 0;
- p->is0 = 2*p->makekernel-1<p->is0 ? 2*p->makekernel-1 : p->is0;
- p->is1 = 2*p->makekernel-1<p->is1 ? 2*p->makekernel-1 : p->is1;
+ hi0 = 2*p->makekernel-1 < isize[0] ? p->ps0/2-p->makekernel : 0;
+ hi1 = 2*p->makekernel-1 < isize[1] ? p->ps1/2-p->makekernel : 0;
+ isize[0] = 2*p->makekernel-1 < isize[0] ? 2*p->makekernel-1 : isize[0];
+ isize[1] = 2*p->makekernel-1 < isize[1] ? 2*p->makekernel-1 : isize[1];
}
else
{
- hi0=(p->ks0-1)/2;
- hi1=(p->ks1-1)/2;
+ hi0 = ( p->kernel->dsize[0] - 1 )/2;
+ hi1 = ( p->kernel->dsize[1] - 1 )/2;
}
- is0=p->is0;
- is1=p->is1;
/* To start with, `start' points to the first pixel in the final
image: */
start=&pimg[hi0*ps1+hi1];
- for(i=0;i<is0;++i)
+ for(i=0;i<isize[0];++i)
{
- o=&input[i*is1];
- df=(d=start+i*ps1)+is1;
+ o = &input[ i * isize[1] ];
+ df = ( d = start + i * ps1 ) + isize[1];
do
*o++ = ( *d<-CONVFLOATINGPOINTERR || *d>CONVFLOATINGPOINTERR )
? *d
@@ -765,60 +768,34 @@ frequencyconvolve(struct convolveparams *p)
void
convolve(struct convolveparams *p)
{
- float *convolved;
- long *meshindexs;
- gal_data_t *data;
- struct gal_mesh_params *mp=&p->mp;
- size_t ndim=2, dsize[]={p->is0, p->is1};
+ size_t ntiles, nch;
+ gal_data_t *tiles, *channels=NULL;
/* Do the convolution. */
if(p->domain==CONVOLVE_DOMAIN_SPATIAL)
{
- /* Prepare the mesh structure: */
- mp->img=p->input; mp->s0=p->is0; mp->s1=p->is1;
- mp->kernel=p->kernel; mp->ks0=p->ks0; mp->ks1=p->ks1;
- mp->numthreads=p->cp.numthreads;
- gal_mesh_make_mesh(mp);
- if(p->meshname)
- {
- /* Allocate the `gal_data_t' with the input imag. */
- gal_mesh_check_mesh_id(mp, &meshindexs);
- data=gal_data_alloc(p->mp.img, GAL_DATA_TYPE_FLOAT32, ndim, dsize,
- p->wcs, 0, p->cp.minmapsize, NULL, NULL, NULL);
- data->name="Input";
- gal_fits_img_write(data, p->meshname, NULL, PROGRAM_STRING);
-
- /* Change the array, type, and name. */
- data->array=meshindexs;
- data->type=GAL_DATA_TYPE_INT64;
- data->name="Mesh indexs";
- gal_fits_img_write(data, p->meshname, NULL, PROGRAM_STRING);
-
- /* Clean up. */
- data->name=NULL;
- free(meshindexs);
- gal_data_free(data);
- }
+ /* Prepare the mesh structure. */
+ gal_tile_all_position_two_layers(p->input, p->channel, p->tile,
+ &channels, &tiles, &ntiles, &nch);
+
+ /* Save it to the output if requested. */
+ if(p->tilesname)
+ error(EXIT_FAILURE, 0, "saving mesh structure to a file is not "
+ "yet implemented");
+
+ /* Do the spatial convolution. */
- /* Do the spatial convolution on the mesh: */
- gal_mesh_spatial_convolve_on_mesh(mp, &convolved);
/* Replace the input image array with the convolved array: */
- free(p->input);
- p->input=convolved;
+
+ /* Clean up */
+ gal_data_array_free(tiles, ntiles, 0);
+ gal_data_array_free(channels, nch, 0);
}
else
frequencyconvolve(p);
- /* Save the output (which is in p->input) array and clean up. Note that
- we will rely on `gal_data_free' to free some of the things that we
- don't need any more. */
- data=gal_data_alloc(p->input, GAL_DATA_TYPE_FLOAT32, ndim, dsize,
- NULL, 0, p->cp.minmapsize, NULL, NULL, NULL);
- data->wcs=p->wcs;
- data->unit=p->unit;
- data->name="Convolved";
- gal_fits_img_write(data, p->cp.output, NULL, PROGRAM_STRING);
- data->name=NULL;
- gal_data_free(data);
+ /* Save the output (which is in p->input) array. */
+ gal_fits_img_write_to_type(p->input, p->cp.output, NULL, PROGRAM_STRING,
+ p->cp.type);
}
diff --git a/bin/convolve/main.c b/bin/convolve/main.c
index 0a927bd..a4483d6 100644
--- a/bin/convolve/main.c
+++ b/bin/convolve/main.c
@@ -36,7 +36,7 @@ int
main(int argc, char *argv[])
{
struct timeval t1;
- struct convolveparams p={{0}, {0}, 0};
+ struct convolveparams p={{0}, 0};
/* Set the starting time.*/
time(&p.rawtime);
diff --git a/bin/convolve/main.h b/bin/convolve/main.h
index 5f3ccbd..55b785f 100644
--- a/bin/convolve/main.h
+++ b/bin/convolve/main.h
@@ -25,7 +25,6 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Include necessary headers */
#include <gnuastro/data.h>
-#include <gnuastro/mesh.h>
#include <options.h>
@@ -72,37 +71,33 @@ enum domain_codes
struct convolveparams
{
/* From command-line */
- struct gal_options_common_params cp; /* Common parameters. */
- struct gal_mesh_params mp; /* gal_mesh_params structure. */
- char *filename; /* Name of input file. */
- char *kernelname; /* File name of kernel. */
- char *khdu; /* HDU of kernel. */
- uint8_t nokernelflip; /* Do not flip the kernel. */
- uint8_t nokernelnorm; /* Do not normalize the kernel. */
- double minsharpspec; /* Deconvolution: min spect. of sharp img.*/
- uint8_t checkfreqsteps; /* View the frequency domain steps. */
- uint8_t checkmesh; /* View the mesh structure. */
- char *domainstr; /* String value specifying domain. */
- uint8_t makekernel; /* ==1: Make a kernel to create input. */
+ struct gal_options_common_params cp; /* Common parameters. */
+ char *filename; /* Name of input file. */
+ char *kernelname; /* File name of kernel. */
+ char *khdu; /* HDU of kernel. */
+ uint8_t nokernelflip; /* Do not flip the kernel. */
+ uint8_t nokernelnorm; /* Do not normalize the kernel. */
+ double minsharpspec; /* Deconvolution: min spect. of sharp img. */
+ uint8_t checkfreqsteps; /* View the frequency domain steps. */
+ size_t *tile; /* Size of tiles along each dim. (C order).*/
+ size_t *numchannels; /* No. of tiles along each dim. (C order). */
+ uint8_t convoverch; /* Convolve over channel borders. */
+ uint8_t checktiles; /* Tile IDs in an img, the size of input. */
+ char *domainstr; /* String value specifying domain. */
+ size_t makekernel; /* Make a kernel to create input. */
/* Internal */
- int domain; /* Frequency or spatial domain conv. */
- float *input; /* Input image array. */
- float *kernel; /* Input Kernel array. */
- time_t rawtime; /* Starting time of the program. */
- double *pimg; /* Padded image array. */
- double *pker; /* Padded kernel array. */
- size_t ps0; /* Padded size along first C axis. */
- size_t ps1; /* Padded size along second C axis. */
- size_t is0; /* Input image size along C's first axis. */
- size_t is1; /* Input image size along C's second axis.*/
- size_t ks0; /* Kernel size along C's first axis. */
- size_t ks1; /* Kernel size along C's second axis. */
- int nwcs; /* Number of WCS headers. */
- struct wcsprm *wcs; /* WCS structure. */
- char *unit; /* The unit string from the input. */
- char *freqstepsname; /* Name of file to check frequency steps. */
- char *meshname; /* Name of file to check mesh tiles. */
+ size_t *channel; /* Size of channels along each dimension. */
+ int domain; /* Frequency or spatial domain conv. */
+ gal_data_t *input; /* Input image array. */
+ gal_data_t *kernel; /* Input Kernel array. */
+ double *pimg; /* Padded image array. */
+ double *pker; /* Padded kernel array. */
+ size_t ps0; /* Padded size along first C axis. */
+ size_t ps1; /* Padded size along second C axis. */
+ char *freqstepsname; /* Name of file to check frequency steps. */
+ char *tilesname; /* Name of file to check tiles. */
+ time_t rawtime; /* Starting time of the program. */
};
#endif
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index 3f1d54e..ebaee5f 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -30,8 +30,12 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/wcs.h>
#include <gnuastro/fits.h>
+#include <gnuastro/tile.h>
#include <gnuastro/blank.h>
#include <gnuastro/table.h>
+#include <gnuastro/threads.h>
+#include <gnuastro/arithmetic.h>
+#include <gnuastro/statistics.h>
#include <gnuastro/linkedlist.h>
#include <timing.h>
@@ -118,23 +122,31 @@ ui_initialize_options(struct convolveparams *p,
/* Set the necessary common parameters structure. */
- cp->poptions = program_options;
+ cp->program_struct = p;
cp->program_name = PROGRAM_NAME;
cp->program_exec = PROGRAM_EXEC;
cp->program_bibtex = PROGRAM_BIBTEX;
cp->program_authors = PROGRAM_AUTHORS;
- cp->coptions = gal_commonopts_options;
+ cp->poptions = program_options;
cp->numthreads = gal_threads_number();
-
+ cp->coptions = gal_commonopts_options;
/* Set the mandatory common options. */
for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
switch(cp->coptions[i].key)
{
case GAL_OPTIONS_KEY_HDU:
+ case GAL_OPTIONS_KEY_TYPE:
case GAL_OPTIONS_KEY_MINMAPSIZE:
cp->coptions[i].mandatory=GAL_OPTIONS_MANDATORY;
break;
+
+ case GAL_OPTIONS_KEY_LOG:
+ case GAL_OPTIONS_KEY_SEARCHIN:
+ case GAL_OPTIONS_KEY_IGNORECASE:
+ case GAL_OPTIONS_KEY_TABLEFORMAT:
+ cp->coptions[i].flags=OPTION_HIDDEN;
+ break;
}
}
@@ -207,6 +219,41 @@ parse_opt(int key, char *arg, struct argp_state *state)
/*************** Sanity Check *******************/
/**************************************************************/
static void
+ui_read_check_only_options(struct convolveparams *p)
+{
+ /* Read the domain from a string into an integer. */
+ if( !strcmp("spatial", p->domainstr) )
+ p->domain=CONVOLVE_DOMAIN_SPATIAL;
+ else if( !strcmp("frequency", p->domainstr) )
+ p->domain=CONVOLVE_DOMAIN_FREQUENCY;
+ else
+ error(EXIT_FAILURE, 0, "domain value `%s' not recognized. Please use "
+ "either `spatial' or `frequency'", p->domainstr);
+
+ /* If we are in the spatial domain, make sure that the necessary
+ parameters are set. */
+ if( p->domain==CONVOLVE_DOMAIN_SPATIAL )
+ if( p->tile==NULL || p->numchannels==NULL )
+ {
+ if( p->tile==NULL && p->numchannels==NULL )
+ error(EXIT_FAILURE, 0, "in spatial convolution, `--numchannels' "
+ "and `--tile' are mandatory");
+ else
+ error(EXIT_FAILURE, 0, "in spatial convolution, `--%s' is "
+ "mandatory: you should use it to set the %s",
+ p->tile ? "numchannels" : "tile",
+ ( p->tile
+ ? "number of channels along each dimension of the input"
+ : "size of tiles to cover the input along each "
+ "dimension" ) );
+ }
+}
+
+
+
+
+
+static void
ui_check_options_and_arguments(struct convolveparams *p)
{
/* Make sure an input file name was given and if it was a FITS file, that
@@ -237,48 +284,25 @@ ui_check_options_and_arguments(struct convolveparams *p)
/**************************************************************/
/*************** Preparations *******************/
/**************************************************************/
-static void
-ui_read_domain(struct convolveparams *p)
-{
- if( !strcmp("spatial", p->domainstr) )
- p->domain=CONVOLVE_DOMAIN_SPATIAL;
- else if( !strcmp("frequency", p->domainstr) )
- p->domain=CONVOLVE_DOMAIN_FREQUENCY;
- else
- error(EXIT_FAILURE, 0, "domain value `%s' not recognized. Please use "
- "either `spatial' or `frequency'", p->domainstr);
-}
-
-
-
-
-
+/* Read the kernel. VERY IMPORTANT: We can't use the `fits_img_read_kernel'
+ because the Convolve program also does de-convolution. */
static void
ui_read_kernel(struct convolveparams *p)
{
float *f, *ff;
- gal_data_t *data;
/* Read the image into file. */
- data=gal_fits_img_read_to_type(p->kernelname, p->khdu,
- GAL_DATA_TYPE_FLOAT32, p->cp.minmapsize);
-
- /* Put its values into the main program structure. */
- p->ks0=data->dsize[0];
- p->ks1=data->dsize[1];
- p->kernel=data->array;
+ p->kernel = gal_fits_img_read_to_type(p->kernelname, p->khdu,
+ GAL_DATA_TYPE_FLOAT32,
+ p->cp.minmapsize);
/* Convert all the NaN pixels to zero if the kernel contains blank
pixels. */
- if(gal_blank_present(data))
+ if(gal_blank_present(p->kernel))
{
- ff=(f=data->array)+data->size;
+ ff = (f=p->kernel->array) + p->kernel->size;
do *f = isnan(*f) ? 0.0f : *f; while(++f<ff);
}
-
- /* Clean up. Note that we need the array, so it must be set to NULL. */
- data->array=NULL;
- gal_data_free(data);
}
@@ -288,17 +312,12 @@ ui_read_kernel(struct convolveparams *p)
static void
ui_preparations(struct convolveparams *p)
{
- double sum;
size_t i, size;
- gal_data_t *data;
- float *f, *ff, tmp;
+ gal_data_t *sum;
+ float *kernel, tmp;
char *outsuffix = p->makekernel ? "_kernel.fits" : "_convolved.fits";
- /* Read the domain */
- ui_read_domain(p);
-
-
/* Set the output name if the user hasn't set it. */
if(p->cp.output==NULL)
p->cp.output=gal_checkset_automatic_output(&p->cp, p->filename,
@@ -306,38 +325,37 @@ ui_preparations(struct convolveparams *p)
if(p->checkfreqsteps)
p->freqstepsname=gal_checkset_automatic_output(&p->cp, p->filename,
"_freqsteps.fits");
- if(p->checkmesh)
- p->meshname=gal_checkset_automatic_output(&p->cp, p->filename,
- "_mesh.fits");
+ if(p->checktiles)
+ p->tilesname=gal_checkset_automatic_output(&p->cp, p->filename,
+ "_tiled.fits");
- /* Read the input image as a float array and its WCS info. */
- data=gal_fits_img_read_to_type(p->filename, p->cp.hdu,
- GAL_DATA_TYPE_FLOAT32, p->cp.minmapsize);
- gal_wcs_read(p->filename, p->cp.hdu, 0, 0, &p->nwcs, &p->wcs);
- p->unit=data->unit;
- data->unit=NULL;
+ /* Read the input image as a float64 array and its WCS info. */
+ p->input=gal_fits_img_read_to_type(p->filename, p->cp.hdu,
+ GAL_DATA_TYPE_FLOAT64, p->cp.minmapsize);
+ gal_wcs_read(p->filename, p->cp.hdu, 0, 0, &p->input->nwcs, &p->input->wcs);
/* See if there are any blank values. */
- if(gal_blank_present(data) && p->domain==CONVOLVE_DOMAIN_FREQUENCY)
- fprintf(stderr, "\n----------------------------------------\n"
- "######## %s WARNING ########\n"
- "There are blank pixels in `%s' (hdu: `%s') and you have asked "
- "for frequency domain convolution. As a result, all the pixels "
- "in `%s' will be blank. Only spatial domain convolution can "
- "account for blank pixels in the input data. You can run %s "
- "again with `--domain=spatial'\n"
- "----------------------------------------\n\n",
- PROGRAM_NAME, p->filename, p->cp.hdu, p->cp.output, PROGRAM_NAME);
-
-
- /* Get all the information from `gal_data_t' then free it. */
- p->is0=data->dsize[0];
- p->is1=data->dsize[1];
- p->input=data->array;
- data->array=NULL;
- gal_data_free(data);
+ if(p->domain==CONVOLVE_DOMAIN_FREQUENCY)
+ {
+ if( gal_blank_present(p->input) )
+ fprintf(stderr, "\n----------------------------------------\n"
+ "######## %s WARNING ########\n"
+ "There are blank pixels in `%s' (hdu: `%s') and you have "
+ "asked for frequency domain convolution. As a result, all "
+ "the pixels in the output (`%s') will be blank. Only "
+ "spatial domain convolution can account for blank pixels "
+ "in the input data. You can run %s again with "
+ "`--domain=spatial'\n"
+ "----------------------------------------\n\n",
+ PROGRAM_NAME, p->filename, p->cp.hdu, p->cp.output,
+ PROGRAM_NAME);
+ }
+ else
+ p->channel=gal_tile_all_sanity_check(p->filename, p->cp.hdu, p->input,
+ p->tile, p->numchannels);
+
/* Read the file specified by --kernel. If makekernel is specified, then
@@ -345,25 +363,29 @@ ui_preparations(struct convolveparams *p)
argument) is the blurry image. */
if(p->makekernel)
{
- /* Read in the kernel array: */
+ /* Read in the kernel array. */
ui_read_kernel(p);
/* Make sure the size of the kernel is the same as the input */
- if(p->ks0!=p->is0 || p->ks1!=p->is1)
+ if( p->input->dsize[0]!=p->kernel->dsize[0]
+ || p->input->dsize[1]!=p->kernel->dsize[1] )
error(EXIT_FAILURE, 0, "with the `--makekernel' (`-m') option, "
"the input image and the image specified with the `--kernel' "
"(`-k') option should have the same size. The lower resolution "
"input image (%s) has %zux%zu pixels while the sharper image "
"(%s) specified with the kernel option has %zux%zu pixels",
- p->filename, p->is1, p->is0, p->kernelname, p->ks1, p->ks0);
-
- /* Divide both images by their sum so their lowest frequency
- becomes 1 (and their division would be meaningful!).*/
- size=p->is0*p->is1;
- sum=0.0f; ff=(f=p->input)+size; do sum+=*f++; while(f<ff);
- f=p->input; do *f++ /= sum; while(f<ff);
- sum=0.0f; ff=(f=p->kernel)+size; do sum+=*f++; while(f<ff);
- f=p->kernel; do *f++ /= sum; while(f<ff);
+ p->filename, p->input->dsize[1], p->input->dsize[0],
+ p->kernelname, p->kernel->dsize[1], p->kernel->dsize[0]);
+
+ /* Divide both images by their sum so their lowest frequency becomes
+ 1 and their division (in the frequency domain) would be
+ meaningful. */
+ sum=gal_statistics_sum(p->input);
+ p->input = gal_arithmetic(GAL_ARITHMETIC_OP_DIVIDE,
+ GAL_ARITHMETIC_FLAGS_ALL, p->input, sum);
+ sum=gal_statistics_sum(p->kernel);
+ p->kernel = gal_arithmetic(GAL_ARITHMETIC_OP_DIVIDE,
+ GAL_ARITHMETIC_FLAGS_ALL, p->kernel, sum);
}
/* Read the kernel. If there is anything particular to Convolve, then
@@ -377,41 +399,40 @@ ui_preparations(struct convolveparams *p)
{
/* Read in the kernel array: */
ui_read_kernel(p);
- size=p->ks0*p->ks1;
/* Check its size (must be odd). */
- if(p->ks0%2==0 || p->ks1%2==0)
+ if(p->kernel->dsize[0]%2==0 || p->kernel->dsize[1]%2==0)
error(EXIT_FAILURE, 0, "the kernel image has to have an odd "
"number of pixels on both sides (there has to be on pixel "
"in the center). %s (hdu: %s) is %zu by %zu",
- p->kernelname, p->khdu, p->ks1, p->ks0);
+ p->kernelname, p->khdu, p->kernel->dsize[1],
+ p->kernel->dsize[0]);
/* Normalize the kernel: */
if( !p->nokernelnorm )
{
- sum=0.0f; ff=(f=p->kernel)+size; do sum+=*f++; while(f<ff);
- f=p->kernel; do *f++ /= sum; while(f<ff);
+ sum=gal_statistics_sum(p->kernel);
+ p->kernel = gal_arithmetic(GAL_ARITHMETIC_OP_DIVIDE,
+ GAL_ARITHMETIC_FLAGS_ALL,
+ p->kernel, sum);
}
/* Flip the kernel: */
if( !p->nokernelflip )
- for(i=0;i<size/2;++i)
- {
- tmp = p->kernel[i];
- p->kernel[i] = p->kernel[size-i-1];
- p->kernel[size-i-1]=tmp;
- }
+ {
+ size=p->kernel->size;
+ kernel=p->kernel->array;
+ for(i=0;i<p->kernel->size/2;++i)
+ {
+ tmp = kernel[ i ];
+ kernel[ i ] = kernel[ size - i - 1 ];
+ kernel[ size - i - 1 ] = tmp;
+ }
+ }
}
else
- {
- data=gal_fits_img_read_kernel(p->kernelname, p->khdu,
- p->cp.minmapsize);
- p->ks0=data->dsize[0];
- p->ks1=data->dsize[1];
- p->kernel=data->array;
- data->array=NULL;
- gal_data_free(data);
- }
+ p->kernel = gal_fits_img_read_kernel(p->kernelname, p->khdu,
+ p->cp.minmapsize);
}
}
@@ -480,6 +501,10 @@ ui_read_check_inputs_setup(int argc, char *argv[], struct
convolveparams *p)
gal_options_read_config_set(&p->cp);
+ /* Do a sanity check only on options. */
+ ui_read_check_only_options(p);
+
+
/* Print the option values if asked. Note that this needs to be done
after the option checks so un-sane values are not printed in the
output state. */
@@ -528,9 +553,10 @@ ui_free_report(struct convolveparams *p, struct timeval
*t1)
{
/* Free the allocated arrays: */
free(p->khdu);
- free(p->kernel);
free(p->cp.hdu);
free(p->cp.output);
+ gal_data_free(p->input);
+ gal_data_free(p->kernel);
/* Print the final message. */
if(!p->cp.quiet)
diff --git a/bin/convolve/ui.h b/bin/convolve/ui.h
index f5f5616..85fe4c1 100644
--- a/bin/convolve/ui.h
+++ b/bin/convolve/ui.h
@@ -29,19 +29,18 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Available letters for short options:
- e f g i j l n r p t u v w x y z
- A B E F G H I J M O Q R T W X Y Z */
+ a b e f g i j l n p r s v w x y z
+ A B E F G J L M O Q R W X Y Z
+*/
enum option_keys_enum
{
/* With short-option version. */
ARGS_OPTION_KEY_KERNEL = 'k',
- ARGS_OPTION_KEY_KHDU = 'U',
- ARGS_OPTION_KEY_MINSHARPSPEC = 'c',
+ ARGS_OPTION_KEY_KHDU = 'u',
+ ARGS_OPTION_KEY_MINSHARPSPEC = 'H',
ARGS_OPTION_KEY_CHECKFREQSTEPS = 'C',
- ARGS_OPTION_KEY_MESHSIZE = 'c',
- ARGS_OPTION_KEY_NCH1 = 'a',
- ARGS_OPTION_KEY_NCH2 = 'b',
- ARGS_OPTION_KEY_LASTMESHFRAC = 'L',
+ ARGS_OPTION_KEY_TILE = 't',
+ ARGS_OPTION_KEY_NUMCHANNELS = 'c',
ARGS_OPTION_KEY_DOMAIN = 'd',
ARGS_OPTION_KEY_MAKEKERNEL = 'm',
@@ -49,8 +48,8 @@ enum option_keys_enum
automatically). */
ARGS_OPTION_KEY_NOKERNELFLIP = 1000,
ARGS_OPTION_KEY_NOKERNELNORM,
- ARGS_OPTION_KEY_CHECKMESH,
- ARGS_OPTION_KEY_FULLCONVOLUTION,
+ ARGS_OPTION_KEY_CHECKTILES,
+ ARGS_OPTION_KEY_CONVOVERCH,
};
diff --git a/bin/fits/fits.c b/bin/fits/fits.c
index bbce9d9..c319ef7 100644
--- a/bin/fits/fits.c
+++ b/bin/fits/fits.c
@@ -222,7 +222,7 @@ fits_print_extension_info(struct fitsparams *p)
/* Print the resutls. */
if(!p->cp.quiet)
{
- printf("%s\nRun on %s-----\n", PACKAGE_STRING, ctime(&p->rawtime));
+ printf("%s\nRun on %s-----\n", PROGRAM_STRING, ctime(&p->rawtime));
printf("HDU (extension) information: `%s'.\n", p->filename);
printf(" Column 1: Index (counting from 0).\n");
printf(" Column 2: Name (`EXTNAME' in FITS standard).\n");
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index 5bb9b58..8caa2a5 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -327,32 +327,29 @@ write_output_table(struct statisticsparams *p, gal_data_t
*table,
{
char *output;
int use_auto_output=0;
- char *fix, *suffix, *tmp;
+ char *fix, *suffix=NULL, *tmp;
struct gal_linkedlist_stll *comments=NULL;
- /* See if we should use automatic output. */
- if(p->cp.output)
- {
- if(p->numoutfiles>1) use_auto_output=1;
- }
- else use_auto_output=1;
+ /* Automatic output should be used when no output name was specified or
+ we have more than one output file. */
+ use_auto_output = p->cp.output ? (p->numoutfiles>1 ? 1 : 0) : 1;
- /* Set the type suffix */
+ /* Set the `fix' and `suffix' strings. Note that `fix' is necessary in
+ every case, even when no automatic output is to be used. Since it is
+ used to determine the format of the output. */
+ fix = ( p->cp.output
+ ? gal_fits_name_is_fits(p->cp.output) ? "fits" : "txt"
+ : "txt" );
if(use_auto_output)
- {
- fix = ( p->cp.output
- ? gal_fits_name_is_fits(p->cp.output) ? "fits" : "txt"
- : "txt" );
- asprintf(&suffix, "%s.%s", suf, fix);
- }
+ asprintf(&suffix, "%s.%s", suf, fix);
/* Make the output name. */
- output= ( use_auto_output
- ? gal_checkset_automatic_output(&p->cp, p->inputname, suffix)
- : p->cp.output );
+ output = ( use_auto_output
+ ? gal_checkset_automatic_output(&p->cp, p->inputname, suffix)
+ : p->cp.output );
/* Write the comments, NOTE: we are writing the first two in reverse of
@@ -379,7 +376,7 @@ write_output_table(struct statisticsparams *p, gal_data_t
*table,
/* Clean up. */
- free(suffix);
+ if(suffix) free(suffix);
if(output!=p->cp.output) free(output);
gal_linkedlist_free_stll(comments, 1);
}
@@ -621,13 +618,17 @@ print_basics(struct statisticsparams *p)
printf(" %-*s %.10g\n", namewidth, "Standard deviation:", std);
/* Ascii histogram. Note that we don't want to force the user to have the
- plotting parameters. */
- printf("-------\nHistogram:\n");
+ plotting parameters. Also, when a reference column is defined, the
+ range shown in the basic information section applies to that, not the
+ range of the histogram. In that case, we want to print the histogram
+ information. */
+ printf("-------");
p->asciiheight = p->asciiheight ? p->asciiheight : 10;
p->numasciibins = p->numasciibins ? p->numasciibins : 70;
bins=gal_statistics_regular_bins(p->input, NULL, p->numasciibins, NAN);
hist=gal_statistics_histogram(p->input, bins, 0, 0);
- print_ascii_plot(p, hist, bins, 1, 0);
+ if(p->refcol==NULL) printf("\nHistogram:\n");
+ print_ascii_plot(p, hist, bins, 1, p->refcol ? 1 : 0);
gal_data_free(bins);
gal_data_free(hist);
}
diff --git a/bin/table/ui.c b/bin/table/ui.c
index 6bd2ce0..383c563 100644
--- a/bin/table/ui.c
+++ b/bin/table/ui.c
@@ -278,7 +278,7 @@ ui_preparations(struct tableparams *p)
/* Print the file information. */
printf("--------\n");
tmp=gal_fits_name_save_as_string(p->filename, p->cp.hdu);
- printf("%s", tmp);
+ printf("%s\n", tmp);
free(tmp);
/* Print each column's information. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 53b0b43..5fdd9cd 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -5020,21 +5020,31 @@ should remove it.
@item Make
@cindex Make
address@hidden GNU Make
Make is a program for building ``targets'' (e.g., files) using ``recipes''
(a set of operations) when their known ``prerequisites'' (other files) have
been updated. It elegantly allows you to define dependency structures for
building your final output and updating it efficiently when the inputs
-change. It is most commonly used to build programs, but it is also very
-useful for scientific research, for example see
+change. It is the most common infra-structure to build software
+today.
+
+Scientific research methodology is very similar to software developement:
+you start by testing a hypothesis on a small sample of objects/targets with
+a simple set of steps. As you are able to get promising results, you
+improve the method and use it on a larger, more general, sample. In the
+process, you will confront many issues that have to be corrected (bugs in
+software development jargon). Make a wonderful tool to manage this style of
+development. It has been used to make reproducible papers, for example see
@url{https://gitlab.com/makhlaghi/NoiseChisel-paper, the reproduction
pipeline} of the paper introducing @ref{NoiseChisel} (one of Gnuastro's
-programs). GNU address@hidden@url{https://www.gnu.org/software/make/}} is
-the most common implementation which (like nearly all GNU programs comes
-with a wonderful
address@hidden@url{https://www.gnu.org/software/make/manual/}}). It is
-very basic and short (the most important parts are in the first roughly 100
-pages).
+programs).
+
address@hidden GNU Make
+GNU address@hidden@url{https://www.gnu.org/software/make/}} is the most
+common implementation which (similar to nearly all GNU programs, comes with
+a wonderful
address@hidden@url{https://www.gnu.org/software/make/manual/}}). Make is
+very basic and simple, and thus the manual is short (the most important
+parts are in the first roughly 100 pages) and easy to read/understand.
Make comes with a @option{--jobs} (@option{-j}) option which allows you to
specify the maximum number of jobs that can be done simultaneously. For
@@ -5049,15 +5059,16 @@ With this command, Make will process your
@file{Makefile} and create all
the targets (can be thousands of FITS images for example) simultaneously on
8 threads, while fully respecting their dependencies (only building a
file/target when its prerequisites are successfully built). Make is thus
-strongly recommended for managing scientific research where reproducibility
-and address@hidden its multi-threaded capabilities, Make will only
-re-build those targets that depend on a change you have made, not the whole
-work. For example, if you have set the prerequisites properly, you can
-easily test the changing of a parameter on your paper's results without
-having to re-do everything (which is much faster). This allows you to be
-much more productive in easily checking various ideas/assumptions of the
-different stages of your research and thus produce a more robust result for
-your exciting science.} are very important.
+strongly recommended for managing scientific research where robustness,
+archivability, reproducibility and address@hidden its
+multi-threaded capabilities, Make will only re-build those targets that
+depend on a change you have made, not the whole work. For example, if you
+have set the prerequisites properly, you can easily test the changing of a
+parameter on your paper's results without having to re-do everything (which
+is much faster). This allows you to be much more productive in easily
+checking various ideas/assumptions of the different stages of your research
+and thus produce a more robust result for your exciting science.} are very
+important.
@end table
@@ -16130,6 +16141,8 @@ and replace @code{PREFIX} with Gnuastro's installation
directory (see
@example
# Compile and correctly link Gnuastro related programs
+# example:
+# $ gnuastro-gcc myprog "-lfirst -lsecond"
function gnuastro-gcc @{
libtool --mode=link gcc $1.c $2 PREFIX/lib/libgnuastro.la -o $1
@}
@@ -16145,7 +16158,7 @@ you want to compile a C program that uses Gnuastro
library (for example
contains the source file:
@example
-gnuastro-gcc myprog
+$ gnuastro-gcc myprog
@end example
@noindent
@@ -16153,7 +16166,7 @@ If your program needs to link with libraries other than
those needed by
Gnuastro, you can specify them as a second argument (within double quotes):
@example
-gnuastro-gcc myprog "-lfirst -lsecond"
+$ gnuastro-gcc myprog "-lfirst -lsecond"
@end example
@noindent
@@ -16163,7 +16176,7 @@ example below. If your program can take arguments, then
you can also
include them afterwards.
@example
-gnuastro-gcc myprog && ./myprog
+$ gnuastro-gcc myprog && ./myprog
@end example
diff --git a/lib/Makefile.am b/lib/Makefile.am
index 5c40f9b..b3b4c73 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -47,10 +47,10 @@ libgnuastro_la_LIBADD =
$(top_builddir)/bootstrapped/lib/libgnu.la
# Specify the library .c files
-libgnuastro_la_SOURCES = arithmetic.c arithmetic-binary.c \
- arithmetic-onlyint.c blank.c box.c checkset.c data.c fits.c git.c \
- linkedlist.c mesh.c options.c polygon.c qsort.c spatialconvolve.c \
- statistics.c table.c threads.c timing.c txt.c wcs.c
+libgnuastro_la_SOURCES = arithmetic.c arithmetic-binary.c \
+ arithmetic-onlyint.c blank.c box.c checkset.c convolve.c data.c fits.c \
+ git.c linkedlist.c options.c polygon.c qsort.c multidim.c statistics.c \
+ table.c threads.c tile.c timing.c txt.c wcs.c
@@ -62,13 +62,13 @@ libgnuastro_la_SOURCES = arithmetic.c arithmetic-binary.c
\
# in the $(headersdir) directory. Some of the header files don't need to be
# installed.
headersdir=$(top_srcdir)/lib/gnuastro
-pkginclude_HEADERS = gnuastro/config.h $(headersdir)/arithmetic.h \
- $(headersdir)/blank.h $(headersdir)/box.h $(headersdir)/data.h \
- $(headersdir)/fits.h $(headersdir)/git.h $(headersdir)/linkedlist.h \
- $(headersdir)/mesh.h $(headersdir)/polygon.h $(headersdir)/qsort.h \
- $(headersdir)/spatialconvolve.h $(headersdir)/statistics.h \
- $(headersdir)/table.h $(headersdir)/threads.h $(headersdir)/wcs.h \
- $(headersdir)/txt.h
+pkginclude_HEADERS = gnuastro/config.h $(headersdir)/arithmetic.h \
+ $(headersdir)/blank.h $(headersdir)/box.h $(headersdir)/convolve.h \
+ $(headersdir)/data.h $(headersdir)/fits.h $(headersdir)/git.h
\
+ $(headersdir)/linkedlist.h $(headersdir)/multidim.h \
+ $(headersdir)/polygon.h $(headersdir)/qsort.h $(headersdir)/statistics.h \
+ $(headersdir)/table.h $(headersdir)/threads.h $(headersdir)/tile.h \
+ $(headersdir)/wcs.h $(headersdir)/txt.h
diff --git a/bin/convolve/main.c b/lib/convolve.c
similarity index 53%
copy from bin/convolve/main.c
copy to lib/convolve.c
index 0a927bd..ad75704 100644
--- a/bin/convolve/main.c
+++ b/lib/convolve.c
@@ -1,11 +1,11 @@
/*********************************************************************
-Convolve - Convolve input data with a given kernel.
-Convolve is part of GNU Astronomy Utilities (Gnuastro) package.
+convolve -- Convolve a dataset with a given kernel.
+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
@@ -23,34 +23,8 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <config.h>
#include <stdio.h>
+#include <errno.h>
+#include <error.h>
#include <stdlib.h>
-#include <timing.h> /* Includes time.h and sys/time.h */
-
-#include "main.h"
-
-#include "convolve.h"
-#include "ui.h" /* Needs convolveparams in main.h */
-
-int
-main(int argc, char *argv[])
-{
- struct timeval t1;
- struct convolveparams p={{0}, {0}, 0};
-
- /* Set the starting time.*/
- time(&p.rawtime);
- gettimeofday(&t1, NULL);
-
- /* Read the input parameters. */
- ui_read_check_inputs_setup(argc, argv, &p);
-
- /* Run Image Crop */
- convolve(&p);
-
- /* Free all non-freed allocations. */
- ui_free_report(&p, &t1);
-
- /* Return successfully. */
- return EXIT_SUCCESS;
-}
+#include <gnuastro/tile.h>
diff --git a/lib/data.c b/lib/data.c
index a73058e..8168ce4 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -608,6 +608,7 @@ gal_data_initialize(gal_data_t *data, void *array, uint8_t
type,
data->next=NULL;
data->ndim=ndim;
data->type=type;
+ data->block=NULL;
data->mmapname=NULL;
data->minmapsize=minmapsize;
gal_checkset_allocate_copy(name, &data->name);
@@ -709,46 +710,6 @@ gal_data_alloc(void *array, uint8_t type, size_t ndim,
size_t *dsize,
-/* Allocate an array of data structures and initialize all the values. */
-gal_data_t *
-gal_data_calloc_dataarray(size_t size)
-{
- size_t i;
- gal_data_t *out;
-
- /* Allocate the array to keep the structures. */
- errno=0;
- out=malloc(size*sizeof *out);
- if(out==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for `out' in "
- "`gal_data_calloc_dataarray'", size*sizeof *out);
-
-
- /* Set the pointers to NULL if they didn't exist and the non-pointers to
- impossible integers (so the caller knows the array is only
- allocated. `minmapsize' should be set when allocating the array and
- should be set when you run `gal_data_initialize'. */
- for(i=0;i<size;++i)
- {
- out[i].array = NULL;
- out[i].type = GAL_DATA_TYPE_INVALID;
- out[i].ndim = 0;
- out[i].dsize = NULL;
- out[i].nwcs = 0;
- out[i].wcs = NULL;
- out[i].mmapname = NULL;
- out[i].name = out[i].unit = out[i].comment = NULL;
- out[i].disp_fmt = out[i].disp_width = out[i].disp_precision = -1;
- }
-
- /* Return the array pointer. */
- return out;
-}
-
-
-
-
-
/* In some contexts, it is necessary for all the strings to have the same
allocated space (when the `strlen' is different). This function will
allocate new copies for all elements to have the same length as the
@@ -886,6 +847,92 @@ gal_data_free(gal_data_t *data)
/*********************************************************************/
+/************* Array of data structures ******************/
+/*********************************************************************/
+/* Allocate an array of data structures and initialize all the values. */
+gal_data_t *
+gal_data_array_calloc(size_t size)
+{
+ size_t i;
+ gal_data_t *out;
+
+ /* Allocate the array to keep the structures. */
+ errno=0;
+ out=malloc(size*sizeof *out);
+ if(out==NULL)
+ error(EXIT_FAILURE, errno, "%zu bytes for `out' in "
+ "`gal_data_calloc_dataarray'", size*sizeof *out);
+
+
+ /* Set the pointers to NULL if they didn't exist and the non-pointers to
+ impossible integers (so the caller knows the array is only
+ allocated. `minmapsize' should be set when allocating the array and
+ should be set when you run `gal_data_initialize'. */
+ for(i=0;i<size;++i)
+ {
+ out[i].array = NULL;
+ out[i].type = GAL_DATA_TYPE_INVALID;
+ out[i].ndim = 0;
+ out[i].dsize = NULL;
+ out[i].nwcs = 0;
+ out[i].wcs = NULL;
+ out[i].mmapname = NULL;
+ out[i].name = out[i].unit = out[i].comment = NULL;
+ out[i].disp_fmt = out[i].disp_width = out[i].disp_precision = -1;
+ }
+
+ /* Return the array pointer. */
+ return out;
+}
+
+
+
+
+
+/* When you have an array of data structures. */
+void
+gal_data_array_free(gal_data_t *data, size_t size, int free_array)
+{
+ size_t i;
+
+ /* If its NULL, don't do anything. */
+ if(data==NULL) return;
+
+ /* First free all the contents. */
+ for(i=0;i<size;++i)
+ {
+ /* See if the array should be freed or not. */
+ if(free_array==0)
+ data[i].array=NULL;
+
+ /* Now clear the contents of the dataset. */
+ gal_data_free_contents(&data[i]);
+ }
+
+ /* Now you can free the whole array. */
+ free(data);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*********************************************************************/
/************* Data structure as a linked list ******************/
/*********************************************************************/
/* Add a new data structure to the top of an existing linked list of data
diff --git a/lib/fits.c b/lib/fits.c
index caf50e8..161c6c9 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -1500,6 +1500,28 @@ gal_fits_img_write(gal_data_t *data, char *filename,
+void
+gal_fits_img_write_to_type(gal_data_t *data, char *filename,
+ struct gal_fits_key_ll *headers,
+ char *program_string, int type)
+{
+ /* If the input dataset is not the correct type, then convert it,
+ otherwise, use the input data structure. */
+ gal_data_t *towrite = (data->type==type
+ ? data
+ : gal_data_copy_to_new_type(data, type));
+
+ /* Write the converted dataset into an image. */
+ gal_fits_img_write(towrite, filename, headers, program_string);
+
+ /* Free the dataset if it was allocated. */
+ if(towrite!=data) gal_data_free(towrite);
+}
+
+
+
+
+
/* This function is mainly useful when you want to make FITS files in
parallel (from one main WCS structure, with just differing CRPIX) for
two reasons:
@@ -1794,7 +1816,7 @@ gal_fits_tab_info(char *filename, char *hdu, size_t
*numcols,
/* Read the total number of fields, then allocate space for the data
structure array and store the information within it. */
fits_read_key(fptr, TINT, "TFIELDS", &tfields, NULL, &status);
- allcols=gal_data_calloc_dataarray(tfields);
+ allcols=gal_data_array_calloc(tfields);
/* See comments of `fits_correct_bin_table_int_types'. Here we are
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index 9243f18..23b513f 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -65,7 +65,9 @@ __BEGIN_C_DECLS /* From C++ preparations */
#define GAL_ARITHMETIC_FREE 2
#define GAL_ARITHMETIC_NUMOK 4
-
+#define GAL_ARITHMETIC_FLAGS_ALL ( GAL_ARITHMETIC_INPLACE \
+ | GAL_ARITHMETIC_FREE \
+ | GAL_ARITHMETIC_NUMOK )
/* Identifiers for each operator. */
diff --git a/bin/convolve/main.c b/lib/gnuastro/convolve.h
similarity index 50%
copy from bin/convolve/main.c
copy to lib/gnuastro/convolve.h
index 0a927bd..1a18d53 100644
--- a/bin/convolve/main.c
+++ b/lib/gnuastro/convolve.h
@@ -1,11 +1,11 @@
/*********************************************************************
-Convolve - Convolve input data with a given kernel.
-Convolve is part of GNU Astronomy Utilities (Gnuastro) package.
+convolve -- Convolve the dataset with a given kernel.
+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
@@ -20,37 +20,34 @@ 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/>.
**********************************************************************/
-#include <config.h>
+#ifndef __GAL_TILE_H__
+#define __GAL_TILE_H__
-#include <stdio.h>
-#include <stdlib.h>
+/* Include other headers if necessary here. Note that other header files
+ must be included before the C++ preparations below */
+#include <gnuastro/data.h>
-#include <timing.h> /* Includes time.h and sys/time.h */
+/* 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 */
-#include "main.h"
+/* Actual header contants (the above were for the Pre-processor). */
+__BEGIN_C_DECLS /* From C++ preparations */
-#include "convolve.h"
-#include "ui.h" /* Needs convolveparams in main.h */
-int
-main(int argc, char *argv[])
-{
- struct timeval t1;
- struct convolveparams p={{0}, {0}, 0};
- /* Set the starting time.*/
- time(&p.rawtime);
- gettimeofday(&t1, NULL);
- /* Read the input parameters. */
- ui_read_check_inputs_setup(argc, argv, &p);
- /* Run Image Crop */
- convolve(&p);
- /* Free all non-freed allocations. */
- ui_free_report(&p, &t1);
- /* Return successfully. */
- return EXIT_SUCCESS;
-}
+__END_C_DECLS /* From C++ preparations */
+
+#endif
diff --git a/lib/gnuastro/data.h b/lib/gnuastro/data.h
index 31a614c..d1f4fa9 100644
--- a/lib/gnuastro/data.h
+++ b/lib/gnuastro/data.h
@@ -105,47 +105,126 @@ enum gal_data_types
/* Main data structure.
- Notes
- -----
-
- - If mmapname==NULL, then the array is allocated (using malloc, in the
- RAM), otherwise its is mmap'd (is actually a file on the ssd/hdd).
-
- - minmapsize is stored in the data structure to allow any derivative
- data structures to follow the same number and decide if they should
- be mmap'd or allocated.
-
- - `minmapsize' == 0: array is definitely mmap'd.
-
- - `minmapsize' == -1: array is definitely in RAM. */
+ mmap (keep data outside of RAM)
+ -------------------------------
+
+ `mmap' is C's facility to keep the data on the HDD/SSD instead of
+ inside the RAM. This can be very useful for large data sets which can
+ be very memory intensive. Ofcourse, the speed of operation greatly
+ decreases when defining not using RAM, but that is worth it because
+ increasing RAM might not be possible. So in `gal_data_t' when the size
+ of the requested array (in bytes) is larger than a certain minimum
+ size (in bytes), Gnuastro won't write the array in memory but on
+ non-volatile memory (like HDDs and SSDs) as a file in the
+ `./.gnuastro' directory of the directory it was run from.
+
+ - If mmapname==NULL, then the array is allocated (using malloc, in
+ the RAM), otherwise its is mmap'd (is actually a file on the
+ ssd/hdd).
+
+ - minmapsize is stored in the data structure to allow any
+ derivative data structures to follow the same number and decide
+ if they should be mmap'd or allocated.
+
+ - `minmapsize' == 0: array is definitely mmap'd.
+
+ - `minmapsize' == -1: array is definitely in RAM.
+
+
+ block (work with only a subset of the data)
+ -------------------------------------------
+
+ In many contexts, it is desirable to slice the data set into subsets
+ or tiles, not necessarily overlapping. In such a way that you can work
+ on each independently. One option is 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 `block' pointer in
+ `gal_data_t' is defined for 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 larger box is a
+ contiguous block of memory that you are interpretting as a 2D
+ array. But you only want to work on the region marked by the smaller
+ box.
+
+ tile->block = larger
+ ---------------------------
+ | |
+ | tile |
+ | ---------- |
+ | | | |
+ | |_ | |
+ | |*| | |
+ | ---------- |
+ |_ |
+ |*| |
+ ---------------------------
+
+ To use `gal_data_t's block concept, you allocate a `gal_data_t *tile'
+ which is initialized with the pointer to the first element in the
+ sub-array (as its `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 non-NULL pointer as `array', `gal_data_initialize' (and thus
+ `gal_data_alloc') do not allocate any space and just uses the given
+ pointer for the new `array' element of the `gal_data_t'. So your
+ `tile' data structure will not be pointing to a separately allocated
+ space.
+
+ After the allocation is done, you just point `tile->block' to the
+ `gal_data_t' which hosts the larger array. Afterwards, the programs
+ that take in the `sub' array can check the `block' pointer to see how
+ to deal with dimensions and increments (strides) while working on the
+ `sub' datastructure. For example to increment along the vertical
+ dimension, the program must look at index i+W (where `W' is the width
+ of the larger array, not the tile).
+
+ Since the block structure is defined as a pointer, arbitrary levels of
+ tesselation/griding are possible. Therefore, just like a linked list,
+ it is important to have the `block' pointer of the largest dataset set
+ to NULL. Normally, you won't have to worry about this, because
+ `gal_data_initialize' (and thus `gal_data_alloc') will set it to NULL
+ by default (just remember not to change it). You can then only change
+ the `block' element for the tiles you define over the allocated space.
+
+ In Gnuastro, there is a separate library for tiling operations called
+ `tile.h', see the functions there for tools to effectively use this
+ feature. This approach to dealing with parts of a larger block was
+ inspired from the way the GNU Scientific Library does it in the
+ "Vectors and Matrices" chapter.
+*/
typedef struct gal_data_t
{
/* Basic information on array of data. */
- void *array; /* Array keeping data elements. */
- uint8_t type; /* Type of data (from `gal_data_alltypes'). */
- size_t ndim; /* Number of dimensions in the array. */
- size_t *dsize; /* Size of array along each dimension. */
- size_t size; /* Total number of data-elements. */
- char *mmapname; /* File name of the mmap. */
- size_t minmapsize; /* Minimum number of bytes to mmap the array. */
+ void *array; /* Array keeping data elements. */
+ uint8_t type; /* Type of data (from `gal_data_alltypes'). */
+ size_t ndim; /* Number of dimensions in the array. */
+ size_t *dsize; /* Size of array along each dimension. */
+ size_t size; /* Total number of data-elements. */
+ char *mmapname; /* File name of the mmap. */
+ size_t minmapsize; /* Minimum number of bytes to mmap the array. */
/* WCS information. */
- int nwcs; /* for WCSLIB: no. coord. representations. */
- struct wcsprm *wcs; /* WCS information for this dataset. */
+ int nwcs; /* for WCSLIB: no. coord. representations. */
+ struct wcsprm *wcs; /* WCS information for this dataset. */
/* Content descriptions. */
- int status; /* Any context-specific status value. */
- char *name; /* e.g., EXTNAME, or column, or keyword. */
- char *unit; /* Units of the data. */
- char *comment; /* A more detailed description of the data. */
+ int status; /* Any context-specific status value. */
+ char *name; /* e.g., EXTNAME, or column, or keyword. */
+ char *unit; /* Units of the data. */
+ char *comment; /* A more detailed description of the data. */
/* For printing */
- int disp_fmt; /* See `gal_table_diplay_formats'. */
- int disp_width; /* Width of space to print in ASCII. */
- int disp_precision; /* Precision to print in ASCII. */
+ int disp_fmt; /* See `gal_table_diplay_formats'. */
+ int disp_width; /* Width of space to print in ASCII. */
+ int disp_precision; /* Precision to print in ASCII. */
- /* As linked list. */
- struct gal_data_t *next; /* To use it as a linked list if necessary. */
+ /* Pointers to other data structures. */
+ struct gal_data_t *next; /* To use it as a linked list if necessary. */
+ struct gal_data_t *block; /* `gal_data_t' of hosting block, see above. */
} gal_data_t;
@@ -202,9 +281,6 @@ gal_data_alloc(void *array, uint8_t type, size_t ndim,
size_t *dsize,
struct wcsprm *wcs, int clear, size_t minmapsize,
char *name, char *unit, char *comment);
-gal_data_t *
-gal_data_calloc_dataarray(size_t size);
-
size_t
gal_data_string_fixed_alloc_size(gal_data_t *data);
@@ -217,6 +293,17 @@ gal_data_free(gal_data_t *data);
/*********************************************************************/
+/************* Array of data structures ******************/
+/*********************************************************************/
+gal_data_t *
+gal_data_array_calloc(size_t size);
+
+void
+gal_data_array_free(gal_data_t *data, size_t num, int free_array);
+
+
+
+/*********************************************************************/
/************* Data structure as a linked list ******************/
/*********************************************************************/
void
diff --git a/lib/gnuastro/fits.h b/lib/gnuastro/fits.h
index c54e8e8..3a7010b 100644
--- a/lib/gnuastro/fits.h
+++ b/lib/gnuastro/fits.h
@@ -228,6 +228,11 @@ gal_fits_img_write(gal_data_t *data, char *filename,
struct gal_fits_key_ll *headers, char *program_string);
void
+gal_fits_img_write_to_type(gal_data_t *data, char *filename,
+ struct gal_fits_key_ll *headers,
+ char *program_string, int type);
+
+void
gal_fits_img_write_corr_wcs_str(gal_data_t *data, char *filename,
char *wcsheader, int nkeyrec, double *crpix,
struct gal_fits_key_ll *headers,
diff --git a/lib/gnuastro/mesh.h b/lib/gnuastro/mesh.h
deleted file mode 100644
index 66ac5ba..0000000
--- a/lib/gnuastro/mesh.h
+++ /dev/null
@@ -1,201 +0,0 @@
-/*********************************************************************
-meshgrid -- Create a mesh grid ready for multithreaded analysis.
-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.
-
-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_MESH_H__
-#define __GAL_MESH_H__
-
-/* Include other headers if necessary here. Note that other header files
- must be included before the C++ preparations below */
-#include <wcslib/wcs.h>
-#include <gnuastro/threads.h>
-
-
-
-/* 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 */
-
-
-
-/* Operations to do on each mesh. If input parameters are needed (for
- example the quantile), they are given as other argument(s) to the
- fillcharray function. */
-#define GAL_MESH_MAX_NUM_CHARRAY 2 /* Maximum number of charrays.
*/
-#define GAL_MESH_INTERP_ALL 1 /* Interpolate over the whole image as
one. */
-#define GAL_MESH_INTERP_CHANNEL 2 /* Interpolate over each channel
individually.*/
-
-
-/* The minimum number of acceptable nearest pixels. */
-#define GAL_MESH_MIN_ACCEPTABLE_NEAREST 3
-
-
-
-struct gal_mesh_thread_params
-{
- /* For convolve: */
- float *conv; /* The convolved array. */
- size_t *chbrd; /* Bordering x and y values all channels. */
-
- /* For all: */
- struct gal_mesh_params *mp; /* Pointer to gal_mesh_params structure. */
- size_t id; /* The thread ID starting from zero. */
-};
-
-
-
-
-struct gal_mesh_params
-{
- /* Image: */
- void *img; /* Input image array. */
- size_t s0; /* Height of input image. */
- size_t s1; /* Width of input image. */
-
- /* Threads: */
- size_t numthreads; /* Number of CPU threads. */
- size_t *indexs; /* 2D array of mesh indexs for each thread. */
- size_t thrdcols; /* The number of columns in indexs array. */
- pthread_barrier_t b; /* pthreads barrier for running threads. */
-
- /* Channel(s): */
- size_t nch; /* Total number of channels. */
- size_t nch1; /* Number of channels along first FITS axis. */
- size_t nch2; /* Number of channels along first FITS axis. */
- size_t gs0; /* Number of meshes on axis 0 in each channel. */
- size_t gs1; /* Number of meshes on axis 1 in each channel. */
-
- /* Meshs: */
- float lastmeshfrac; /* Last mesh fraction of size, to add new. */
- size_t meshsize; /* Size of each mesh. */
- size_t nmeshc; /* Number of meshes in each channel. */
- size_t nmeshi; /* Number of meshes in all image. */
- size_t *start; /* Starting pixel for each mesh. */
- size_t *types; /* The type of each mesh. */
- size_t *chindex; /* The index of each mesh in its channel. */
- size_t *imgindex; /* The index of each mesh in the image. */
- size_t maxs0; /* Maximum number of rows in all types. */
- size_t maxs1; /* Maximum number of columns in all types. */
-
- /* garrays: */
- int ngarrays; /* Number of garrays in this run. */
- float *garray1; /* Either equal to cgarray1 or fgarray1. */
- float *garray2; /* Either equal to cgarray2 or fgarray2. */
- float *cgarray1; /* In cgarray1 or cgarray2, the meshs in each */
- float *cgarray2; /* channel are contiguous. */
- float *fgarray1; /* In fgarray1 or fgarray2, we have contiguous */
- float *fgarray2; /* meshs in the full image. Ignore channels.*/
-
- /* Operate on each mesh: */
- void *params; /* Pointer to parameters structure of caller. */
- void *oneforall; /* One array that can contain all the meshs. */
-
- /* Interpolation: */
- float mirrordist; /* For finding the mode. Distance after mirror.*/
- float minmodeq; /* Minimum acceptable quantile for the mode. */
- unsigned char *byt; /* To keep track of pixels already checked. */
- size_t numnearest; /* Number of the nearest pixels for interp. */
- float *nearest1; /* Array keeping nearest pixels for garray1. */
- float *nearest2; /* Array keeping nearest pixels for garray2. */
- int interponlyblank; /* Only interpolate over blank pixels. */
- float *outgarray1; /* The interpolated garray1. */
- float *outgarray2; /* The interpolated garray2. */
- int fullinterpolation; /* ==1: Ignore channels in interpolation. */
- char *errstart; /* First sentence for error message. */
-
- /* Smoothing: */
- size_t smoothwidth; /* Width of smoothing kernel. */
- int fullsmooth; /* ==1: Ignore channels in smoothing. */
-
- /* Convolution: */
- float *kernel; /* Convolution kernel. */
- size_t ks0; /* Size of kernel along first C axis. */
- size_t ks1; /* Size of kernel along second C axis. */
- unsigned char fullconvolution; /* ==1: Convove over all channels. */
-
- /* Mesh types and information: */
- int meshbasedcheck; /* ==1: use one pixel for each mesh in checks. */
- size_t ts0[4]; /* Size (along first FITS axis) of mesh types. */
- size_t ts1[4]; /* Size (along second FITS axis) of mesh types.*/
-};
-
-size_t
-gal_mesh_ch_based_id_from_gid(struct gal_mesh_params *mp, size_t gid);
-
-size_t
-gal_mesh_gid_from_ch_based_id(struct gal_mesh_params *mp, size_t chbasedid);
-
-size_t
-gal_mesh_img_xy_to_mesh_id(struct gal_mesh_params *mp, size_t x, size_t y);
-
-void
-gal_mesh_check_mesh_id(struct gal_mesh_params *mp, long **out);
-
-void
-gal_mesh_check_garray(struct gal_mesh_params *mp, float **out1, float **out2);
-
-void
-gal_mesh_value_file(struct gal_mesh_params *mp, char *filename, char *extname1,
- char *extname2, struct wcsprm *wcs, char *spack_string);
-
-void
-gal_mesh_full_garray(struct gal_mesh_params *mp, int reverse);
-
-void
-gal_mesh_make_mesh(struct gal_mesh_params *mp);
-
-void
-gal_mesh_free_mesh(struct gal_mesh_params *mp);
-
-void
-gal_mesh_operate_on_mesh(struct gal_mesh_params *mp, void *(*meshfunc)(void *),
- size_t oneforallsize, int makegarray2, int
initialize);
-
-void
-gal_mesh_interpolate(struct gal_mesh_params *mp, char *errstart);
-
-void
-gal_mesh_smooth(struct gal_mesh_params *mp);
-
-void
-gal_mesh_spatial_convolve_on_mesh(struct gal_mesh_params *mp, float **conv);
-
-void
-gal_mesh_change_to_full_convolution(struct gal_mesh_params *mp, float *conv);
-
-
-
-__END_C_DECLS /* From C++ preparations */
-
-#endif /* __GAL_MESH_H__ */
diff --git a/lib/gnuastro/multidim.h b/lib/gnuastro/multidim.h
new file mode 100644
index 0000000..df51d4d
--- /dev/null
+++ b/lib/gnuastro/multidim.h
@@ -0,0 +1,68 @@
+/*********************************************************************
+multidim -- Functions for multi-dimensional operations.
+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_MULTIDIM_H__
+#define __GAL_MULTIDIM_H__
+
+/* Include other headers if necessary here. Note that other header files
+ must be included before the C++ preparations below */
+#include <gnuastro/data.h>
+
+/* 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 */
+
+
+/************************************************************************/
+/******************** Info **********************/
+/************************************************************************/
+size_t
+gal_multidim_total_size(size_t ndim, size_t *dsize);
+
+
+
+
+/************************************************************************/
+/******************** Coordinates **********************/
+/************************************************************************/
+size_t
+gal_multidim_coord_to_index(size_t ndim, size_t *dsize, size_t *coord);
+
+void
+gal_multidim_index_to_coord(size_t ind, size_t ndim, size_t *dsize,
+ size_t *coord);
+
+
+__END_C_DECLS /* From C++ preparations */
+
+#endif
diff --git a/lib/gnuastro/spatialconvolve.h b/lib/gnuastro/spatialconvolve.h
deleted file mode 100644
index c249b5c..0000000
--- a/lib/gnuastro/spatialconvolve.h
+++ /dev/null
@@ -1,96 +0,0 @@
-/*********************************************************************
-SpatialConvolve - Convolve an image in the spatial domain.
-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.
-
-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_SPATIALCONVOLVE_H__
-#define __GAL_SPATIALCONVOLVE_H__
-
-/* Include other headers if necessary here. Note that other header files
- must be included before the C++ preparations below */
-#include <gnuastro/threads.h> /* For pthread_barrier_t: */
-
-
-
-/* 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 */
-
-
-
-/* Main structure: */
-struct gal_spatialconvolve_params
-{
- /* General input parameters: */
- float *input; /* Input image array. */
- float *kernel; /* Kernel array. */
- float *out; /* Output image. */
- size_t is0; /* Image size along first C axis. */
- size_t is1; /* Image size along second C axis. */
- size_t ks0; /* Kernel size along first C axis. */
- size_t ks1; /* Kernel size along second C axis. */
- int edgecorrection; /* Correct the edges of the image. */
- long fpixel_i[2]; /* First pixel in input image. */
- long lpixel_i[2]; /* Last pixel in input image. */
- long fpixel_o[2]; /* First pixel in kernel. */
- long lpixel_o[2]; /* Last pixel in kernel. */
-
- /* Thread parameters. */
- size_t numthreads; /* Number of threads. */
- size_t *indexs; /* Indexs to be used in this thread. */
- pthread_barrier_t *b; /* Barrier to keep threads waiting. */
-};
-
-
-
-/* Functions: */
-void
-gal_spatialconvolve_pparams(float *input, size_t is0, size_t is1,
- float *kernel, size_t ks0, size_t ks1,
- size_t nt, int edgecorrection, float *out,
- size_t *indexs,
- struct gal_spatialconvolve_params *scp);
-
-void *
-gal_spatialconvolve_thread(void *inparam);
-
-void
-gal_spatialconvolve_convolve(float *input, size_t is0, size_t is1,
- float *kernel, size_t ks0, size_t ks1,
- size_t nt, int edgecorrection, float **out);
-
-
-
-__END_C_DECLS /* From C++ preparations */
-
-#endif /* __GAL_SPATIALCONVOLVE_H__ */
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
new file mode 100644
index 0000000..b53d683
--- /dev/null
+++ b/lib/gnuastro/tile.h
@@ -0,0 +1,81 @@
+/*********************************************************************
+tile -- work with tesselations over a host dataset.
+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_TILE_H__
+#define __GAL_TILE_H__
+
+/* Include other headers if necessary here. Note that other header files
+ must be included before the C++ preparations below */
+#include <gnuastro/data.h>
+
+/* 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 */
+
+
+
+
+
+/***********************************************************************/
+/************** About block **********************/
+/***********************************************************************/
+size_t *
+gal_tile_block_size(gal_data_t *input);
+
+void
+gal_tile_block_tile_start_coord(gal_data_t *tile, size_t *start_coord);
+
+
+/***********************************************************************/
+/************** Tile full dataset ********************/
+/***********************************************************************/
+size_t *
+gal_tile_all_sanity_check(char *filename, char *hdu, gal_data_t *input,
+ size_t *tile, size_t *numchannels);
+
+size_t
+gal_tile_all_position(gal_data_t *input, size_t *regular, gal_data_t **out,
+ size_t multiple);
+
+void
+gal_tile_all_position_two_layers(gal_data_t *input, size_t *channel_size,
+ size_t *tile_size, gal_data_t **channels,
+ gal_data_t **tiles, size_t *numchannels,
+ size_t *numtiles);
+
+
+
+
+__END_C_DECLS /* From C++ preparations */
+
+#endif
diff --git a/lib/mesh.c b/lib/mesh.c
deleted file mode 100644
index 0b81c7c..0000000
--- a/lib/mesh.c
+++ /dev/null
@@ -1,1868 +0,0 @@
-/*********************************************************************
-meshgrid -- Create a mesh grid ready for multithreaded analysis.
-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.
-
-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/>.
-**********************************************************************/
-#include <config.h>
-
-#include <math.h>
-#include <stdio.h>
-#include <errno.h>
-#include <error.h>
-#include <float.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include <gnuastro/fits.h>
-#include <gnuastro/mesh.h>
-#include <gnuastro/qsort.h>
-#include <gnuastro/linkedlist.h>
-#include <gnuastro/statistics.h>
-#include <gnuastro/spatialconvolve.h>
-
-#include "neighbors.h"
-
-
-
-/*
- Basic idea
- ==========
-
- The meshes in an image are completely independent. There is no
- overlap. So when we want to do any sort of processing on separate
- meshes, they can be done completely independently. Whether they are
- positioned in the same channel or not. Usually, after the initial
- processing is complete, and some meshes pass the test and some
- don't, we need to do interpolation to find values for all the
- meshes and finally smooth the result.
-
- The programs outside of mesh.c, should only be using garray (1 or
- 2). Not cgarray or fgarray. The reason for this is that the user
- might choose to ignore channels on one step and use them in
- another. All the functions in mesh.c will set garray and use
- cgarray and fgarray such that when they return, garray points to
- their output. So the caller to any of these functions doesn't have
- to worry about which one was last used, they can just ignorantly
- use garray and everything will be fine.
-
- Maybe in the future all these arrays can be removed and calculated
- immediately for each mesh during the processing. This might
- increase the speed and decrease the amount of necessary
- memory. However, it will add to the complexity of the code. So at
- this early stage, I am doing this job separately here for now so
- the main functions can be cleaner and more easily
- understandable. Once they have become mature enough, all these
- calculations should be done during the processing of each mesh.
-*/
-
-
-
-
-
-
-
-
-
-/*********************************************************************/
-/************ Finding the proper mesh ID ******************/
-/*********************************************************************/
-/* This is useful for when you are going over the elements in garray (and
-you are completley ignorant to which one of cgarrays or fgarrays garray
-points to) and you need the channel based IDs to get basic mesh information
-like the mesh type and size. */
-size_t
-gal_mesh_ch_based_id_from_gid(struct gal_mesh_params *mp, size_t gid)
-{
- if(mp->nch==1 || mp->garray1==mp->cgarray1)
- return gid;
- else
- {
- /* The X and Y of this mesh in the full image: */
- size_t f0=gid/(mp->gs1*mp->nch1), f1=gid%(mp->gs1*mp->nch1);
-
- /* The channel ID: */
- size_t chid = (f0/mp->gs0) * mp->nch1 + f1/mp->gs1;
-
- /* The ID of this mesh in this channel: */
- size_t inchannelid = (f0%mp->gs0) * mp->gs1 + f1%mp->gs1;
-
- /* For a check:
- printf("%zu:\n\t(f0, f1): (%zu, %zu)\t chid: %zu\tinchannelid: %zu\n",
- gid, f0, f1, chid, inchannelid);
- */
-
- /* Return the channel-based ID: */
- return chid * mp->nmeshc + inchannelid;
- }
-
- /* This function should not reach here! If it does there is problem,
- so just return an impossible value so the root of the issue can
- be found easily: */
- return (size_t)(-1);
-}
-
-
-
-
-
-/* Get the garray-based ID from the channel-based ID. See the comments above
- gal_mesh_ch_based_id_from_gid for a complete explanation. */
-size_t
-gal_mesh_gid_from_ch_based_id(struct gal_mesh_params *mp, size_t chbasedid)
-{
- if(mp->nch==1 || mp->garray1==mp->cgarray1)
- return chbasedid;
- else
- {
- /* The X and Y positions of this channel in the channels array: */
- size_t chx=(chbasedid/mp->nmeshc)/mp->nch1;
- size_t chy=(chbasedid/mp->nmeshc)%mp->nch1;
-
- /* The X and Y of this mesh in this channel: */
- size_t mx=(chbasedid%mp->nmeshc)/mp->gs1;
- size_t my=(chbasedid%mp->nmeshc)%mp->gs1;
-
- /* For a check:
- printf("%zu:\n\t(chx, chy): (%zu, %zu)"
- "\n\t(mx, my): (%zu, %zu)"
- "\n\t%zu\n\n",
- chbasedid, chx, chy, mx, my,
- (chx*mp->gs0+mx) * mp->nch1 + (chy*mp->gs1+my));
- */
-
- /* Return the : */
- return (chx*mp->gs0+mx) * mp->nch1 + (chy*mp->gs1+my);
- }
-
- /* This function should not reach here! If it does there is problem,
- so just return an impossible value so the root of the issue can
- be found easily: */
- return (size_t)(-1);
-}
-
-
-
-
-
-/* The user has a pixel index in the final image and wants to know the
- id it should plug into the garrays to get a value for this
- pixel. This is the job of this function. So to find the value on
- the mesh grid for a pixel at index `ind', the user should just run:
-
- mp->garray[imgindextomeshid(mp, ind)]
- */
-size_t
-gal_mesh_img_xy_to_mesh_id(struct gal_mesh_params *mp, size_t x, size_t y)
-{
- /* Take the proper action. The ternary conditional is here because
- when the meshsize is not an exact multiple of the the channel
- (image) size, there might be some extra pixels in the last mesh
- in each dimension which will cause trouble in the end. So without
- these checks, a pixel lying in those extra regions will be
- thought of as belongin to another mesh (that doesn't exist). We
- have to make sure that doesn't happen. */
- if(mp->nch==1)
- return ( (x/mp->meshsize<mp->gs0 ? x/mp->meshsize : x/mp->meshsize -1)
- * mp->gs1
- + (y/mp->meshsize<mp->gs1 ? y/mp->meshsize : y/mp->meshsize -1));
- else
- {
- /* Number of pixels along each axis in all channels: */
- size_t cps0=mp->s0/mp->nch2, cps1=mp->s1/mp->nch1;
-
- /* The X and Y positions of this channel in the channels array: */
- size_t chx=x/cps0, chy=y/cps1;
-
- /* The X and Y of this mesh in this channel: */
- size_t mx = (x%cps0)/mp->meshsize;
- size_t my = (y%cps1)/mp->meshsize;
-
- /* If the last mesh doesn't have the same size as the rest, mx
- or my might become one larger. Note that we have already made
- sure that this pixel is in the channel specified by chx and
- chy. */
- mx = mx<mp->gs0 ? mx : mx-1;
- my = my<mp->gs1 ? my : my-1;
-
-
- /* Return the proper id to input into garray. */
- if(mp->garray1==mp->cgarray1)
- return mp->nmeshc * (chx*mp->nch1+chy) + mx * mp->gs1 + my;
- else
- return (chx*mp->gs0+mx) * mp->gs1 + (chy * mp->gs1 + my);
- }
-
- /* This function should not reach here! So we will just return a
- value that will always be problematic ;-). */
- return (size_t) -1;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/*********************************************************************/
-/******************** Full garray ********************/
-/*********************************************************************/
-/* By default, the garray1 and garray2 arrays keep the meshes of each
- channel contiguously. So in practice, each channel is like a small
- independent image. This will cause problems when we want to work on the
- meshs irrespective of which channel they belong to. This function
- allocates and fills in the fgarray1 and fgarray2 arrays.
-
- The explanation above is for the case when reverse==0. If it is set
- equal to 1 (or any non-zero number), then
-*/
-void
-gal_mesh_full_garray(struct gal_mesh_params *mp, int reverse)
-{
- size_t nch1=mp->nch1;
- size_t ind, gs1=mp->gs1, gs0=mp->gs0;
- float *fgarray1=NULL, *fgarray2=NULL;
- float *cgarray1=mp->cgarray1, *cgarray2=mp->cgarray2;
- size_t g0, g1, f0, f1, fmeshind, chid, is1=mp->nch1*mp->gs1;
-
- /* First allocate the fgarrays if they were not already
- allocated. */
- if(mp->fgarray1==NULL)
- {
- /* A simple sanity check */
- if(reverse)
- error(EXIT_FAILURE, 0, "a bug! Please contact us at %s so we can "
- "fix this problem. For some reason, gal_mesh_full_garray "
- "has been called with the `reverse' flag set to true while "
- "fgarray is not allocated! This should not happen",
- PACKAGE_BUGREPORT);
-
- /* Allocate the fgarrays */
- errno=0; mp->fgarray1=malloc(mp->nmeshi*sizeof *mp->fgarray1);
- if(mp->fgarray1==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for mp->fgarray1 (mesh.c)",
- mp->nmeshi*sizeof *mp->fgarray1);
- }
- if(mp->ngarrays==2 && mp->fgarray2==NULL)
- {
- errno=0;
- mp->fgarray2=malloc(mp->nmeshi*sizeof *mp->fgarray2);
- if(mp->fgarray2==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for mp->fgarray2 (mesh.c)",
- mp->nmeshi*sizeof *mp->fgarray2);
- }
- fgarray1=mp->fgarray1;
- fgarray2=mp->fgarray2;
-
- /* Fill the fgarrays or cgarrays with the proper values: */
- for(chid=0;chid<mp->nch;++chid)
- {
- /* The first pixel ID in this channel is: chid*mp->nmeshi. */
- f0=(chid/nch1)*gs0; /* Position of first channel mesh */
- f1=(chid%nch1)*gs1; /* in full image. */
- fmeshind=chid*mp->nmeshc;
-
- /* Go over all the meshes in this channel and put them in their
- righful place. */
- for(ind=0;ind<mp->nmeshc;++ind)
- {
- g0=ind/gs1; /* Position of this mesh in this channel. */
- g1=ind%gs1;
- if(reverse)
- cgarray1[ind+fmeshind]=fgarray1[(g0+f0)*is1+g1+f1];
- else
- fgarray1[(g0+f0)*is1+g1+f1]=cgarray1[ind+fmeshind];
- if(mp->ngarrays==2)
- {
- if(reverse)
- cgarray2[ind+fmeshind]=fgarray2[(g0+f0)*is1+g1+f1];
- else
- fgarray2[(g0+f0)*is1+g1+f1]=cgarray2[ind+fmeshind];
- }
- }
- }
-
- /* Just for a check:
- gal_fits_array_to_file("nochannels.fits", "fgarray1", FLOAT_IMG,
- fgarray1, mp->nch2*mp->gs0, mp->nch1*mp->gs1, 1,
- NULL, NULL, "mesh");
- if(mp->ngarrays==2)
- gal_fits_array_to_file("nochannels.fits", "fgarray2", FLOAT_IMG,
- fgarray2, mp->nch2*mp->gs0, mp->nch1*mp->gs1,
- 1, NULL, NULL, "mesh");
- */
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/*********************************************************************/
-/******************** Checking functions ********************/
-/*********************************************************************/
-/* Save the meshid of each pixel into an array the size of the image. */
-void
-gal_mesh_check_mesh_id(struct gal_mesh_params *mp, long **out)
-{
- long i, *l, *lp;
- size_t row, start, *types=mp->types;
- size_t s0, s1, is1=mp->s1, *ts0=mp->ts0, *ts1=mp->ts1;
-
- /* Allocate the array to keep the mesh indexs. Calloc is used so we
- can add all the indexes to the existing value to make sure that
- there is no overlap. */
- errno=0;
- *out=calloc(mp->s0*mp->s1, sizeof **out);
- if(*out==NULL)
- error(EXIT_FAILURE, errno,
- "the array to show mesh labels in checkmesh");
-
- /* Fill the indexs: */
- for(i=0;i<mp->nmeshi;++i)
- {
- row=0;
- s0=ts0[types[i]];
- s1=ts1[types[i]];
- start=mp->start[i];
- do
- {
- lp=(l = *out + start + row++ * is1 ) + s1;
- do *l++ += i; while(l<lp);
- }
- while(row<s0);
- }
-}
-
-
-
-
-
-/* Put the values of the check array(s) into an array the size of the
- input image. Note that the check arrays are only the size of the
- number of meshs, not the actual input image size. */
-void
-gal_mesh_check_garray(struct gal_mesh_params *mp, float **out1,
- float **out2)
-{
- int ngarrays=mp->ngarrays;
- size_t gid, row, start, chbasedid, *types=mp->types;
- size_t s0, s1, is1=mp->s1, *ts0=mp->ts0, *ts1=mp->ts1;
- float *f, *fp, *ff=NULL, garray1=FLT_MAX, garray2=FLT_MAX;
-
- /* Allocate the array to keep the mesh indexs. Calloc is used so we
- can add all the indexes to the existing value to make sure that
- there is no overlap. */
- errno=0; *out1=malloc(mp->s0*mp->s1*sizeof **out1);
- if(*out1==NULL)
- error(EXIT_FAILURE, errno,
- "%zu bytes for out1 in gal_mesh_check_garray (mesh.c)",
- mp->s0*mp->s1*sizeof **out1);
- if(ngarrays==2)
- {
- errno=0; *out2=malloc(mp->s0*mp->s1*sizeof **out2);
- if(*out2==NULL)
- error(EXIT_FAILURE, errno,
- "%zu bytes for out2 in gal_mesh_check_garray (mesh.c)",
- mp->s0*mp->s1*sizeof **out2);
- }
-
- /* Fill the array: */
- for(gid=0;gid<mp->nmeshi;++gid)
- {
- /* Set the proper meshid depending on what garray points to, see
- the explanation above setmeshid. */
- chbasedid = gal_mesh_ch_based_id_from_gid(mp, gid);
-
- /* Fill the output array with the value in this mesh. It is
- really important that `i' should be used for the garrays, not
- cmeshid. cmeshid is only for the basic mesh parameters that
- it is used for. */
- row=0;
- s0=ts0[types[chbasedid]];
- s1=ts1[types[chbasedid]];
- start=mp->start[chbasedid];
- garray1 = mp->garray1[gid];
- if(ngarrays==2)
- garray2 = mp->garray2[gid];
- do
- {
- fp= ( f = *out1 + start + row * is1 ) + s1;
- if(ngarrays==2) ff= *out2 + start + row * is1;
- do
- {
- *f++ = garray1;
- if(ngarrays==2) *ff++ = garray2;
- }
- while(f<fp);
- ++row;
- }
- while(row<s0);
- }
-}
-
-
-
-
-
-/* Save the mesh grid values into an output file. */
-void
-gal_mesh_value_file(struct gal_mesh_params *mp, char *filename,
- char *extname1, char *extname2, struct wcsprm *wcs,
- char *spack_string)
-{
- size_t dsize[2];
- gal_data_t data;
- float *tmp1=NULL, *tmp2=NULL;
-
- /* Set the basic shared properties of the dataset. Note that the float
- image is written as a float FITS file, so anyblank is irrelevant, so
- is `size'.*/
- data.ndim=2;
- data.dsize=dsize;
- data.type=GAL_DATA_TYPE_FLOAT32;
-
- if(mp->meshbasedcheck)
- {
- /* We want one pixel per mesh. If the last operation was on
- cgarray, then first the fgarray has to be filled. Note that
- when more than one channel is present, only fgarray can be
- used for this job. In cgarray the meshs are ordered
- differently. */
- if(mp->garray1==mp->cgarray1) gal_mesh_full_garray(mp, 0);
- data.wcs=NULL; /* Not the original image size, to have same WCS */
- data.name=extname1;
- data.array=mp->fgarray1;
- data.dsize[0]=mp->gs0*mp->nch2;
- data.dsize[1]=mp->gs1*mp->nch1;
- gal_fits_img_write(&data, filename, NULL, spack_string);
- if(mp->ngarrays==2)
- {
- /* Note that gal_mesh_full_garray will correct both the meshs if
- there are two.*/
- data.name=extname2;
- data.array=mp->fgarray2;
- gal_fits_img_write(&data, filename, NULL, spack_string);
- }
- }
- else
- {
- gal_mesh_check_garray(mp, &tmp1, &tmp2);
- data.wcs=wcs;
- data.array=tmp1;
- data.name=extname1;
- data.dsize[0]=mp->s0;
- data.dsize[1]=mp->s1;
- gal_fits_img_write(&data, filename, NULL, spack_string);
- if(mp->ngarrays==2)
- {
- data.array=tmp2;
- data.name=extname2;
- gal_fits_img_write(&data, filename, NULL, spack_string);
- }
- free(tmp1);
- free(tmp2);
- }
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/*********************************************************************/
-/******************** Creating the mesh ********************/
-/*********************************************************************/
-/* In each channel we can have four types (sizes of meshs):
-
- 0. If meshsize is chosen nicely in relation to the channel size,
- the majority (or ideally all) of the mesh tiles are going to
- have an area of meshsize*meshsize. These are built from pixel
- (0,0) until ((s0/mesh_size)*mesh_size,
- (s1/mesh_size)*mesh_size), note that these are integer
- divisions, not floating point, so 5/2=2.
-
- 1. The top row of meshes. If s0%meshsize is larger than a certain
- fraction of meshsize, they will have s0%meshsize pixels in each
- row. If not, they will have meshsize+s0%meshsize rows. In any
- case, all but the last (type 3) will have meshsize columns. The
- number of columns of the last mesh in the top row, will depend
- on the third type of meshes.
-
- 2. Just like type 2, but for the width or the number of columns in
- the image.
-
- 3. The size of the last (top left) mesh in the grid will get its
- height (number of rows) from type 2s and its width from type 3s.
-
- Note that the numbers for the types used here are only for when all
- four types exist, when there is only one or two types, the first
- (in the order above) gets type id of zero and the second gets type
- id of 1. In the actual function, the indexs are given in the order
- the types are found.
- */
-static void
-fillmeshinfo(struct gal_mesh_params *mp, size_t chs0, size_t chs1,
- size_t lasts0, size_t lasts1)
-{
- size_t i, j, chi, chj, gs0=mp->gs0, gs1=mp->gs1;
- size_t meshsize=mp->meshsize, *chindex=mp->chindex;
- size_t numtypes=0, totalmeshcount=0, chid, s1=mp->s1;
- size_t meshid, typeind, lasti, lastj, nmeshc=mp->nmeshc;
- size_t nch1=mp->nch1, nch2=mp->nch2, *ts0=mp->ts0, *ts1=mp->ts1;
- size_t *types=mp->types, *start=mp->start, *imgindex=mp->imgindex;
-
- /* Initialize the maximum number of rows and columns: */
- mp->maxs1=mp->maxs0=0;
-
- /* Main meshs (type 0) in each channel. If we have more than one row
- or column of meshes in each channel, then we will have some type1
- meshes.*/
- if(gs0>1 || gs1>1)
- {
- typeind=numtypes++;
- ts1[typeind]=ts0[typeind]=meshsize;
- lasti = lasts0==meshsize ? gs0 : gs0-1;
- lastj = lasts1==meshsize ? gs1 : gs1-1;
- if(ts0[typeind]>mp->maxs0) mp->maxs0=ts0[typeind];
- if(ts1[typeind]>mp->maxs1) mp->maxs1=ts1[typeind];
-
- for(chi=0;chi<nch2;++chi) /* Don't forget that C and FITS */
- for(chj=0;chj<nch1;++chj) /* axises have different order. */
- {
- chid=chi*nch1+chj;
- for(i=0;i<lasti;++i)
- for(j=0;j<lastj;++j)
- {
- ++totalmeshcount;
- meshid=chid*nmeshc + i*gs1+j;
- types[meshid]=typeind;
- chindex[meshid]=i*gs1+j;
- imgindex[meshid]=( (chi*gs0+i) * nch1*gs1 + chj*gs1+j );
- start[meshid]=( (chi*chs0+i*meshsize) * s1
- + (chj*chs1+j*meshsize) );
- }
- }
- }
-
- /* Top row of meshes (type 1) in each channel. Note that when there
- is only one column of meshes, then the only one mesh remaining
- will be type 4, not type 1 (here).*/
- if(gs1>1 && lasts0!=meshsize)
- {
- typeind=numtypes++;
- ts0[typeind] = lasts0;
- ts1[typeind] = meshsize;
- lastj = lasts1==meshsize ? gs1 : gs1-1;
- if(ts0[typeind]>mp->maxs0) mp->maxs0=ts0[typeind];
- if(ts1[typeind]>mp->maxs1) mp->maxs1=ts1[typeind];
-
- for(chi=0;chi<nch2;++chi)
- for(chj=0;chj<nch1;++chj)
- {
- chid=chi*nch1+chj;
- for(j=0;j<lastj;++j)
- {
- ++totalmeshcount;
- meshid=chid*nmeshc + (gs0-1)*gs1+j;
- types[meshid]=typeind;
- chindex[meshid]=(gs0-1)*gs1+j;
- imgindex[meshid]=(chi*gs0+gs0-1) * nch1*gs1 + (chj*gs1+j);
- start[meshid] = ( (chi*chs0+(gs0-1)*meshsize) * s1
- + (chj*chs1+j*meshsize) );
- }
- }
- }
-
- /* Left column of meshes (type 2) in each channel. When there is
- only one row of meshs, then the last single remaining mesh will
- be type 3, not type 2 (here).*/
- if(gs0>1 && lasts1!=meshsize)
- {
- typeind=numtypes++;
- ts0[typeind] = meshsize;
- ts1[typeind] = lasts1;
- lasti = lasts0==meshsize ? gs0 : gs0-1;
- if(ts0[typeind]>mp->maxs0) mp->maxs0=ts0[typeind];
- if(ts1[typeind]>mp->maxs1) mp->maxs1=ts1[typeind];
-
- for(chi=0;chi<nch2;++chi)
- for(chj=0;chj<nch1;++chj)
- {
- chid=chi*nch1+chj;
- for(i=0;i<lasti;++i)
- {
- ++totalmeshcount;
- meshid=chid*nmeshc + i*gs1+(gs1-1);
- types[meshid]=typeind;
- chindex[meshid]=i*gs1+(gs1-1);
- imgindex[meshid]=( (chi*gs0+i) * nch1*gs1
- + (chj*gs1+(gs1-1)) );
- start[meshid] = ( (chi*chs0+i*meshsize) * s1
- + (chj*chs1+(gs1-1)*meshsize) );
- }
- }
- }
-
- /* Top left mesh or only mesh (type 3) in each channel. Note that
- this might only happen for one mesh in each channel (if there are
- 4 types of meshes). */
- if(mp->nmeshi-totalmeshcount==mp->nch)
- {
- typeind=numtypes++;
- ts0[typeind] = lasts0;
- ts1[typeind] = lasts1;
- if(ts0[typeind]>mp->maxs0) mp->maxs0=ts0[typeind];
- if(ts1[typeind]>mp->maxs1) mp->maxs1=ts1[typeind];
-
- for(chi=0;chi<nch2;++chi)
- for(chj=0;chj<nch1;++chj)
- {
- ++totalmeshcount;
- meshid=(chi*nch1+chj)*nmeshc + gs0*gs1-1;
- types[meshid]=typeind;
- chindex[meshid]=gs0*gs1-1;
- imgindex[meshid]=( (chi*gs0+(gs0-1)) * nch1*gs1
- + (chj*gs1+(gs1-1)) );
- start[meshid]=( (chi*chs0+(gs0-1)*meshsize) * s1
- + (chj*chs1+(gs1-1)*meshsize) );
- }
- }
-
- /* Just for a check: */
- if(totalmeshcount!=mp->nmeshi)
- error(EXIT_FAILURE, 0, "a bug! Please contact us at %s so we can fix "
- "it. The basic information for some meshes has not been found "
- "(in fillmeshinfo of mesh.c)", PACKAGE_BUGREPORT);
-}
-
-
-
-
-
-/* In the explanation below, all the parameters named are from the
- `struct gal_mesh_params` of `mesh.h`.
-
- An image contians `s0*s1` pixels in s0 rows and s1 columns. A mesh
- is a box of pixels within the image. A channel is defined as large
- regions of the image that should contain the meshes discussed
- above. So each channel contains a certain number of meshs and each
- mesh contains a certain number of pixels.
-
- In general, for any channel size, there can be four types of
- meshes. See the explanations above fillmeshinfo() for a more
- detailed discussion. So when all the channels have the same size,
- in the whole image, there can only be four types of meshes.
-
- The idea behind this classification is that CCDs often have more
- than one reading channel and essentially all the properties, vary
- from channel to channel. Therefore it is very important to do
- convolution, thresholding, detection and finding the sky background
- and noise properties differently for each channel.
-
- Each mesh is treated independently in its thread, but when
- interpolation over the meshes is desired, only meshes in each
- channel should be interpolated over.
-*/
-void
-gal_mesh_make_mesh(struct gal_mesh_params *mp)
-{
- size_t meshsize=mp->meshsize, lasts0, lasts1;
- size_t i, chs0=mp->s0/mp->nch2, chs1=mp->s1/mp->nch1;
-
- /* Set the number of channels. */
- mp->nch=mp->nch1*mp->nch2;
-
- /* Incase meshsize is larger than the channel sizes, make it equal
- to the smaller of the channel axis sizes. */
- if(meshsize>chs0 || meshsize>chs1)
- mp->meshsize = meshsize = (chs0<chs1
- ? (chs0%2 ? chs0-1 : chs0)
- : (chs1%2 ? chs1-1 : chs1));
-
- /* Initialize the necessary arrays: */
- for(i=0;i<4;++i) mp->ts0[i]=mp->ts1[i]=0;
-
- /* Set the mesh grid height (gs0, number of rows) and width (gs1,
- number of columns) on each channel. Note that the channels are
- assumed to be identical in size, so these values are the same for
- all the channels. In case the remainder of the image side size
- and meshsize is larger than lastmeshfrac, add a new mesh. If it
- is smaller than that fraction, then the remainder should be
- merged with the last mesh in the row or column. */
- if( (float)(chs0%meshsize) > mp->lastmeshfrac*(float)(meshsize) )
- {
- mp->gs0 = chs0/meshsize + 1;
- lasts0 = chs0%meshsize;
- }
- else
- {
- mp->gs0 = chs0/meshsize;
- lasts0 = meshsize+chs0%meshsize;
- }
- if( (float)(chs1%meshsize) > mp->lastmeshfrac*(float)(meshsize) )
- {
- mp->gs1 = chs1/meshsize + 1;
- lasts1 = chs1%meshsize;
- }
- else
- {
- mp->gs1 = chs1/meshsize;
- lasts1 = meshsize+chs1%meshsize;
- }
-
- mp->nmeshc=mp->gs0*mp->gs1;
- mp->nmeshi=mp->nmeshc*mp->nch;
-
- /* Allocate the arrays to keep all the mesh starting points and
- types. Irrespective of which channel they lie in. */
- mp->cgarray1=mp->cgarray2=mp->fgarray1=mp->fgarray2=NULL;
- errno=0; mp->start=malloc(mp->nmeshi*sizeof *mp->start);
- if(mp->start==NULL) error(EXIT_FAILURE, errno, "mesh starting points");
- errno=0; mp->types=malloc(mp->nmeshi*sizeof *mp->types);
- if(mp->types==NULL) error(EXIT_FAILURE, errno, "mesh types");
- errno=0; mp->chindex=malloc(mp->nmeshi*sizeof *mp->chindex);
- if(mp->chindex==NULL) error(EXIT_FAILURE, errno, "mesh in channel index");
- errno=0; mp->imgindex=malloc(mp->nmeshi*sizeof *mp->imgindex);
- if(mp->imgindex==NULL) error(EXIT_FAILURE, errno, "mesh in image index");
-
- /* Distribute the meshes in all the threads. */
- gal_threads_dist_in_threads(mp->nmeshi, mp->numthreads, &mp->indexs,
- &mp->thrdcols);
-
- /* Fill in the information for each mesh and each type. */
- fillmeshinfo(mp, chs0, chs1, lasts0, lasts1);
-}
-
-
-
-
-
-void
-gal_mesh_free_mesh(struct gal_mesh_params *mp)
-{
- free(mp->start);
- free(mp->types);
- free(mp->indexs);
- free(mp->chindex);
- free(mp->cgarray1);
- free(mp->cgarray2);
- free(mp->fgarray1);
- free(mp->fgarray2);
- free(mp->imgindex);
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/*********************************************************************/
-/******************** Mesh operations ********************/
-/*********************************************************************/
-/* The purpose of this function is to prepare the (possibly)
- multi-threaded environemnt, spin-off each thread and wait for them
- to finish for any operation that is to be done on the mesh
- grid.
-
- The arguments are:
-
- 1. A pointer to the gal_mesh_params structure that keeps all the
- information.
-
- 2. A pointer to a function that returns and gets a `void *' as its only
- argument. This function will be directly given to
- pthread_create. Through this argument, you can choose the function to
- operate on the mesh grid.
-
- 3. The size of each element to copy the mesh grid into, this has to
- be type size of the same type that constitutes `img' in
- gal_mesh_params. If the value to this argument is non-zero, an array
- will be allocated that can contain all the pixels in all the
- meshs and can be used by threads to manipute the pixels (for
- example sort them) in each mesh.
-
- 4. If the value of this argument is 1, then a second garray will be
- allocated in case your operation needs it.
-
- 5. Wether the allocated garrays should be initialized or not.
-*/
-void
-gal_mesh_operate_on_mesh(struct gal_mesh_params *mp,
- void *(*meshfunc)(void *), size_t oneforallsize,
- int makegarray2, int initialize)
-{
- int err;
- size_t i, nb;
- float *f, *fp;
- pthread_t t; /* We don't use the thread id, so all are saved here. */
- pthread_attr_t attr;
- struct gal_mesh_thread_params *mtp;
- size_t numthreads=mp->numthreads;
-
- /* Allocate the arrays to keep the thread and parameters for each
- thread. */
- errno=0; mtp=malloc(numthreads*sizeof *mtp);
- if(mtp==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes in fillmesh (mesh.c) for mtp",
- numthreads*sizeof *mtp);
-
- /* Set the number of garrays to operate on: */
- mp->ngarrays = makegarray2 ? 2 : 1;
-
- /* If cgarrays have not been allocated before, allocate them. */
- if(mp->cgarray1==NULL)
- {
- errno=0; mp->cgarray1=malloc(mp->nmeshi*sizeof *mp->cgarray1);
- if(mp->cgarray1==NULL) error(EXIT_FAILURE, errno, "mp->cgarray1");
- }
- if(mp->ngarrays==2 && mp->cgarray2==NULL)
- {
- errno=0; mp->cgarray2=malloc(mp->nmeshi*sizeof *mp->cgarray2);
- if(mp->cgarray2==NULL) error(EXIT_FAILURE, errno, "mp->cgarray2");
- }
-
- /* Initialize the cgarrays to NaN:*/
- if(initialize)
- {
- fp=(f=mp->cgarray1)+mp->nmeshi; do *f++=NAN; while(f<fp);
- if(mp->ngarrays==2)
- { fp=(f=mp->cgarray2)+mp->nmeshi; do *f++=NAN; while(f<fp); }
- mp->garray1=mp->cgarray1;
- mp->garray2=mp->cgarray2;
- }
-
- /* `oneforall' is an array with the sides equal to the maximum side
- of the meshes in the image. The purpose is to enable manipulating
- the mesh pixels (for example sorting them and so on.). One such
- array is allocated for each thread in this one full
- allocation. Each thread can then use its own portion of this
- array through the following declaration:
-
- float *oneforall=&mp->oneforall[mtp->id*mp->maxs0*mp->maxs1];
-
- In gal_mesh_params, `oneforall' is defined as a `void *', so the
- caller function, can cast it to any type it wants. The size of
- each type is given to `fillmesh' through the `oneforallsize'
- argument.
- */
- if(oneforallsize)
- {
- errno=0;
- mp->oneforall=malloc(numthreads*mp->maxs0*mp->maxs1*oneforallsize);
- if(mp->oneforall==NULL)
- error(EXIT_FAILURE, errno, "unable to allocate %zu bytes for"
- "mtp->oneforall in fillmesh of mesh.c",
- numthreads*mp->maxs0*mp->maxs1*oneforallsize);
- }
-
- /* Spin off the threads: */
- if(numthreads==1)
- {
- mtp[0].id=0;
- mtp[0].mp=mp;
- (*meshfunc)(&mtp[0]);
- }
- else
- {
- /* Initialize the attributes. Note that this running thread
- (that spinns off the nt threads) is also a thread, so the
- number the barrier should be one more than the number of
- threads spinned off. */
- if(mp->nmeshi<numthreads) nb=mp->nmeshi+1;
- else nb=numthreads+1;
- gal_threads_attr_barrier_init(&attr, &mp->b, nb);
-
- /* Spin off the threads: */
- for(i=0;i<numthreads;++i)
- if(mp->indexs[i*mp->thrdcols]!=GAL_THREADS_NON_THRD_INDEX)
- {
- mtp[i].id=i;
- mtp[i].mp=mp;
- err=pthread_create(&t, &attr, meshfunc, &mtp[i]);
- if(err) error(EXIT_FAILURE, 0, "can't create thread %zu", i);
- }
-
- /* Wait for all threads to finish and free the spaces. */
- pthread_barrier_wait(&mp->b);
- pthread_attr_destroy(&attr);
- pthread_barrier_destroy(&mp->b);
- }
-
- free(mtp);
- if(oneforallsize) free(mp->oneforall);
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/*********************************************************************/
-/******************** Interpolate ********************/
-/*********************************************************************/
-static void
-preparemeshinterparrays(struct gal_mesh_params *mp)
-{
- size_t bs0=mp->gs0, bs1=mp->gs1;
- size_t numthreads=mp->numthreads;
-
-
- /* If full-interpolation is to be done (ignoring channels) we will
- need two arrays to temporarily keep the actual positions of the
- meshs. Note that by default, the meshs in each channels are
- placed contiguously. In these arrays, the meshs are placed as
- they would in the image (channels are ignored).*/
- if(mp->fullinterpolation)
- {
- /* In case the previous operation was on cgarrays, then you have
- to fill in fgarray. */
- if(mp->garray1==mp->cgarray1)
- gal_mesh_full_garray(mp, 0);
- bs0=mp->nch2*mp->gs0;
- bs1=mp->nch1*mp->gs1;
- mp->garray1=mp->fgarray1;
- if(mp->ngarrays==2) mp->garray2=mp->fgarray2;
- }
- else
- {
- mp->garray1=mp->cgarray1;
- if(mp->ngarrays==2) mp->garray2=mp->cgarray2;
- }
-
-
- /* Allocate the output arrays to keep the final values. Note that we
- cannot just do the interpolaton on the same input grid, because
- the newly filled interpolated pixels will affect the later
- ones. In the end, these copies are going to replace the
- garrays. */
- errno=0; mp->outgarray1=malloc(mp->nmeshi*sizeof *mp->outgarray1);
- if(mp->outgarray1==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for outgarray1 (mesh.c)",
- mp->nmeshi*sizeof *mp->outgarray1);
- if(mp->ngarrays==2)
- {
- errno=0; mp->outgarray2=malloc(mp->nmeshi*sizeof *mp->outgarray2);
- if(mp->outgarray2==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for outgarray2 (mesh.c)",
- mp->nmeshi*sizeof *mp->outgarray2);
- }
-
-
- /* Allocate the array for the values of the nearest pixels. For each
- pixel in each thread, the nearest pixels will be put here and
- sorted to find the median value. If garray2 is not NULL, then it
- should be interpolated. Note that both arrays have blank pixels
- on the same places.*/
- errno=0;
- mp->nearest1=malloc(numthreads*mp->numnearest*sizeof *mp->nearest1);
- if(mp->nearest1==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for the array to keep the "
- "nearest1 values for interpolation (mesh.c)",
- numthreads*mp->numnearest*sizeof *mp->nearest1);
- if(mp->ngarrays==2)
- {
- errno=0;
- mp->nearest2=malloc(numthreads*mp->numnearest*sizeof *mp->nearest2);
- if(mp->nearest2==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for the array to keep the "
- "nearest2 values for interpolation (mesh.c)",
- numthreads*mp->numnearest*sizeof *mp->nearest2);
- }
-
-
- /* Allocate space for the byte array keeping a record of which
- pixels were used to search. Note that each thread needs one byt
- array. */
- errno=0; mp->byt=malloc(numthreads*bs0*bs1*sizeof *mp->byt);
- if(mp->byt==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for mp->byt array in "
- "mesh.c", numthreads*bs0*bs1*sizeof *mp->byt);
-}
-
-
-
-
-
-/* Ideally, this should be a radial distance using the square root of
- the differences. But here we are not dealing with subpixel
- distances, so manhattan distance is fine. It also avoids the square
- root function which is much slower. */
-static float
-manhattandistance(long ind, long xc, long yc, long s1)
-{
- return labs(ind%s1 - xc) + labs(ind/s1 - yc);
-}
-
-
-
-
-/* Some of the variables have different names than the gal_mesh_params
- structure because they are to be fed into the
- GAL_NEIGHBORS_FILL_4_ALLIMG macro. */
-static void *
-meshinterponthread(void *inparams)
-{
- struct gal_mesh_thread_params *mtp=(struct gal_mesh_thread_params *)inparams;
- struct gal_mesh_params *mp=mtp->mp;
-
- /* Basic variables used in other definitions: */
- size_t numnearest=mp->numnearest;
- size_t is0=mp->fullinterpolation ? mp->gs0*mp->nch2 : mp->gs0;
- size_t is1=mp->fullinterpolation ? mp->gs1*mp->nch1 : mp->gs1;
-
- /* Variables for this function: */
- int ngarrays=mp->ngarrays;
- struct gal_linkedlist_tosll *lQ, *sQ;
- size_t xc, yc, *n, *nf, currentnum, thisind;
- unsigned char *byt=&mp->byt[mtp->id*is0*is1];
- float *nearest1=&mp->nearest1[mtp->id*numnearest];
- size_t i, *indexs=&mp->indexs[mtp->id*mp->thrdcols];
- float mdist, *garray1=mp->garray1, *garray2=mp->garray2;
- size_t fmeshid, position, *ind=&position, numngb, ngb[4];
- float *outgarray1=mp->outgarray1, *outgarray2=mp->outgarray2;
- float *nearest2 = ngarrays==2 ? &mp->nearest2[mtp->id*numnearest] : NULL;
-
- /* Go over all the meshes for this thread. */
- for(i=0;indexs[i]!=GAL_THREADS_NON_THRD_INDEX;++i)
- {
- /* Get the index of this NaN mesh: */
- thisind=indexs[i];
-
- /* If this mesh is not blank and the user has only asked to
- interpolate blank pixels, then set the final values and go
- onto the next mesh. */
- if(mp->interponlyblank && !isnan(garray1[thisind]))
- {
- outgarray1[thisind]=garray1[thisind];
- if(ngarrays==2) outgarray2[thisind]=garray2[thisind];
- continue;
- }
-
- /* Reset the byt array: */
- memset(byt, 0, is0*is1);
-
- /*
- `fmeshid' (first-mesh-id) keeps the index of the first mesh
- in this channel. We are going to consider each channel as a
- sparate image for the neighbor finding job. This is why the
- pixels in each channel were designed to be contiguous. When
- we want to look at the value on each mesh, we simply sum its
- channel index with fmeshid. From this point on, *ind
- corresponds to the index of the mesh in the channel.
-
- We will also only be concerened with the portion of `mp->byt'
- array that is within this channel. So the `byt' pointer is
- set to point to the first mesh index in this channel.
- */
- lQ=sQ=NULL;
- fmeshid = ( mp->fullinterpolation
- ? 0
- : (thisind/mp->nmeshc) * mp->nmeshc);
-
- /* Initialize the necessary parameters. */
- *ind=thisind-fmeshid;
- xc=*ind%is1;
- yc=*ind/is1;
- byt[*ind]=1;
- currentnum=0;
- gal_linkedlist_add_to_tosll_end( &lQ, &sQ, *ind, 0 );
-
- /* Start finding the nearest filled pixels. */
- while(sQ)
- {
- /* Pop out a pixel index (p) from the queue: */
- gal_linkedlist_pop_from_tosll_start(&lQ, &sQ, ind, &mdist);
-
- /* If it isn't a NaN, then put it in the `nearest1' and
- `nearest2' arrays. */
- if(!isnan(garray1[*ind+fmeshid]))
- {
- nearest1[currentnum]=garray1[*ind+fmeshid];
- if(ngarrays==2) nearest2[currentnum]=garray2[*ind+fmeshid];
- if(++currentnum>=numnearest) break;
- }
-
- /* Check the four neighbors and if they have not already
- been checked, put them into the queue. */
- GAL_NEIGHBORS_FILL_4_ALLIMG;
-
- /* It might happen that there are no neighbors (e.g., that there
- is only one mesh). In that case, we shouldn't look into the
- neighbors.*/
- n=ngb;
- nf = numngb ? n+numngb : n;
- while(n<nf)
- {
- if(byt[*n]==0)
- {
- byt[*n]=1;
- gal_linkedlist_add_to_tosll_end(&lQ, &sQ, *n,
- manhattandistance(*n, xc,
- yc, is1));
- }
- ++n;
- }
-
- /* If there are no more meshes to add to the queue, then
- this shows, there were not enough points for
- interpolation. Normally, this loop should only be exited
- through the `currentnum>=numnearest' check above. */
- if(sQ==NULL)
- error(EXIT_FAILURE, 0, "%s: only %zu mesh(s) are filled for "
- "interpolation in channel %zu. Either set less "
- "restrictive requirements to get more interpolation "
- "points or decrease the number of nearest points to "
- "use for interpolation. Problem encountered on thread "
- "%zu, for pixel %zu. When running on multiple threads, "
- "This message might be repeated for different threads",
- mp->errstart, currentnum, thisind/mp->nmeshc, mtp->id,
- thisind);
- }
- gal_linkedlist_tosll_free(lQ); /* Rest of the queue not needed. */
-
-
- /* Find the median of the nearest neighbors and put it in: */
- qsort(nearest1, numnearest, sizeof *nearest1,
- gal_qsort_float32_increasing);
- outgarray1[thisind] = ( numnearest%2 ?
- nearest1[numnearest/2] : /* Odd. */
- (nearest1[numnearest/2] /* Even. */
- +nearest1[numnearest/2-1])/2 );
- if(ngarrays==2)
- {
- qsort(nearest2, numnearest, sizeof *nearest2,
- gal_qsort_float32_increasing);
- outgarray2[thisind] = ( numnearest%2 ?
- nearest2[numnearest/2] : /* Odd. */
- (nearest2[numnearest/2] /* Even. */
- +nearest2[numnearest/2-1])/2 );
- }
- }
-
- /* If there is more than one thread, wait until the others
- finish. */
- if(mp->numthreads>1)
- pthread_barrier_wait(&mp->b);
- return NULL;
-}
-
-
-
-
-
-void
-gal_mesh_interpolate(struct gal_mesh_params *mp, char *errstart)
-{
- int err;
- pthread_t t; /* We don't use the thread id, so all are saved here. */
- size_t i, nb;
- pthread_attr_t attr;
- struct gal_mesh_thread_params *mtp;
- size_t numthreads=mp->numthreads;
-
- /* Prepare all the gal_mesh_params arrays: */
- mp->errstart=errstart;
- preparemeshinterparrays(mp);
-
- /* Allocate the arrays to keep the thread and parameters for each
- thread. */
- errno=0; mtp=malloc(numthreads*sizeof *mtp);
- if(mtp==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes in fillmesh (mesh.c) for mtp",
- numthreads*sizeof *mtp);
-
- /* Spin off the threads: */
- if(numthreads==1)
- {
- mtp[0].id=0;
- mtp[0].mp=mp;
- meshinterponthread(&mtp[0]);
- }
- else
- {
- /* Initialize the attributes. Note that this running thread
- (that spinns off the nt threads) is also a thread, so the
- number the barrier should be one more than the number of
- threads spinned off. */
- if(mp->nmeshi<numthreads) nb=mp->nmeshi+1;
- else nb=numthreads+1;
- gal_threads_attr_barrier_init(&attr, &mp->b, nb);
-
- /* Spin off the threads: */
- for(i=0;i<numthreads;++i)
- if(mp->indexs[i*mp->thrdcols]!=GAL_THREADS_NON_THRD_INDEX)
- {
- mtp[i].id=i;
- mtp[i].mp=mp;
- err=pthread_create(&t, &attr, meshinterponthread, &mtp[i]);
- if(err) error(EXIT_FAILURE, 0, "can't create thread %zu", i);
- }
-
- /* Wait for all threads to finish and free the spaces. */
- pthread_barrier_wait(&mp->b);
- pthread_attr_destroy(&attr);
- pthread_barrier_destroy(&mp->b);
- }
-
- /* Replace garray1 and garray2 with outgarray1 and outgarray2 */
- if(mp->fullinterpolation)
- { free(mp->fgarray1); mp->garray1=mp->fgarray1=mp->outgarray1; }
- else
- { free(mp->cgarray1); mp->garray1=mp->cgarray1=mp->outgarray1; }
- if(mp->ngarrays==2)
- {
- if(mp->fullinterpolation)
- { free(mp->fgarray2); mp->garray2=mp->fgarray2=mp->outgarray2; }
- else
- { free(mp->cgarray2); mp->garray2=mp->cgarray2=mp->outgarray2; }
- }
- mp->outgarray1=mp->outgarray2=NULL;
-
- /* For a check
- system("rm test.fits");
- gal_fits_array_to_file("test.fits", "garray1", FLOAT_IMG, mp->garray1,
- mp->nch2*mp->gs0, mp->nch1*mp->gs1, 1, NULL, NULL,
- "mesh");
- gal_fits_array_to_file("test.fits", "garray2", FLOAT_IMG, mp->garray2,
- mp->nch2*mp->gs0, mp->nch1*mp->gs1, 1, NULL, NULL,
- "mesh");
- */
-
- /* Clean up. */
- free(mtp);
- free(mp->byt);
- free(mp->nearest1);
- if(mp->ngarrays==2) free(mp->nearest2);
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/*********************************************************************/
-/******************** Smooth ********************/
-/*********************************************************************/
-void
-gal_mesh_smooth(struct gal_mesh_params *mp)
-{
- float *charray;
- float *f, *o, *fp, *tmp, *kernel, *sgarray1, *sgarray2;
- size_t numthreads=mp->numthreads, gs0=mp->gs0, gs1=mp->gs1;
- size_t chid, nmeshc=mp->nmeshc, smoothwidth=mp->smoothwidth;
-
- /* Make the smoothing kernel and set all its elements to 1. */
- errno=0;
- kernel=malloc(smoothwidth*smoothwidth*sizeof *kernel);
- if(kernel==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for kernel in mesh.c",
- mp->smoothwidth*sizeof *kernel);
- fp=(f=kernel)+smoothwidth*smoothwidth; do *f++=1; while(f<fp);
-
- /* Smooth the garrays. */
- if(mp->fullsmooth)
- {
- /* In case the previous operation was on the cgarrays. garray2
- should not be checked. garray1 and garray2 are always in sync
- during one operation. However, if a previous operation had
- set garray2 and this operation only uses garray1, then things
- are going to go wrong. So Based on the principle that garray1
- and garray2 are always going to be in sync with each other,
- we don't the mp->garray2==mp->cgarray2 check should not be
- done. */
- if(mp->garray1==mp->cgarray1)
- gal_mesh_full_garray(mp, 0);
-
- /* Do the spatial convolution */
- gal_spatialconvolve_convolve(mp->fgarray1, gs0*mp->nch2, gs1*mp->nch1,
- kernel, smoothwidth, smoothwidth,
- numthreads, 1, &sgarray1);
-
- free(mp->fgarray1);
- mp->garray1=mp->fgarray1=sgarray1;
- if(mp->ngarrays==2)
- {
- gal_spatialconvolve_convolve(mp->fgarray2, gs0*mp->nch2,
- gs1*mp->nch1, kernel, smoothwidth,
- smoothwidth, mp->numthreads,
- 1, &sgarray2);
- free(mp->fgarray2);
- mp->garray2=mp->fgarray2=sgarray2;
- }
- }
- else
- for(chid=0;chid<mp->nch;++chid) /* Note that in this mode, we do not */
- { /* Allocate anything :-). */
-
- /* If the last operation was done on the full image, then
- mp->garray1==mp->fgarray1. Therefore, the mesh boxes in
- each channel will not be congituous. So we have to update
- cgarray and set mp->garray1=mp->cgarray1. */
- if(mp->garray1==mp->fgarray1)
- gal_mesh_full_garray(mp, 1);
- mp->garray1=mp->cgarray1;
- mp->garray2=mp->cgarray2;
-
- charray=&mp->cgarray1[chid*nmeshc];
- gal_spatialconvolve_convolve(charray, gs0, gs1, kernel, smoothwidth,
- smoothwidth, numthreads, 1, &tmp);
- o=tmp; fp=(f=charray)+gs0*gs1; do *f=*o++; while(++f<fp);
- free(tmp);
- if(mp->ngarrays==2)
- {
- charray=&mp->cgarray2[chid*nmeshc];
- gal_spatialconvolve_convolve(charray, gs0, gs1, kernel,
- smoothwidth, smoothwidth,
- numthreads, 1, &tmp);
- o=tmp; fp=(f=charray)+gs0*gs1; do *f=*o++; while(++f<fp);
- free(tmp);
- }
- }
- /* Clean up: */
- free(kernel);
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/*********************************************************************/
-/******************** Convolve in mesh ********************/
-/*********************************************************************/
-static void*
-meshspatialconvonthreads(void *inparam)
-{
- struct gal_mesh_thread_params *mtp = (struct gal_mesh_thread_params
*)inparam;
- struct gal_mesh_params *mp = mtp->mp;
- const size_t ks0=mp->ks0, ks1=mp->ks1, is1=mp->s1;
-
- double sum, ksum;
- const size_t *types=mp->types;
- size_t chid, *chbrd=mtp->chbrd, ystart;
- const size_t *ts0=mp->ts0, *ts1=mp->ts1;
- size_t i, meshid, ind, *start=mp->start;
- size_t x, y, fx, fy, a, b, k0h=ks0/2, k1h=ks1/2;
- size_t xmin, xmax, ymin, ymax, nmeshc=mp->nmeshc;
- size_t *indexs=&mp->indexs[mtp->id*mp->thrdcols];
- float *img=mp->img, *conv=mtp->conv, *kernel=mp->kernel;
-
- /* For each mesh in this thread, convolve the mesh */
- for(i=0; (meshid=indexs[i])!=GAL_THREADS_NON_THRD_INDEX; ++i)
- {
- /* Get the channel ID and the relevant information for this mesh. */
- chid=meshid/nmeshc;
- xmin=chbrd[chid*4]+k0h; /* Within these four points, the */
- xmax=chbrd[chid*4+2]-k0h; /* convolution does not have to worry */
- ymin=chbrd[chid*4+1]+k1h; /* about the edges. */
- ymax=chbrd[chid*4+3]-k1h;
-
- /* Get the starting and ending positions of this mesh. NOTE:
- ystart is needed because for every `x', the value of y needs
- to change. */
- fx = ( x = start[meshid]/is1 ) + ts0[types[meshid]];
- fy = ( ystart = start[meshid]%is1 ) + ts1[types[meshid]];
-
- /* If the mesh is on the edge of the channel it should be
- treated differently compared to when it is not. */
- if( x>=xmin && ystart>=ymin && fx<xmax && fy<ymax )
- { /* All pixels in this mesh are distant */
- for(;x<fx;++x) /* enough from the edge of the channel. */
- for(y=ystart;y<fy;++y)
- {
- if(isnan(img[x*is1+y]))
- conv[x*is1+y]=NAN;
- else
- {
- ksum=sum=0.0f;
- for(a=0;a<ks0;++a)
- for(b=0;b<ks1;++b)
- if( !isnan( img[ ind = (x+a-k0h) * is1 + y+b-k1h ] ) )
- {
- ksum+=kernel[a*ks1+b];
- sum+=img[ind] * kernel[a*ks1+b];
- }
- conv[x*is1+y] = sum/ksum;
- }
- }
- }
- else
- { /* Some pixels in this mesh are too */
- for(;x<fx;++x) /* close to the edge. */
- for(y=ystart;y<fy;++y)
- {
- if( isnan(img[x*is1+y]) )
- conv[x*is1+y]=NAN;
- else
- {
- ksum=sum=0.0f;
- if(x>=xmin && y>=ymin && fx<xmax && fy<ymax)
- for(a=0;a<ks0;++a)
- for(b=0;b<ks1;++b)
- { /* {} needed, because of next `else'. */
- if( !isnan(img[ ind=(x+a-k0h) * is1 + y+b-k1h ]))
- {
- ksum+=kernel[a*ks1+b];
- sum+=img[ind] * kernel[a*ks1+b];
- }
- }
- else
- {
- for(a=0;a<ks0;++a)
- if(x+a>=xmin && x+a<chbrd[chid*4+2]+k0h)
- for(b=0;b<ks1;++b)
- if(y+b>=ymin && y+b<chbrd[chid*4+3]+k1h)
- {
- if( !isnan( img[ind=(x+a-k0h)
- * is1
- + y+b-k1h]) )
- {
- ksum+=kernel[a*ks1+b];
- sum+=img[ind] * kernel[a*ks1+b];
- }
- }
- }
- conv[x*is1+y] = sum/ksum;
- }
- }
- }
- }
-
- /* Free alltype and if multiple threads were used, wait until all
- other threads finish. */
- if(mp->numthreads>1)
- pthread_barrier_wait(&mp->b);
- return NULL;
-}
-
-
-
-
-
-void
-gal_mesh_spatial_convolve_on_mesh(struct gal_mesh_params *mp, float **conv)
-{
- int err;
- pthread_t t; /* We don't use the thread id, so all are saved here. */
- pthread_attr_t attr;
- size_t i, nb, *chbrd;
- struct gal_mesh_thread_params *mtp;
- size_t numthreads=mp->numthreads;
-
-
- /* Allocate the arrays to keep the thread and parameters for each
- thread. */
- errno=0; mtp=malloc(numthreads*sizeof *mtp);
- if(mtp==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes in fillmesh (mesh.c) for mtp",
- numthreads*sizeof *mtp);
-
-
- /* Allocate space for the convolved array. */
- errno=0; *conv=malloc(mp->s0*mp->s1* sizeof **conv);
- if(*conv==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for convolution on mesh output",
- mp->s0*mp->s1* sizeof **conv);
-
-
- /* Channel borders so we can see if a mesh is inside or outside a
- channel. chbrd contians for number for each channel: in order
- they are:
-
- 1. The channel's first pixel's first C axis value.
- 2. The channel's first pixel's second C axis value.
- 3. The channel's last pixel's first C axis value.
- 4. The channel's last pixel's second C axis value.
- */
- errno=0; chbrd=malloc(mp->nch*4*sizeof *chbrd);
- if(chbrd==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for chbrd in "
- "gal_mesh_spatial_convolve_on_mesh", mp->nch*4*sizeof *chbrd);
- for(i=0;i<mp->nch;++i)
- {
- if(mp->fullconvolution)
- {
- chbrd[i*4+2]=mp->s0;
- chbrd[i*4+3]=mp->s1;
- chbrd[i*4]=chbrd[i*4+1]=0;
- }
- else
- {
- chbrd[i*4+0]=mp->start[i*mp->nmeshc]/mp->s1;
- chbrd[i*4+1]=mp->start[i*mp->nmeshc]%mp->s1;
- chbrd[i*4+2]=chbrd[i*4+0]+mp->s0/mp->nch2;
- chbrd[i*4+3]=chbrd[i*4+1]+mp->s1/mp->nch1;
- }
- }
-
-
- /* Spin off the threads: */
- if(numthreads==1)
- {
- mtp[0].id=0;
- mtp[0].mp=mp;
- mtp[0].conv=*conv;
- mtp[0].chbrd=chbrd;
- meshspatialconvonthreads(&mtp[0]);
- }
- else
- {
- /* Initialize the attributes. Note that this running thread
- (that spinns off the nt threads) is also a thread, so the
- number the barrier should be one more than the number of
- threads spinned off. */
- if(mp->nmeshi<numthreads) nb=mp->nmeshi+1;
- else nb=numthreads+1;
- gal_threads_attr_barrier_init(&attr, &mp->b, nb);
-
- /* Spin off the threads: */
- for(i=0;i<numthreads;++i)
- if(mp->indexs[i*mp->thrdcols]!=GAL_THREADS_NON_THRD_INDEX)
- {
- mtp[i].id=i;
- mtp[i].mp=mp;
- mtp[i].conv=*conv;
- mtp[i].chbrd=chbrd;
- err=pthread_create(&t, &attr, meshspatialconvonthreads, &mtp[i]);
- if(err) error(EXIT_FAILURE, 0, "can't create thread %zu", i);
- }
-
- /* Wait for all threads to finish and free the spaces. */
- pthread_barrier_wait(&mp->b);
- pthread_attr_destroy(&attr);
- pthread_barrier_destroy(&mp->b);
- }
-
- free(mtp);
- free(chbrd);
-}
-
-
-
-
-
-
-/* The indexs array for correcting the convolution on inner channel edges
- has been allocated. Note that gal_threads_dist_in_threads will
- distribute indexs from zero to numpix-1. After it, we should fill in
- all the channels.
-
- The method of filling in the indexs array with the proper indexs to
- re-convolve is very similar to the method explained below in
- gal_mesh_change_to_full_convolution, where it is explained how to count
- the number of pixels that should be re-convolved. */
-static void
-corrconvindexs(struct gal_mesh_params *mp, size_t **indexs,
- size_t *numpix, size_t *thrdcols)
-{
- size_t i, j, a, b;
- size_t numthreads=mp->numthreads;
- size_t alow, ahigh, blow, bhigh, *ind;
- size_t npch0=mp->s0/mp->nch2, npch1=mp->s1/mp->nch1;
- size_t nch1=mp->nch1, nch2=mp->nch2, ks0=mp->ks0, ks1=mp->ks1;
- size_t s0=mp->s0, s1=mp->s1, hk0=mp->ks0/2+1, hk1=mp->ks1/2+1;
-
- /* Find the number of pixels that must be convolved. Note that when
- we don't care about the edges, on each dimension, we have one
- less border than the number of channels in that dimention.
-
- On each channel's border, we want to re-convolve the ks0/2+1
- pixels before the channel edge along the first dimension. So in
- total, for each channel edge, we want ks0+1 pixels on its two
- sides. Don't forget that the kernel sides are odd.
-
- So on the first dimension, we have (ks0+1)*(nch2-1) pixels that
- should be re-convolved. Multiplying it by s1, we get the total
- number of pixels in 2D around the first axis internal edges.
-
- Along the second dimension, we have (ks1+1)*(nch1-1) pixels. But
- this time we can't just multiply by s0 to get the total number of
- 2D pixels. Because of the overlap. So we only have to multily by
- the number of rows that were not accounted in the first run,
- which is s0-(ks0+1)*(nch2-1). */
- *numpix = ( (ks0+1)*(nch2-1)*s1
- + (ks1+1)*(nch1-1)*(s0-(ks0+1)*(nch2-1)) );
-
-
- /* Distribute the indexs of the desired pixels into the indexs
- array. */
- gal_threads_dist_in_threads(*numpix, numthreads, indexs, thrdcols);
-
- ind=*indexs;
- for(i=1;i<nch2;++i) /* FIRST LOOP. */
- {
- alow = i*npch0<hk0 ? 0 : i*npch0-hk0;
- ahigh = i*npch0+hk0>s0 ? s0 : i*npch0+hk0;
-
- /* For a check
- printf("1: alow: %-4zu ahigh: %-4zu.\t(0 -- %zu)\n",
- alow, ahigh, s1);
- */
- for(a=alow; a<ahigh; ++a)
- for(b=0;b<s1;++b)
- {
- while(*ind==GAL_THREADS_NON_THRD_INDEX) ++ind;
- *ind++=a*s1+b;
- }
- }
- for(j=1;j<nch1;++j) /* SECOND LOOP */
- {
- blow = j*npch1<hk1 ? 0 : j*npch1-hk1;
- bhigh = j*npch1+hk1>s1 ? s1 : j*npch1+hk1;
- for(b=blow; b<bhigh; ++b)
- {
- /* Since there might be multiple channels along the first C axis
- and we want the spaces between their borders, alow is
- initiated with 0. */
- alow=0;
-
- /* Check the areas under the bordering area for each
- vertical channel row: */
- for(i=1;i<nch2;++i)
- {
- /* This `ahigh' is actually the alow in the FIRST LOOP. */
- ahigh = i*npch0<hk0 ? 0 : i*npch0-hk0;
-
- /* For a check:
- printf("2: alow: %-4zu ahigh: %-4zu. b: %-10zu "
- "(%zu -- %zu)\n", alow, ahigh, b, blow, bhigh);
- */
- for(a=alow;a<ahigh;++a)
- {
- while(*ind==GAL_THREADS_NON_THRD_INDEX) ++ind;
- *ind++=a*s1+b;
- }
- /* This alow is actually the ahigh in the FIRST LOOP. */
- alow = i*npch0+hk0>s0 ? s0 : i*npch0+hk0;
- }
-
- /* All the areas under were checked, now we have to add the
- area ontop of the highest middle channel edge. Note that
- when there is no channel along the first C axis
- (nch2==1), then it will immediately come here and scan
- all the rows between the blow and bhigh columns. */
- ahigh=s0;
-
- /* For a check:
- printf("3: alow: %-4zu ahigh: %-4zu. b: %-10zu (%zu -- %zu)\n",
- alow, ahigh, b, blow, bhigh);
- */
- for(a=alow;a<ahigh;++a)
- {
- while(*ind==GAL_THREADS_NON_THRD_INDEX) ++ind;
- *ind++=a*s1+b;
- }
- }
- }
-}
-
-
-
-
-
-/* Convolution is a very expensive (time consuming) operation. The
- problem is that sometimes, convolution was done on each channel
- independently, but later on, a program might need convolution over
- the full image (for example, during the program the differences
- between the channels has been calculated and removed), so it is
- very important to remove the discontinuities that the initial
- convolution on each channel can cause.
-
- This is the job of this function. It recieves a mesh structure and
- an already convolved image. Only the pixels lying within half of
- the PSF width of channel borders are chosen and only they are
- convolved (this time over the full image). Most of the image pixels
- (whose distance from the channel edges is more than half the PSF),
- do not need to undergo convolution again.
-
- Note that the pixels on the edges of the image do not need to undergo
- this correction. Basically this function is very similar to
- gal_spatialconvolve_convolve (spatialconvolve.c), other than the fact
- that the indexs are not over the full image but only a select number of
- pixels. */
-void
-gal_mesh_change_to_full_convolution(struct gal_mesh_params *mp, float *conv)
-{
- int err;
- pthread_t t; /* All thread ids saved in this, not used. */
- pthread_attr_t attr;
- pthread_barrier_t b;
- struct gal_spatialconvolve_params *scp;
- size_t i, nb, *indexs, numpix, thrdcols;
-
- /* If convolution was done over the full image, then there is
- nothing this function should do so just return. After this
- function, the image is fully convolved, so fullconvolution should
- be set to 1. */
- if(mp->nch==1 || mp->fullconvolution) return;
- mp->fullconvolution=1;
-
-
- /* Array keeping thread parameters for each thread.*/
- errno=0;
- scp=malloc(mp->numthreads*sizeof *scp);
- if(scp==NULL)
- error(EXIT_FAILURE, errno,
- "%zu bytes for scp in gal_mesh_change_to_full_convolution "
- "(mesh.c)", mp->numthreads*sizeof *scp);
-
-
- /* Put the indexs of the pixels to re-convolve here. */
- corrconvindexs(mp, &indexs, &numpix, &thrdcols);
-
-
- /* Start the convolution on the desired pixels. */
- if(mp->numthreads==1)
- {
- gal_spatialconvolve_pparams(mp->img, mp->s0, mp->s1, mp->kernel,
- mp->ks0, mp->ks1, mp->numthreads, 1,
- conv, indexs, &scp[0]);
- gal_spatialconvolve_thread(&scp[0]);
- }
- else
- {
- /* Initialize the attributes. Note that this running thread
- (that spinns off the nt threads) is also a thread, so the
- number the barrier should be one more than the number of
- threads spinned off. */
- if(numpix<mp->numthreads) nb=numpix+1;
- else nb=mp->numthreads+1;
- gal_threads_attr_barrier_init(&attr, &b, nb);
-
- /* Spin off the threads: */
- for(i=0;i<mp->numthreads;++i)
- if(indexs[i*thrdcols]!=GAL_THREADS_NON_THRD_INDEX)
- {
- scp[i].b=&b;
- gal_spatialconvolve_pparams(mp->img, mp->s0, mp->s1, mp->kernel,
- mp->ks0, mp->ks1, mp->numthreads, 1,
- conv, &indexs[i*thrdcols], &scp[i]);
- err=pthread_create(&t, &attr, gal_spatialconvolve_thread,
- &scp[i]);
- if(err)
- error(EXIT_FAILURE, 0, "can't create thread %zu", i);
- }
-
- /* Wait for all threads to finish and free the spaces. */
- pthread_barrier_wait(&b);
- pthread_attr_destroy(&attr);
- pthread_barrier_destroy(&b);
- }
-
- free(scp);
- free(indexs);
-}
diff --git a/lib/multidim.c b/lib/multidim.c
new file mode 100644
index 0000000..837f1f4
--- /dev/null
+++ b/lib/multidim.c
@@ -0,0 +1,170 @@
+/*********************************************************************
+multidim -- Functions for multi-dimensional operations.
+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/>.
+**********************************************************************/
+#include <config.h>
+
+#include <stdio.h>
+#include <errno.h>
+#include <error.h>
+#include <stdlib.h>
+
+#include <gnuastro/multidim.h>
+
+
+
+
+
+/************************************************************************/
+/******************** Info **********************/
+/************************************************************************/
+size_t
+gal_multidim_total_size(size_t ndim, size_t *dsize)
+{
+ size_t i, num=1;
+ for(i=0;i<ndim;++i) num *= dsize[i];
+ return num;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/************************************************************************/
+/******************** Coordinates **********************/
+/************************************************************************/
+/* Return the index of an element from its coordinates. The index is the
+ position in the contiguous array (assuming it is a 1D arrray). */
+size_t
+gal_multidim_coord_to_index(size_t ndim, size_t *dsize, size_t *coord)
+{
+ size_t i, d, ind=0, in_all_faster_dim;
+
+ switch(ndim)
+ {
+ case 0:
+ error(EXIT_FAILURE, 0, "`gal_multidim_coord_to_index' doesn't accept "
+ "0 dimensional arrays");
+
+ case 1:
+ ind=coord[0];
+ break;
+
+ case 2:
+ ind=coord[0]*dsize[1]+coord[1];
+ break;
+
+ default:
+ for(d=0;d<ndim;++d)
+ {
+ /* First, find the number of elements in all dimensions faster
+ than this one. */
+ in_all_faster_dim=1;
+ for(i=d+1;i<ndim;++i)
+ in_all_faster_dim *= dsize[i];
+
+ /* Multiply it by the coordinate value of this dimension and add
+ to the index. */
+ ind += coord[d] * in_all_faster_dim;
+ }
+ }
+
+ /* Return the derived index. */
+ return ind;
+}
+
+
+
+
+
+/* You know the index (`ind') of a point/tile in an n-dimensional (`ndim')
+ array which has `dsize[i]' elements along dimension `i'. You want to
+ know the coordinates of that point along each dimension. The output is
+ not actually returned, it must be allocated (`ndim' elements) before
+ calling this function. This function will just fill it. The reason for
+ this is that this function will often be called with a loop and a single
+ allocated space would be enough for the whole loop. */
+void
+gal_multidim_index_to_coord(size_t ind, size_t ndim, size_t *dsize,
+ size_t *coord)
+{
+ size_t d, i, in_all_faster_dim;
+
+ switch(ndim)
+ {
+ case 0:
+ error(EXIT_FAILURE, 0, "a 0-dimensional dataset is not defined in "
+ "`gal_multidim_ind_to_coord'");
+
+ /* One dimensional dataset. */
+ case 1:
+ coord[0] = ind;
+ break;
+
+ /* 2D dataset. */
+ case 2:
+ coord[0] = ind / dsize[1];
+ coord[1] = ind % dsize[1];
+ break;
+
+ /* Higher dimensional datasets. */
+ default:
+ /* We start with the slowest dimension (first in the C standard). */
+ for(d=0;d<ndim;++d)
+ {
+ /* First, find the number of elements in all dimensions faster
+ than this one. */
+ in_all_faster_dim=1;
+ for(i=d+1;i<ndim;++i)
+ in_all_faster_dim *= dsize[i];
+
+ /* If we are on the fastest dimension (last in the C standard,
+ just before the loop finishes), then no division must be
+ done. */
+ if(d==ndim-1)
+ coord[d]=ind;
+ else
+ {
+ /* Set the coordinate value for this dimension. */
+ coord[d] = ind / in_all_faster_dim;
+
+ /* Replace the index with its remainder with the number of
+ elements in all faster dimensions. */
+ ind %= in_all_faster_dim;
+ }
+ }
+ }
+}
diff --git a/lib/options.c b/lib/options.c
index 1953a3b..01d5393 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -298,7 +298,7 @@ gal_options_parse_list_of_numbers(char *string, char
*filename, size_t lineno)
/**********************************************************************/
-/************ Convert values ***************/
+/************ Parser functions for common options ***************/
/**********************************************************************/
void *
gal_options_read_type(struct argp_option *option, char *arg,
@@ -415,6 +415,92 @@ gal_options_read_tableformat(struct argp_option *option,
char *arg,
+/* Parse the given string into a series of size values (integers, stored as
+ an array of size_t). The ouput array will be stored in the `value'
+ element of the option. The last element of the array is `-1' to allow
+ finding the number of elements within it later (similar to a string
+ which terminates with a '\0' element). */
+#define PARSE_SIZES_STATICSTR_LEN 2000
+void *
+gal_options_parse_sizes_reverse(struct argp_option *option, char *arg,
+ char *filename, size_t lineno, void *params)
+{
+ int i;
+ double *v;
+ gal_data_t *values;
+ size_t nc, num, *array;
+ char *str, sstr[PARSE_SIZES_STATICSTR_LEN];
+
+ /* We want to print the stored values. */
+ if(lineno==-1)
+ {
+ /* Find the number of elements within the array. */
+ array = *(size_t **)(option->value);
+ for(i=0; array[i]!=-1; ++i);
+ num=i;
+
+ /* Write all the dimensions into the static string. */
+ nc=0;
+ for(i=num-1;i>=0;--i)
+ {
+ if( nc > PARSE_SIZES_STATICSTR_LEN-100 )
+ error(EXIT_FAILURE, 0, "a bug! please contact us at %s so we "
+ "can address the problem. The number of necessary "
+ "characters in the statically allocated string of "
+ "`gal_options_parse_sizes' has become too close to %d",
+ PACKAGE_BUGREPORT, PARSE_SIZES_STATICSTR_LEN);
+ nc += sprintf(sstr+nc, "%zu,", array[i]);
+ }
+ sstr[nc-1]='\0';
+
+ /* Copy the string into a dynamically allocated space, because it
+ will be freed later.*/
+ gal_checkset_allocate_copy(sstr, &str);
+ return str;
+ }
+
+ /* We want to read the user's string. */
+ else
+ {
+ /* If the option is already set, just return. */
+ if(option->set) return NULL;
+
+ /* Read the values. */
+ values=gal_options_parse_list_of_numbers(arg, filename, lineno);
+
+ /* Check if the values are an integer. */
+ v=values->array;
+ for(i=0;i<values->size;++i)
+ {
+ if(v[i]<0)
+ error(EXIT_FAILURE, 0, "the given value in `%s' (%g) is not 0 "
+ "or positive. The values to the `--%s' option must be "
+ "positive", arg, v[i], option->name);
+
+ if(ceil(v[i]) != v[i])
+ error(EXIT_FAILURE, 0, "the given value in `%s' (%g) is not an "
+ "integer. The values to the `--%s' option must be "
+ "integers", arg, v[i], option->name);
+ }
+
+ /* Write the values into an allocated size_t array and finish it with
+ a `-1' so the total number can be found later.*/
+ num=values->size;
+ array=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, num+1);
+ for(i=0;i<num;++i) array[num-1-i]=v[i];
+ array[num] = (size_t)(-1);
+
+ /* Put the array of size_t into the option, clean up and return.*/
+ *(size_t **)(option->value) = array;
+ gal_data_free(values);
+ return NULL;
+ }
+}
+
+
+
+
+
diff --git a/lib/options.h b/lib/options.h
index d19b401..1685052 100644
--- a/lib/options.h
+++ b/lib/options.h
@@ -232,6 +232,10 @@ void *
gal_options_read_tableformat(struct argp_option *option, char *arg,
char *filename, size_t lineno, void *junk);
+void *
+gal_options_parse_sizes_reverse(struct argp_option *option, char *arg,
+ char *filename, size_t lineno, void *params);
+
@@ -247,6 +251,7 @@ gal_options_common_argp_parse(int key, char *arg, struct
argp_state *state);
+
/**********************************************************************/
/************ Configuration files ***************/
/**********************************************************************/
diff --git a/lib/spatialconvolve.c b/lib/spatialconvolve.c
deleted file mode 100644
index 10426a6..0000000
--- a/lib/spatialconvolve.c
+++ /dev/null
@@ -1,237 +0,0 @@
-/*********************************************************************
-SpatialConvolve - Convolve an image in the spatial domain.
-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.
-
-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/>.
-**********************************************************************/
-#include <config.h>
-
-#include <math.h>
-#include <stdio.h>
-#include <errno.h>
-#include <error.h>
-#include <string.h>
-#include <stdlib.h>
-
-#include <gnuastro/box.h>
-#include <gnuastro/spatialconvolve.h>
-
-
-
-
-
-
-
-
-
-
-/*******************************************************************/
-/************** Each thread ****************/
-/*******************************************************************/
-void
-gal_spatialconvolve_pparams(float *input, size_t is0, size_t is1,
- float *kernel, size_t ks0, size_t ks1,
- size_t nt, int edgecorrection, float *out,
- size_t *indexs,
- struct gal_spatialconvolve_params *scp)
-{
- /* Put the simple values in: */
- scp->is0=is0;
- scp->is1=is1;
- scp->ks0=ks0;
- scp->ks1=ks1;
- scp->out=out;
- scp->input=input;
- scp->kernel=kernel;
- scp->indexs=indexs;
- scp->numthreads=nt;
- scp->edgecorrection=edgecorrection;
-}
-
-
-
-
-
-
-void *
-gal_spatialconvolve_thread(void *inparam)
-{
- struct gal_spatialconvolve_params *scp=
- (struct gal_spatialconvolve_params *)inparam;
-
- double sum, ksum;
- long naxes[2]={scp->is1, scp->is0};
- float *f, *fp, *k, *istart, *kstart;
- int edgecorrection=scp->edgecorrection;
- size_t is1=scp->is1, ks0=scp->ks0, ks1=scp->ks1;
- long *fpixel_i=scp->fpixel_i, *lpixel_i=scp->lpixel_i;
- long *fpixel_o=scp->fpixel_o, *lpixel_o=scp->lpixel_o;
- size_t i, j, ind, *indexs=scp->indexs, numrows, numcols;
- float *input=scp->input, *kernel=scp->kernel, *out=scp->out;
-
- /* Go over all the pixels associated with this thread. */
- for(i=0;indexs[i]!=GAL_THREADS_NON_THRD_INDEX;++i)
- {
- /* Set the index, if it is a NaN pixel, then set the output to
- be NaN too. */
- ind=indexs[i];
- if(isnan(input[ind]))
- {
- out[ind]=NAN;
- continue;
- }
-
- /* Set the starting and ending pixels on the kernel, note that
- the overlap function in box.c uses FITS coordinates. */
- fpixel_o[0]=1; fpixel_o[1]=1;
- lpixel_o[0]=ks1; lpixel_o[1]=ks0;
- fpixel_i[0]=ind%is1+1-ks1/2; fpixel_i[1]=ind/is1+1-ks0/2;
- lpixel_i[0]=ind%is1+1+ks1/2; lpixel_i[1]=ind/is1+1+ks0/2;
- gal_box_overlap(naxes, fpixel_i, lpixel_i, fpixel_o, lpixel_o);
-
- /* fpixels and lpixels now point to the overlap's starting and
- ending both on the image and on the kernel. */
- sum=0.0f;
- numcols=lpixel_i[0]-fpixel_i[0]+1; /* lpixel is inside the box. */
- numrows=lpixel_i[1]-fpixel_i[1]+1; /* lpixel is inside the box. */
- ksum = edgecorrection ? 0.0f : 1.0f;
- istart = &input[ (fpixel_i[1]-1) * is1 + fpixel_i[0]-1 ];
- kstart = &kernel[ (fpixel_o[1]-1) * ks1 + fpixel_o[0]-1 ];
- for(j=0;j<numrows;++j)
- {
- k = kstart + j*ks1;
- fp = (f = istart + j*is1) + numcols;
- do
- {
- if( isnan(*f)==0 )
- {
- sum += *k * *f;
- if(edgecorrection) ksum+=*k;
- }
- ++k;
- }
- while(++f<fp);
- out[ind]=sum/ksum;
- }
- }
-
- /* Wait until all other threads finish. */
- if(scp->numthreads>1)
- pthread_barrier_wait(scp->b);
-
- return NULL;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/*******************************************************************/
-/************** Outside function ****************/
-/*******************************************************************/
-void
-gal_spatialconvolve_convolve(float *input, size_t is0, size_t is1,
- float *kernel, size_t ks0, size_t ks1,
- size_t numthreads, int edgecorrection,
- float **out)
-{
- int err;
- pthread_t t; /* All thread ids saved in this, not used. */
- pthread_attr_t attr;
- pthread_barrier_t b;
- struct gal_spatialconvolve_params *scp;
- size_t i, nb, *indexs, thrdcols;
-
-
- /* Array keeping thread parameters for each thread. */
- errno=0;
- scp=malloc(numthreads*sizeof *scp);
- if(scp==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes in gal_spatialconvolve_convolve "
- "(spatialconvolve.c) for scp", numthreads*sizeof *scp);
-
-
- /* Allocate the output array: */
- errno=0;
- *out=malloc(is0*is1*sizeof **out);
- if(*out==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for convolution output",
- is0*is1*sizeof **out);
-
-
- /* Distribute the image pixels into the threads: */
- gal_threads_dist_in_threads(is0*is1, numthreads, &indexs, &thrdcols);
-
- /* Start the convolution. */
- if(numthreads==1)
- {
- gal_spatialconvolve_pparams(input, is0, is1, kernel, ks0, ks1,
- numthreads, edgecorrection, *out, indexs,
- &scp[0]);
- gal_spatialconvolve_thread(&scp[0]);
- }
- else
- {
- /* Initialize the attributes. Note that this running thread (that
- spinns off the numthreads threads) is also a thread, so the number
- the barrier should be one more than the number of threads spinned
- off. */
- if(is0*is1<numthreads) nb=is0*is1+1;
- else nb=numthreads+1;
- gal_threads_attr_barrier_init(&attr, &b, nb);
-
- /* Spin off the threads: */
- for(i=0;i<numthreads;++i)
- if(indexs[i*thrdcols]!=GAL_THREADS_NON_THRD_INDEX)
- {
- scp[i].b=&b;
- gal_spatialconvolve_pparams(input, is0, is1, kernel, ks0,
- ks1, numthreads, edgecorrection,
- *out, &indexs[i*thrdcols], &scp[i]);
- err=pthread_create(&t, &attr, gal_spatialconvolve_thread,
- &scp[i]);
- if(err)
- error(EXIT_FAILURE, 0, "can't create thread %zu", i);
- }
-
- /* Wait for all threads to finish and free the spaces. */
- pthread_barrier_wait(&b);
- pthread_attr_destroy(&attr);
- pthread_barrier_destroy(&b);
- }
-
- /* Free the allocated spaces: */
- free(scp);
- free(indexs);
-}
diff --git a/lib/tile.c b/lib/tile.c
new file mode 100644
index 0000000..7bd9bb1
--- /dev/null
+++ b/lib/tile.c
@@ -0,0 +1,471 @@
+/*********************************************************************
+tile -- work with tesselations over a host dataset.
+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/>.
+**********************************************************************/
+#include <config.h>
+
+#include <stdio.h>
+#include <errno.h>
+#include <error.h>
+#include <stdlib.h>
+
+#include <gnuastro/tile.h>
+#include <gnuastro/multidim.h>
+
+
+
+
+
+
+/***********************************************************************/
+/************** Allocated block of memory ******************/
+/***********************************************************************/
+/* When you are working on an array, it important to know the size of the
+ allocated space in each dimension. This simple function will just follow
+ the block pointer and return the `dsize' element of lowest-level
+ structure. */
+gal_data_t *
+gal_tile_block(gal_data_t *input)
+{
+ while(input->block!=NULL) input=input->block;
+ return input;
+}
+
+
+
+
+
+/* Calculate the starting coordinates of a tile in the allocated block of
+ memory. */
+void
+gal_tile_block_tile_start_coord(gal_data_t *tile, size_t *start_coord)
+{
+ size_t *s, *sf, ind;
+ gal_data_t *block=gal_tile_block(tile);
+
+ /* If the input tile is actually the same as the block, then the
+ reference is all zeros. */
+ if(block==tile)
+ {
+ sf = (s=start_coord) + tile->ndim;
+ do *s++=0; while(s<sf);
+ return;
+ }
+
+ /* Calculate the coordinates of the first pixel of the tile. */
+ ind = tile->array - block->array;
+ gal_multidim_index_to_coord(ind, tile->ndim, block->dsize, start_coord);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/***********************************************************************/
+/************** Tile full dataset ********************/
+/***********************************************************************/
+/* This option is mainly used by Gnuastro's programs to see if the user's
+ given values for the tile size (`tile') and number of channels
+ (`channels') corresponds to the input dataset. */
+size_t *
+gal_tile_all_sanity_check(char *filename, char *hdu, gal_data_t *input,
+ size_t *tile, size_t *numchannels)
+{
+ double d;
+ size_t i, *out;
+
+ /* Check the tile's dimensions. */
+ for(i=0;tile[i]!=-1;++i)
+ {
+ /* Not equal to zero. */
+ if(tile[i]==0)
+ error(EXIT_FAILURE, 0, "the tile size must be larger than zero");
+
+ /* If the tile size is larger than the dataset size in this
+ dimension, then change the tile size to the dataset size. */
+ if(tile[i]>input->dsize[i]) tile[i]=input->dsize[i];
+ }
+
+
+ /* Make sure the number of tile sizes (tile dimensions) are the same as
+ the dataset's dimensions). */
+ if(i!=input->ndim)
+ error(EXIT_FAILURE, 0, "%s (hdu: %s): has %zu dimensions, but only %zu "
+ "value(s) given for the tile size", filename, hdu, input->ndim, i);
+
+
+ /* Check the channel's dimensions. */
+ for(i=0;numchannels[i]!=-1;++i)
+ if(tile[i]==0)
+ error(EXIT_FAILURE, 0, "the number of channels must be larger than "
+ "zero");
+ if(i!=input->ndim)
+ error(EXIT_FAILURE, 0, "%s (hdu: %s): has %zu dimensions, but only %zu "
+ "value(s) given for the number of channels", filename, hdu,
+ input->ndim, i);
+
+
+ /* Allocate space for the channel sizes. */
+ out=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim);
+
+
+ /* Check if the channels are exactly divisible by the input's size along
+ each dimension and set the correct size. */
+ for(i=0;i<input->ndim;++i)
+ {
+ /* Check if the number of channels is not more than the size of the
+ image. Note that the reported dimension must be in FITS format.*/
+ if(input->dsize[i]<numchannels[i])
+ error(EXIT_FAILURE, 0, "the number of channels in dimension %zu "
+ "(%zu) is more than the size of the `%s' (hdu: %s) in that "
+ "dimension", input->ndim-i, numchannels[i], filename, hdu);
+
+ /* Also check the tile size. */
+ if(input->dsize[i]<tile[i])
+ error(EXIT_FAILURE, 0, "the tile size in dimension %zu (%zu) is "
+ "more than the size of the `%s' (hdu: %su) in that dimension",
+ input->ndim-i, tile[i], filename, hdu);
+
+ /* First check. */
+ d=(double)input->dsize[i]/(double)numchannels[i];
+ if(ceil(d)!=d)
+ error(EXIT_FAILURE, 0, "%zu (number of channels along dimension "
+ "%zu) is not exactly divisible by %zu (the length of `%s' "
+ "(hdu: %s) that dimension). The channels cover the input "
+ "dataset, hence, they must be identical", numchannels[i],
+ input->ndim-i, input->dsize[i], filename, hdu);
+
+ /* Put the channel size into the output. */
+ out[i]=d;
+ }
+
+ /* Return the output array (size of channel along each dimension). */
+ return out;
+}
+
+
+
+
+
+/* The user's specified tile size might not be an exact multiple of the
+ parent's size. This function is useful in such cases. It will give the
+ starting tile's size along each dimension. The line below can be the
+ length along any dimension and the tile size along that dimension. You
+ see that when we start with the tile size, we will end up with a last
+ tile that contains the remainder elements.
+
+ | tile size | tile size | tile size | tile size | remainder
+ | | | | | |
+ ---------------------------------------------------------
+
+ The remainder will always be smaller than `tile size'. So, we will merge
+ the last tile size with the remainder and move that tile to the start.
+ In this way, the size of the first tile will always be between between
+ one and two times the size of the regular tile:
+
+ | first tile | tile size | tile size | tile size |
+ | | | | |
+ ---------------------------------------------------------
+
+ When there is only a small remainder (for example one or two pixels),
+ then this layout is fine. But when the remainder is significant compared
+ to the regular tile size (like the example above), then it will make
+ more sense to cut the first tile into two halfs (`f-half' and `l-half')
+ and put them at the start and end of the full length:
+
+
+ | f-half | tile size | tile size | tile size | l-half |
+ | | | | | |
+ ---------------------------------------------------------
+
+ So in any case, knowing the size of the first tile, will allow us to
+ parse all the tiles. We just have to make sure we don't go over the full
+ input's length. */
+static void
+gal_tile_all_regular_first(gal_data_t *parent, size_t *regular,
+ float significance, size_t *first, size_t *last,
+ size_t *number)
+{
+ size_t i, remainder, *dsize=parent->dsize;;
+
+ /* For each dimension, set the size of the first tile. */
+ for(i=0;i<parent->ndim;++i)
+ {
+ /* Calculate the remainder in this dimension. */
+ remainder=dsize[i] % regular[i];
+
+ /* Depending on the remainder, set the first tile size and number. */
+ if(remainder)
+ {
+ if( remainder > significance * regular[i] )
+ {
+ first[i] = ( remainder + regular[i] )/2;
+ number[i] = dsize[i]/regular[i] + 1 ;
+ last[i] = dsize[i] - ( first[i] + regular[i]*(number[i]-2) );
+ }
+ else
+ {
+ first[i] = remainder + regular[i];
+ number[i] = dsize[i]/regular[i];
+ last[i] = regular[i];
+ }
+ }
+ else
+ {
+ first[i] = last[i] = regular[i];
+ number[i] = dsize[i]/regular[i];
+ }
+ }
+}
+
+
+
+
+
+/* Cover the full dataset with (mostly) identical tiles. The regular tile
+ size is determined from the `size' array. If the input data's size is
+ not an exact multiple of `size' 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. For a full description of tiling in
+ `gal_data_t', please see `data.h'.
+
+ Inputs
+ ------
+
+ `input' is the gal_data_t which you want to tile (only used for its
+ sizes).
+
+ `regular' is 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 `input'.
+
+ `out' is the pointer to the array of data structures that is to keep
+ the tile parameters. If `*out==NULL', then the necessary space will
+ be allocated. If it is not NULL, then all the tile information will
+ be filled from the given element, see `multiple' for more.
+
+ `multiple': When the `*out' array is to be allocated, allocate
+ `multiple' times the necessary space. This can be very useful when
+ you have several more identically sized `inputs', and you want all
+ their tiles to be allocated (and thus indexed) together, even
+ though they have different `block' datasets (that then link to one
+ allocated space). See the `gal_tile_all_position_channel_tile'
+ below.
+
+ Output
+ ------
+
+ The returned output is an array of `gal_data_t' with `numtiles'
+ elements, note that you have to pass the pointer to `numtiles' as one
+ of the arguments.
+
+
+ Implementation
+ --------------
+
+ In the most general case, to set the starting pointers for each tile
+ we need the following sizes. If the input array has no parent/block,
+ then both these sizes are equal to it's own size:
+
+ 1. block-size (or `bsize'), which is the size of the allocated
+ 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).
+*/
+size_t
+gal_tile_all_position(gal_data_t *input, size_t *regular, gal_data_t **out,
+ size_t multiple)
+{
+ size_t i, d, tind, numtiles, *start=NULL;
+ gal_data_t *tiles, *block=gal_tile_block(input);
+ size_t *last = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim);
+ size_t *tsize = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim);
+ size_t *first = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim);
+ size_t *coord = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim);
+ size_t *tcoord = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim);
+
+
+ /* Set the first tile size and total number of tiles along each
+ dimension, then allocate the array of tiles. */
+ gal_tile_all_regular_first(input, regular, 0.3, first, last, tsize);
+ numtiles=gal_multidim_total_size(input->ndim, tsize);
+
+
+ /* Allocate the necessary space for all the tiles (if necessary). */
+ if(*out) tiles = *out;
+ else *out = tiles = gal_data_array_calloc(numtiles*multiple);
+
+
+ /* It is possible that the `input' dataset is its-self a larger tile over
+ a region of the allocated block. In that case, we need to account for
+ the block's dimensions when calculating the position of this block. */
+ if(input->block)
+ {
+ start=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim);
+ gal_tile_block_tile_start_coord(input, start);
+ }
+
+
+ /* Initialize each tile. */
+ for(i=0;i<numtiles;++i)
+ {
+ /* Specify the coordinates of the tile between the other tiles. Note
+ that we are dealing with tiles here, not pixels. */
+ gal_multidim_index_to_coord(i, input->ndim, tsize, tcoord);
+
+ /* The coordinates are currently in units of tiles, not
+ pixels. Convert them to the coordinates of the first pixel in each
+ tile. */
+ for(d=0;d<input->ndim;++d)
+ {
+ /* Convert the tile coordinates to pixel coordinates within
+ `input'. See the comments above `gal_tile_all_regular_first':
+ The first tile in every dimension can be different from the
+ regular tile size. */
+ coord[d] = tcoord[d] ? first[d] + (tcoord[d]-1)*regular[d] : 0;
+
+ /* When the `input' data structure (that is to be tiled here) was
+ itself a tile over a larger allocated array, a `start' array
+ has been allocated to correct the coordinates so they refer to
+ a physical position on the allocated block of memory. */
+ if(start)
+ coord[d] += start[d];
+ }
+
+ /* Convert the coordinates (that are now in element/pixel units on
+ the allocated block of memory) into an index. */
+ tind=gal_multidim_coord_to_index(block->ndim, block->dsize, coord);
+
+ /* Now that we have the index of this tile's starting point compared
+ to the allocated block, put it in to the tile's `array'
+ pointer. */
+ tiles[i].array=block->array+tind;
+
+ /* Set the sizes of the tile. */
+ tiles[i].ndim=input->ndim;
+ tiles[i].dsize=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T,input->ndim);
+ for(d=0;d<input->ndim;++d)
+ {
+ /* The size of the first and last tiles can be different from the
+ majority of the `regular' tiles that have the same size. When
+ a tile is on the edge in one of the dimensions, then its
+ `coord[d]' will be either 0 or the last. */
+ if( first[d] != regular[d]
+ && ( tcoord[d]==0 || tcoord[d]==tsize[d]-1 ) )
+ {
+ if( tcoord[d] == 0 ) tiles[i].dsize[d] = first[d];
+ if( tcoord[d] == tsize[d]-1 ) tiles[i].dsize[d] = last[d];
+ }
+ else
+ tiles[i].dsize[d]=regular[d];
+ }
+
+ /* Set the block structure for this tile to the `input'. */
+ tiles[i].block=input;
+
+ /* For a check:
+ printf("%zu:\n\tStart index: %zu\n\tsize: %zu x %zu\n", i, tind,
+ tiles[i].dsize[1], tiles[i].dsize[0]);
+ exit(0);
+ */
+ }
+
+
+ /* Clean up and return. */
+ free(last);
+ free(tsize);
+ free(first);
+ free(coord);
+ free(tcoord);
+ if(start) free(start);
+ return numtiles;
+}
+
+
+
+
+
+/* A dataset can be tiled with two layers that are related:
+
+ Channels: A tesselation of larger tile sizes that all have the same
+ size (`channel_size' must be an exact multiple of `input's size
+ along every dimension. In astronomy images, this can be seen as
+ CCD amplifiers, that cover large parts of the image. If
+ `*channels!=NULL' then it is assumed to be already present and
+ will not be allocated.
+
+ Tiles: A combined tesselation of each channel with smaller
+ tiles. These tiles can be used to calculate things like
+ gradients over each channel and thus over the whole image. */
+void
+gal_tile_all_position_two_layers(gal_data_t *input, size_t *channel_size,
+ size_t *tile_size, gal_data_t **channels,
+ gal_data_t **tiles, size_t *numchannels,
+ size_t *numtiles)
+{
+ gal_data_t *ch, *t;
+ size_t i, nch, ntiles_in_ch;
+
+ /* First allocate the channels. Note that the channels tesselation might
+ have already been set. */
+ if(*channels)
+ {
+ nch=1;
+ for(i=0;i<input->ndim;++i)
+ nch *= input->dsize[i]/channel_size[i];
+ }
+ else
+ /* Note that the actual allocated input array will be the direct
+ `block' of each channel. */
+ nch=gal_tile_all_position(input, channel_size, channels, 1);
+
+
+ /* Now, tile each channel. While tiling the first channel, we are also
+ going to allocate the space for the other channels. Then pass those
+ pointers. */
+ *tiles=NULL;
+ ch=*channels;
+ ntiles_in_ch = gal_tile_all_position(ch, tile_size, tiles, nch);
+ for(i=1;i<nch;++i)
+ {
+ t = *tiles + i*ntiles_in_ch;
+ gal_tile_all_position(&ch[i], tile_size, &t, 1);
+ }
+
+ /* Return the total number of channels and tiles */
+ *numchannels = nch;
+ *numtiles = nch * ntiles_in_ch;
+}
diff --git a/lib/txt.c b/lib/txt.c
index 450fb46..1d8c5a6 100644
--- a/lib/txt.c
+++ b/lib/txt.c
@@ -458,7 +458,7 @@ txt_infoll_to_array(gal_data_t *datall, size_t *numdata)
if(numc>1)
{
/* Now, allocate the array and put in the values. */
- dataarr=gal_data_calloc_dataarray(numc);
+ dataarr=gal_data_array_calloc(numc);
/* Put each dataset/column into its proper place in the array. */
while(datall!=NULL)
- [gnuastro-commits] master 0bc96c9 122/125: Clumps and their S/N tables over detections, (continued)
- [gnuastro-commits] master 0bc96c9 122/125: Clumps and their S/N tables over detections, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 80cdd4f 118/125: Correction of convolution over channel borders, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e90b490 116/125: Internal headers moved to special directory, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 29d7e46 091/125: Statistics now successfully uses gal_data_t, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 60c126f 120/125: New implementation of over-segmentation to find clumps, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master c250a6d 123/125: NoiseChisel's clumps grown and objects defined, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master ffaa059 101/125: Tessellation options now common to all programs, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master dec2b22 078/125: ImageCrop works with gal_data_t features, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master f8b1820 061/125: Other configuration files read immediately when seen, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 46b4651 085/125: Image prefix removed from some program names, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master df716bc 097/125: Tesselation with the new data structure: work started,
Mohammad Akhlaghi <=
- [gnuastro-commits] master db39cee 088/125: Statistics returns basic results, no array lib, new blank lib, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 60cd76e 090/125: Statistics mostly complete, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 15c524d 124/125: New implementation of NoiseChisel complete, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 3b26d5c 114/125: SubtractSky program removed, Statistics estimates Sky, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 99ac7dc 083/125: Gnuastro now uses fixed width (from stdint.h) types, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master eaf19a3 115/125: NoiseChisel's initial detection step implemented, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 3339b92 125/125: New implementation of MakeCatalog now complete, Mohammad Akhlaghi, 2017/04/23