[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 95ba655 100/125: Spatial convolution with new
From: |
Mohammad Akhlaghi |
Subject: |
[gnuastro-commits] master 95ba655 100/125: Spatial convolution with new tessellation |
Date: |
Sun, 23 Apr 2017 22:36:47 -0400 (EDT) |
branch: master
commit 95ba655a24257c66adb32124f21e1d2fa92d6146
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Spatial convolution with new tessellation
Spatial convolution has been fully implemented using the new tessellation
system defined in `lib/tile.c'. It can also account for individual
convolution on each channel as before, but not using any extra structure!
It was fully done with `gal_data_t'.
---
bin/convolve/args.h | 12 +--
bin/convolve/astconvolve.conf | 2 +-
bin/convolve/convolve.c | 5 +-
bin/convolve/main.h | 4 +-
bin/convolve/ui.c | 14 +--
bin/convolve/ui.h | 2 +-
doc/gnuastro.texi | 157 ++++++++++++++-------------
lib/convolve.c | 243 +++++++++++++++++++++++++++++++++---------
lib/tile.c | 38 +++++--
9 files changed, 323 insertions(+), 154 deletions(-)
diff --git a/bin/convolve/args.h b/bin/convolve/args.h
index ab310b4..98bd78e 100644
--- a/bin/convolve/args.h
+++ b/bin/convolve/args.h
@@ -122,17 +122,17 @@ struct argp_option program_options[] =
/* Tile grid options. */
{
0, 0, 0, 0,
- "Tile grid (only for spatial domain):",
+ "Tessellation (only for spatial domain):",
ARGS_GROUP_MESH_GRID
},
{
- "tile",
- ARGS_OPTION_KEY_TILE,
+ "tilesize",
+ ARGS_OPTION_KEY_TILESIZE,
"INT[,INT]",
0,
- "Size of tiles along each dim. (FITS order).",
+ "Size of regular tiles on each dim. (FITS order).",
ARGS_GROUP_MESH_GRID,
- &p->tile,
+ &p->tilesize,
GAL_DATA_TYPE_SIZE_T,
GAL_OPTIONS_RANGE_GT_0,
GAL_OPTIONS_NOT_MANDATORY,
@@ -158,7 +158,7 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_REMAINDERFRAC,
"FLT",
0,
- "Fraction of remainers in each dim to chop.",
+ "Fraction of remainder to split last tile.",
ARGS_GROUP_MESH_GRID,
&p->remainderfrac,
GAL_DATA_TYPE_FLOAT32,
diff --git a/bin/convolve/astconvolve.conf b/bin/convolve/astconvolve.conf
index 941e714..da3d0f5 100644
--- a/bin/convolve/astconvolve.conf
+++ b/bin/convolve/astconvolve.conf
@@ -25,7 +25,7 @@
type float32
# Mesh grid:
- tile 30,30
+ tilesize 30,30
numchannels 1,1
remainderfrac 0.1
convoverch 0
diff --git a/bin/convolve/convolve.c b/bin/convolve/convolve.c
index ff523d1..f179440 100644
--- a/bin/convolve/convolve.c
+++ b/bin/convolve/convolve.c
@@ -261,7 +261,7 @@ removepaddingcorrectroundoff(struct convolveparams *p)
size_t ps1=p->ps1;
size_t i, hi0, hi1;
size_t *isize=p->input->dsize;
- double *o, *input=p->input->array;
+ float *o, *input=p->input->array;
double *d, *df, *start, *pimg=p->pimg;
/* Set all the necessary parameters to crop the desired region. hi0
@@ -776,14 +776,13 @@ convolve(struct convolveparams *p)
if(p->domain==CONVOLVE_DOMAIN_SPATIAL)
{
/* Prepare the mesh structure. */
- gal_tile_all_position_two_layers(p->input, p->channel, p->tile,
+ gal_tile_all_position_two_layers(p->input, p->channelsize, p->tilesize,
p->remainderfrac, &channels, &tiles);
/* Save the tile IDs if they are requested. */
if(p->tilesname)
gal_tile_block_check_tiles(tiles, p->tilesname, PROGRAM_NAME);
-
/* Do the spatial convolution. */
out=gal_convolve_spatial(tiles, p->kernel, p->cp.numthreads,
p->convoverch);
diff --git a/bin/convolve/main.h b/bin/convolve/main.h
index ab8d157..4d08549 100644
--- a/bin/convolve/main.h
+++ b/bin/convolve/main.h
@@ -79,7 +79,7 @@ struct convolveparams
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 *tilesize; /* Size of tiles along each dim. (C order).*/
size_t *numchannels; /* No. of tiles along each dim. (C order). */
float remainderfrac; /* Frac. of remainers in each dim to cut. */
uint8_t convoverch; /* Convolve over channel borders. */
@@ -88,7 +88,7 @@ struct convolveparams
size_t makekernel; /* Make a kernel to create input. */
/* Internal */
- size_t *channel; /* Size of channels along each dimension. */
+ size_t *channelsize; /* 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. */
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index 513c3ee..eb51aff 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -233,16 +233,16 @@ ui_read_check_only_options(struct convolveparams *p)
/* 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->tilesize==NULL || p->numchannels==NULL )
{
- if( p->tile==NULL && p->numchannels==NULL )
+ if( p->tilesize==NULL && p->numchannels==NULL )
error(EXIT_FAILURE, 0, "in spatial convolution, `--numchannels' "
- "and `--tile' are mandatory");
+ "and `--tilesize' 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
+ p->tilesize ? "numchannels" : "tilesize",
+ ( p->tilesize
? "number of channels along each dimension of the input"
: "size of tiles to cover the input along each "
"dimension" ) );
@@ -360,8 +360,8 @@ ui_preparations(struct convolveparams *p)
PROGRAM_NAME);
}
else
- p->channel=gal_tile_all_sanity_check(p->filename, p->cp.hdu, p->input,
- p->tile, p->numchannels);
+ p->channelsize=gal_tile_all_sanity_check(p->filename, p->cp.hdu, p->input,
+ p->tilesize, p->numchannels);
diff --git a/bin/convolve/ui.h b/bin/convolve/ui.h
index a4df9f5..53413ad 100644
--- a/bin/convolve/ui.h
+++ b/bin/convolve/ui.h
@@ -39,7 +39,7 @@ enum option_keys_enum
ARGS_OPTION_KEY_KHDU = 'u',
ARGS_OPTION_KEY_MINSHARPSPEC = 'H',
ARGS_OPTION_KEY_CHECKFREQSTEPS = 'C',
- ARGS_OPTION_KEY_TILE = 't',
+ ARGS_OPTION_KEY_TILESIZE = 't',
ARGS_OPTION_KEY_NUMCHANNELS = 'c',
ARGS_OPTION_KEY_REMAINDERFRAC = 'r',
ARGS_OPTION_KEY_DOMAIN = 'd',
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 2623302..9d01695 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -10793,26 +10793,57 @@ are ultimately poor approximations.
@cindex Tile
@cindex Gradient
@cindex Mesh grid
-Some of the programs in Gnuastro will need to divide the pixels in an
-image into individual tiles or a mesh grid to be able to deal with
-gradients. In this section we will explain the concept in detail and
-how the user can check the grid. In the case of SubtractSky, if the
-image is completely uniform then one sky value will suffice for the
-whole image. See @ref{Sky value} for the definition of the sky value.
-Unfortunately though, as discussed in @ref{SubtractSky}, in most
-images taken with ground or space-based telescopes the sky value is
-not uniform. So we have to break the image up into small tiles on a
-mesh grid, assume the sky value is constant over them and find the sky
-value on those tiles.
-
-Meshes are considered to be a square with a side of
address@hidden pixels. The best mesh size is directly related to
-the gradient on the image. In practice we assume that the gradient is
-not significant over each mesh. So if there is a strong gradient (for
-example in long wavelength ground based images) or the image is of a
-crowded area where there isn't too much blank area, you have to choose
-a smaller mesh size. A larger mesh will give more pixels and so the
-scatter in the results will be less.
address@hidden Tile grid
address@hidden Tessellation
+Some of the programs in Gnuastro will need to divide the pixels in an image
+into a grid of individual, non-overlapping tiles. In the arts and
+mathematics, this process is formally known as
address@hidden://en.wikipedia.org/wiki/Tessellation, tessellation}. Defining a
+tessellation over a dataset (image) is necessary in many contexts. Their
+most common application is to account for local variations over the dataset
+(for example background sky gradients, which are more common as we go to
+longer/redder wavelengths/filters).
+
+In this section we will explain the concept in detail and how you can check
+the grid. In the case of SubtractSky, if the image is completely uniform
+then one sky value will suffice for the whole image. See @ref{Sky value}
+for the definition of the sky value. Unfortunately though, as discussed in
address@hidden, in most images taken with ground or space-based
+telescopes the sky value is not uniform. So we have to break the image up
+into small tiles on a grid, assume the Sky value is constant over them and
+find the sky value on those tiles.
+
+The size of the regular tiles (in units of data-elements, or pixels in an
+image) can be defined with the @option{--tilesize} option. It takes
+multiple numbers (separated by a comma) which will be the length along the
+respective dimension (in Fortran/FITS order). For example
address@hidden,40} can be used for an image (a 2D dataset). The
+regular tile size along the first FITS axis (horizontal when viewed in SAO
+ds9) will be 30 pixels and along the second it will be 40 pixels. Ideally,
address@hidden should be selected such that all tiles in the image
+have exactly the same size (in other words, that the dataset length in each
+dimension is divisible by the tile size in that dimension).
+
+However, this is not always possible: the dataset can be any size and every
+pixel in it is valuable. In such cases, Gnuastro will look at the
+significance of the remainder length, if it is not significant (for example
+one or two pixels), then it will just increase the size of the first tile
+in the respective dimension and allow the rest of the tiles to have the
+required size. When the remainder is significant (for example one pixel
+less than the size along that dimension), the remainder will be added to
+one regular tile's size and the large tile will be cut in half and put in
+the two ends of the grid/tessellation. In this way, all the tiles in the
+central regions of the dataset will have the regular tile sizes and the
+tiles on the edge will be slightly smaller. You can set the significance
+fraction with the @option{--remainderfrac} option.
+
+The best tile size is directly related to the spatial properties of the
+gradient on the image. In practice we assume that the gradient is not
+present over each tile. So if there is a strong gradient (for example in
+long wavelength ground based images) or the image is of a crowded area
+where there isn't too much blank area, you have to choose a smaller tile
+size. A larger mesh will give more pixels and so the scatter in the results
+will be less.
@cindex CCD
@cindex Amplifier
@@ -10820,69 +10851,45 @@ scatter in the results will be less.
@cindex Subaru Telescope
@cindex Hyper Suprime-Cam
@cindex Hubble Space Telescope
-For raw image processing, a simple mesh grid is not sufficient. Raw
-images are the unprocessed outputs of the camera detectors. Large
+For raw image processing, a single tesselation/grid is not sufficient. Raw
+images are the unprocessed outputs of the camera detectors. Modern
detectors usually have multiple readout channels each with its own
amplifier. For example the Hubble Space Telescope Advanced Camera for
-Surveys (ACS) has four amplifiers over its full detector area dividing
-the square field of view to four smaller squares. Ground based image
-detectors are not exempt, for example each CCD of Subaru Telescope's
-Hyper Suprime-Cam camera (which has 104 CCDs) has four amplifiers, but
-they have the same height of the CCD and divide the width by four
-parts.
+Surveys (ACS) has four amplifiers over its full detector area dividing the
+square field of view to four smaller squares. Ground based image detectors
+are not exempt, for example each CCD of Subaru Telescope's Hyper
+Suprime-Cam camera (which has 104 CCDs) has four amplifiers, but they have
+the same height of the CCD and divide the width by four parts.
@cindex Channel
-The bias current on each amplifier is different, and normally bias
-subtraction is not accurately done. So even after subtracting the measured
-bias current, you can usually still identify the boundaries of different
+The bias current on each amplifier is different, and initial bias
+subtraction is not perfect. So even after subtracting the measured bias
+current, you can usually still identify the boundaries of different
amplifiers by eye. See Figure 11(a) in Akhlaghi and Ichikawa (2015) for an
example. This results in the final reduced data to have non-uniform
amplifier-shaped regions with higher or lower background flux values. Such
systematic biases will then propagate to all subsequent measurements we do
on the data (for example photometry and subsequent stellar mass and star
-formation rate measurements in the case of galaxies). Therefore an accurate
-sky subtraction routine should also be able to account for such biases.
-
-To get an accurate result, the mesh boundaries should be located
-exactly on the amplifier boundaries. Otherwise, some meshes will
-contain pixels that have been read from two or four different
-amplifiers. These meshes are going to give very biased results and the
-amplifier boundary will still be present after sky subtraction. So we
-define `channel's. A channel is an independent mesh grid that covers
-one amplifier to ensure that the meshes do not pass the amplifier
-boundary. They can also be used in subsequent steps as the area used
-to identify nearby neighbors to interpolate and smooth the final grid,
-see @ref{Grid interpolation and smoothing}. The number of channels
-along each axis can be specified by the user at run time through the
-command-line @option{--nch1} and @option{--nch2} options or in the
-configuration files, see @ref{Configuration files}. The area of each
-channel will then be tiled by meshes of the given size and subsequent
-processing will be done on those meshes. If the image is processed or
-the detector only has one amplifier, you can set the number of
-channels in both axes to 1.
-
-Unlike the channel size, that has to be an exact multiple of the image
-size, the mesh size can be any number. If it is not an exact multiple of
-the image side, the last mesh (rightest, for the first FITS dimension, and
-highest for the second when viewed in SAO ds9) will have a different size
-than the rest. If the remainder of the image size divided by mesh size is
-larger than a certain fraction (value to @option{--lastmeshfrac}) of the
-mesh size along each axis, a new (smaller) mesh will be put there instead
-of a larger mesh. This is done to avoid the last mesh becoming too large
-compared to the other meshes in the grid. Generally, it is best practice to
-choose the mesh size such that the last mesh is only a few (negligible)
-pixels wider or thinner than the rest.
-
-The final mesh grid can be seen on the image with the
address@hidden option that is available to all programs which
-use the mesh grid for localized operations. When this option is
-called, a multi-extension FITS file with a @file{_mesh.fits} suffix
-will be created along with the outputs, see @ref{Automatic
-output}. The first extension will be the input image. For each mesh
-grid the image produces, there will be a subsequent extension. Each
-pixel in the grid extensions is labeled to the mesh that it is part
-of. You can flip through the extensions to check the mesh sizes and
-locations compared to the input image.
+formation rate measurements in the case of galaxies).
+
+Therefore an accurate analysis requires a lower-level tesselation: with
+larger tiles, each large tile covering one amplifier channel. For clarity
+we'll call these larger tiles ``channels''. The number of channels along
+each dimension is defined through the @option{--numchannels}. Each channel
+will then be covered by its own individual smaller tessellation (with tile
+sizes determined by the @option{--tilesize} option). This will allow two
+adjacent pixels from different channels to never be considered together in
+the higher-level analysis. If the image is processed or the detector only
+has one amplifier, you can set the number of channels in both axes to 1.
+
+The final tesselation can be seen on the image with the
address@hidden option that is available to all programs which use
+tessellation for localized operations. When this option is called, a FITS
+file with a @file{_mesh.fits} suffix will be created along with the
+outputs, see @ref{Automatic output}. Each pixel in this image has the ID of
+the tile that covers it. If the number of channels in any dimension are
+larger than unity, you will notice that the tile IDs are defined such that
+the first channels is covered first, then the second and so on.
@menu
* Quantifying signal in a mesh:: Good meshs for gradients.
@@ -11124,7 +11131,7 @@ the `Invoking ProgramName' section within the list of
options.
@table @option
@item -s INT
address@hidden --meshsize=INT
address@hidden --tilesize=INT
The size of each mesh, see @ref{Tiling an image}. If the width of all
channels are not an exact multiple of the specified size, then the last
mesh on each axis will have a different size to cover the full channel.
diff --git a/lib/convolve.c b/lib/convolve.c
index 47cc530..81bf27e 100644
--- a/lib/convolve.c
+++ b/lib/convolve.c
@@ -56,9 +56,9 @@ convolve_tile_is_on_edge(size_t *h, size_t *start_end_coord,
size_t *k,
length along that dimension, then the tile is on the edge.
If the end of the tile in this dimension added with the half-kernel
- length along that dimension is larger than the host's size, then the
- tile is on the edge. */
- do if( (*s++ < *k/2) || (*e++ + *k/2 > *h++) ) return 1; while(++k<kf);
+ length along that dimension is equal to or larger than the host's
+ size, then the tile is on the edge. */
+ do if( (*s++ < *k/2) || (*e++ + *k/2 >= *h++) ) return 1; while(++k<kf);
/* If control reaches here, then this is a central tile. */
return 0;
@@ -97,6 +97,143 @@ struct spatial_params
+/* Define the overlap of the kernel and image over this part of the image,
+ the necessary input image parameters are stored in `overlap' (its
+ `array' and `dsize' elements). The pointer to the kernel array that the
+ overlap starts will be returned. */
+int
+convolve_spatial_overlap(int on_edge, size_t *parse_coords, size_t *hsize,
+ gal_data_t *block, gal_data_t *kernel,
+ gal_data_t *overlap, size_t *k_start_inc,
+ size_t tile_ind)
+{
+ size_t overlap_inc, ndim=kernel->ndim;
+ int full_overlap=1, dim_full_overlap, is_start, is_end;
+
+ size_t *pix = parse_coords; /* needs 2*ndim, also for end. */
+ size_t *overlap_start = parse_coords + ndim * 3;
+ size_t *kernel_start = parse_coords + ndim * 4;
+ size_t *host_start = parse_coords + ndim * 5;
+
+ size_t *p=pix, *pf=pix+ndim, *os=overlap_start, *h=hsize;
+ size_t *k=kernel->dsize, *ks=kernel_start, *od=overlap->dsize;
+ /*
+ if(tile_ind==2053)
+ printf("pix: %zu, %zu\n", pix[0], pix[1]);
+ */
+ /* Coordinate to start convolution for this pixel. */
+ do
+ {
+ /* Initialize the overlap for this dimension (we'll assume it
+ overlaps because this is the most common case usually). */
+ dim_full_overlap=1;
+
+
+ /* When the tile is on the edge, still some pixels in it can have
+ full overlap. So using the `dim_full_overlap', we will do the same
+ thing we do for the tiles that don't overlap for them. */
+ if(on_edge)
+ {
+ /* See if this pixel is on the start and/or end of the dimension
+ relative to the kernel. */
+ is_start = *p < *k/2;
+ is_end = *p + *k/2 >= *h;
+ if( is_start || is_end )
+ {
+ /* Overlapping with the outside.
+
+ With the start: assume that in this dimension, the pixel
+ is at position 2, while the kernel is 11 pixels wide (or 5
+ pixels in half-width). As seen below, the kernel should
+ start from pixel `5-2=3' in this dimension and the overlap
+ size should decrease by the same amount.
+
+ image: 0 1 2 3 4 5 6 7 8 9 ...
+ pixel: p
+ kernel: 0 1 2 3 4 5 6 7 8 9 10
+
+ With the end: Similar to above, but assume the pixel is
+ two pixels away from the edge of a 100-pixel image. We are
+ no longer worried about the overlap or kernel starting
+ point, it is the width that we need to decrease it by:
+
+ 97 + 5 - 100 + 1 : The `1' is because we want the pixel
+ immediately after the end.
+
+ image: ... 92 93 94 95 96 97 98 99 | 100 101 102
+ pixel: p
+ kernel: 0 1 2 3 4 5 6 7 8 | 9 10 */
+ *ks++ = is_start ? *k/2 - *p : 0;
+ *os++ = is_start ? 0 : *p - *k/2;
+
+ /* We will start with the full kernel width, then decrease it
+ if the pixel is too close to the start or end along this
+ dimension. Note that the host array/image might actually
+ be smaller than kernel, so both cases might occur. */
+ *od = *k;
+ if(is_start) *od -= *k/2 - *p;
+ if(is_end) *od -= *p + *k/2 - *h + 1;
+
+ /* Increment and finalize. */
+ ++h;
+ ++k;
+ ++od;
+ full_overlap=0;
+ dim_full_overlap=0;
+ }
+ }
+
+ /* There is full overlap for this pixel or tile over this
+ dimension. */
+ if(dim_full_overlap)
+ {
+ /* Set the values. */
+ *ks++ = 0;
+ *od++ = *k;
+ *os++ = *p - *k/2;
+
+ /* Increment. */
+ ++h;
+ ++k;
+ }
+ }
+ while(++p<pf);
+ /*
+ if(tile_ind==2053)
+ {
+ printf("\tk/2: %zu, %zu\n", kernel->dsize[0]/2, kernel->dsize[1]/2);
+ printf("\toverlap_start: %zu, %zu\n", overlap_start[0],
+ overlap_start[1]);
+ printf("\toverlap->dsize: %zu, %zu\n", overlap->dsize[0],
+ overlap->dsize[1]);
+ exit(0);
+ }
+ */
+ /* Set the increment to start working on the kernel. */
+ *k_start_inc = ( full_overlap
+ ? 0
+ : gal_multidim_coord_to_index(ndim, kernel->dsize,
+ kernel_start) );
+
+ /* Add the host's starting location (necessary when convolution over the
+ host/channel is treated independently). Until now we worked as if the
+ the host/channel is the full image so the edges don't get mixed. But
+ from now on we will be working over the allocated block to look at
+ pixel values, so we need to convert the location to the proper place
+ within the allocated array. */
+ gal_multidim_add_coords(overlap_start, host_start, overlap_start, ndim);
+
+ /* Set the increment to start working on the overlap region and use that
+ to set the starting pointer of the overlap region. */
+ overlap_inc=gal_multidim_coord_to_index(ndim, block->dsize, overlap_start);
+ overlap->array=gal_data_ptr_increment(block->array, overlap_inc,
+ block->type);
+ return full_overlap;
+}
+
+
+
+
/* Convolve over one tile that is not touching the edge. */
static void
@@ -104,17 +241,20 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t
tile_ind,
size_t *parse_coords, gal_data_t *overlap)
{
gal_data_t *tile=&cprm->tiles[tile_ind];
+ gal_data_t *block=gal_tile_block(tile), *kernel=cprm->kernel;
- int on_edge;
double sum, ksum;
+ int on_edge, full_overlap;
size_t i, ndim=tile->ndim, csize=tile->dsize[ndim-1];
- gal_data_t *block=gal_tile_block(tile), *kernel=cprm->kernel;
+ gal_data_t *host=cprm->convoverch ? block : tile->block;
- size_t i_inc, i_ninc, i_st_en[2], o_inc, o_ninc, o_st_en[2];
- size_t *p, *pf, *k, *o, *e, o_st_inc, k_start, start_fastdim;
+ /* Variables for scanning a tile (`i_*') and the region around every
+ pixel of a tile (`o_*'). */
+ size_t start_fastdim, k_start_inc;
+ size_t k_inc, i_inc, i_ninc, i_st_en[2], o_inc, o_ninc, o_st_en[2];
/* These variables depend on the type of the input. */
- float *kv, *iv, *ivf, *i_start, *o_start;
+ float *kv, *iv, *ivf, *i_start, *o_start, *k_start;
float *in_v, *in=block->array, *out=cprm->out->array;
/* The `parse_coords' array was allocated once before this function for
@@ -122,10 +262,9 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t
tile_ind,
necessary coordinates. At first, `pix' will be the starting element of
the tile, but then it will be incremented as we carpet the tile. So we
aren't calling it `start'. */
- size_t *pix = parse_coords;
- size_t *end = parse_coords + ndim;
- size_t *overlap_start = parse_coords + ndim * 2;
- size_t *o_c = parse_coords + ndim * 3;
+ size_t *pix = parse_coords; /* needs 2*ndim (also for end). */
+ size_t *o_c = parse_coords + ndim * 2;
+ size_t *host_start = parse_coords + ndim * 5;
/* Set the starting and ending coordinates of this tile (recall that the
@@ -136,16 +275,21 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t
tile_ind,
gal_tile_start_end_coord(tile, parse_coords, cprm->convoverch);
start_fastdim = pix[ndim-1];
+ /* Set the starting coordinate of the host for this tile. */
+ gal_tile_start_coord(host, host_start);
/* See if this tile is on the edge or not. */
- on_edge=convolve_tile_is_on_edge( ( cprm->convoverch
- ? block->dsize
- : tile->block->dsize ),
- parse_coords, kernel->dsize, ndim);
-
-
- /* Go over the tile, first find its limits, then loop over the fastest
- dimension. */
+ on_edge=convolve_tile_is_on_edge(host->dsize, parse_coords, kernel->dsize,
+ ndim);
+ /*
+ if(tile_ind==2053)
+ {
+ printf("\ntile %zu...\n", tile_ind);
+ printf("\tpix: %zu, %zu\n", pix[0], pix[1]);
+ exit(0);
+ }
+ */
+ /* Go over the tile. */
i_inc=0; i_ninc=1;
i_start=gal_tile_start_end_ind_inclusive(tile, block, i_st_en);
while( i_st_en[0] + i_inc <= i_st_en[1] )
@@ -167,38 +311,32 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t
tile_ind,
out[ in_v - in ]=NAN;
else
{
- /* Coordinate to start convolution for this pixel. */
- pf=(p=pix)+ndim; e=end;
- o=overlap_start; k=kernel->dsize;
- do
- {
- if(on_edge)
- {
- return;
- }
- else { *o++ = *p - *k++/2; k_start=0; }
- }
- while(++p<pf);
- o_st_inc=gal_multidim_coord_to_index(ndim, block->dsize,
- overlap_start);
-
- /* Set the overlap parameters, then set the region. */
- overlap->dsize=kernel->dsize;
- overlap->array=gal_data_ptr_increment(block->array, o_st_inc,
- block->type);
+ /* Define the overlap region. */
+ full_overlap=convolve_spatial_overlap(on_edge, parse_coords,
+ host->dsize, block,
+ kernel, overlap,
+ &k_start_inc, tile_ind);
+
+ /* Set the starting pixel over the image (`o_start'). */
o_start=gal_tile_start_end_ind_inclusive(overlap, block,
o_st_en);
+ /* Set the starting kernel pixel. Note that `kernel_array' is
+ `void *' (pointer arithmetic is not defined on it). So we
+ will first put it in `k_start, and then increment that. */
+ k_start=kernel->array; k_start+=k_start_inc;
+
/* Go over the kernel-overlap region. */
- kv = ( k_start
- ? gal_data_ptr_increment(kernel->array, k_start,
- kernel->type)
- : kernel->array );
- ksum=sum=0.0f; o_inc=0; o_ninc=1;
+ ksum=sum=0.0f; k_inc=0; o_inc=0; o_ninc=1; kv=k_start;
while( o_st_en[0] + o_inc <= o_st_en[1] )
{
- /* Go over the contiguous region. */
+ /* Go over the contiguous region. When there is full
+ overlap, we don't need to calculate incrementation
+ over the kernel, it is always a single
+ incrementation. But when we have partial overlap,
+ we'll need to calculate a different incrementation. */
ivf = ( iv = o_start + o_inc ) + overlap->dsize[ndim-1];
+ if(full_overlap==0) kv = k_start + k_inc;
do
{
if( !isnan(*iv) )
@@ -208,9 +346,13 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t
tile_ind,
while(++iv<ivf);
/* Update the incrementation to the next contiguous
- region of memory over this tile. */
+ region of memory over this tile. Note that the
+ contents of `o_c' are irrelevant here.*/
o_inc += gal_tile_block_increment(block, overlap->dsize,
o_ninc++, o_c);
+ if(full_overlap==0)
+ k_inc += gal_tile_block_increment(kernel, overlap->dsize,
+ o_ninc-1, NULL);
}
/* Set the output value. */
@@ -225,6 +367,10 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t
tile_ind,
contiguous patch. */
i_inc += gal_tile_block_increment(block, tile->dsize, i_ninc++, pix);
}
+ /*
+ if(tile_ind==2053)
+ printf("... done.\n");
+ */
}
@@ -241,7 +387,7 @@ convolve_spatial_on_thread(void *inparam)
size_t i;
gal_data_t overlap, *block=gal_tile_block(cprm->tiles);
size_t *parse_coords=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T,
- 4*block->ndim);
+ 6*block->ndim);
/* Initialize the necessary overlap parameters. */
overlap.block=block;
@@ -256,6 +402,8 @@ convolve_spatial_on_thread(void *inparam)
/* Clean up, wait until all other threads finish, then return. In a
single thread situation, `tprm->b==NULL'. */
+ free(parse_coords);
+ free(overlap.dsize);
if(tprm->b) pthread_barrier_wait(tprm->b);
return NULL;
}
@@ -286,7 +434,6 @@ gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel,
out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, block->ndim, block->dsize,
block->wcs, 0, block->minmapsize, "CONVOLVED",
block->unit, NULL);
- { float *f=out->array, *ff=f+out->size; do *f=NAN; while(++f<ff); }
/* Set the pointers in the parameters structure. */
params.out=out;
diff --git a/lib/tile.c b/lib/tile.c
index b2cad46..1fceecf 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -89,21 +89,36 @@ gal_tile_start_coord(gal_data_t *tile, size_t *start_coord)
void
gal_tile_start_end_coord(gal_data_t *tile, size_t *start_end, int rel_block)
{
- size_t start_ind, ndim=tile->ndim;
+ size_t *s, *sf, *h;
gal_data_t *block=gal_tile_block(tile);
gal_data_t *host=rel_block ? block : tile->block;
+ size_t *hcoord, start_ind, ndim=tile->ndim, *end=start_end+ndim;
/* Get the starting index. Note that for the type we need the allocated
block dataset and can't rely on the tiles. */
- start_ind=gal_data_ptr_dist(host->array, tile->array, block->type);
+ start_ind=gal_data_ptr_dist(block->array, tile->array, block->type);
- /* Get the coordinates of the starting point. */
- gal_multidim_index_to_coord(start_ind, tile->ndim, host->dsize,
- start_end);
+ /* Get the coordinates of the starting point relative to the allocated
+ block. */
+ gal_multidim_index_to_coord(start_ind, ndim, block->dsize, start_end);
+
+ /* When the host is different from the block, the tile's starting
+ position needs to be corrected. */
+ if(host!=block)
+ {
+ /* Get the host's starting coordinates. */
+ start_ind=gal_data_ptr_dist(block->array, host->array, block->type);
+
+ /* Temporarily put the host's coordinates in the place held for the
+ ending coordinates. */
+ hcoord=end;
+ gal_multidim_index_to_coord(start_ind, ndim, block->dsize, hcoord);
+ sf=(s=start_end)+ndim; h=hcoord; do *s++ -= *h++; while(s<sf);
+ }
/* Add the dimensions of the tile to the starting coordinate. Note that
the ending coordinates are stored immediately after the start.*/
- gal_multidim_add_coords(start_end, tile->dsize, start_end+ndim, ndim);
+ gal_multidim_add_coords(start_end, tile->dsize, end, ndim);
}
@@ -264,13 +279,13 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
/* 1D: the increment is just the tile size. */
case 1:
increment=t[0];
- coord[0]+=increment;
+ if(coord) coord[0]+=increment;
break;
/* 2D: the increment is the block's number of fastest axis pixels. */
case 2:
increment=b[1];
- ++coord[0];
+ if(coord) ++coord[0];
break;
/* Higher dimensions. */
@@ -278,13 +293,13 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
if(num_increment % t[n-2])
{
increment=b[n-1];
- ++coord[n-2];
+ if(coord) ++coord[n-2];
}
else
{
increment=(b[n-2] * b[n-1]) - ( (t[n-2]-1) * b[n-1] );
++coord[n-3];
- coord[n-2]=coord[n-1]=0;
+ if(coord) coord[n-2]=coord[n-1]=0;
}
break;
}
@@ -422,7 +437,8 @@ gal_tile_all_sanity_check(char *filename, char *hdu,
gal_data_t *input,
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);
+ "value(s) given for the tile size (`--tilesize' option).",
+ filename, hdu, input->ndim, i);
/* Check the channel's dimensions. */
- [gnuastro-commits] master 5eb817a 052/125: Option values stored in argp_option, (continued)
- [gnuastro-commits] master 5eb817a 052/125: Option values stored in argp_option, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 7e888ce 073/125: First implementation of MakeProfiles with gal_data_t, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e8d9751 087/125: MakeNoise now uses the new gal_data_t infrastructure, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e978d44 068/125: Numeric option sanity checks, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e74eb25 093/125: Changed name of Header program to Fits, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e353932 109/125: NaN magnitude will only blank the profile's area, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 96ab804 098/125: Working on tiles implemented in tile checking, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 31b3b76 099/125: Spatial convolution now using new tessellation, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master c0dd6db 102/125: The noedgecorrection option in spatial domain convolution, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0dfdc34 069/125: Option values directly set in library, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 95ba655 100/125: Spatial convolution with new tessellation,
Mohammad Akhlaghi <=
- [gnuastro-commits] master 0dec9ac 110/125: Allow for floating point errors in Crop comparisons, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master db4c16d 076/125: Pulled in all the recent work in master, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 64f8db5 080/125: ConvertType working except for text inputs and outputs, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 4d853b4 112/125: Old neighbors.h header removed, used new in MakeProfiles, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 9959e45 004/125: Arithmetic using new data structure, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0a2fd0e 089/125: New implementations for histogram and CFP, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 213d1c3 106/125: gal_data_copy_* functions can now work on tiles, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 4f6a8c4 095/125: Fits prints HDU info of FITS file with no options, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 63043e5 105/125: Statistical functions can accept tiled input, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 4839a11 104/125: Preparations for tile interpolation, Mohammad Akhlaghi, 2017/04/23