[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master a1ff97b 111/125: Permutations and neighbor fin
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master a1ff97b 111/125: Permutations and neighbor finding algorithms |
|
Date: |
Sun, 23 Apr 2017 22:36:49 -0400 (EDT) |
branch: master
commit a1ff97be6ff237e0b574df4fc0e09aed25fb368e
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Permutations and neighbor finding algorithms
As part of finalizing the interpolation, it was necessary to rebuild some
infra-structure which has been done in this commit. The two main necessary
infrastructures were: Neighbor-finding (for each data-element/pixel) and
permutations. They are both described below:
Until now, the neighbor finding algorithm would need to actually check
every neighbor to see if it is outside the dataset or not. With the new
implementation, it is only necessary to check if the pixel is on the edge
of the image or not just once. The result is then stored as bits, so more
distant neighbors can be checked using them (instead of checking 8-byte
indexs). Like before, it was also defined as a macro, in the
`gnuastro/dimension.h' (new name for `gnuastro/multidim.h') header.
Permutations are necessary because tiles are defined contiguously over each
channel (not over the whole input dataset). So when it was necessary to do
operations over the whole image (for example interpolating when the
channels are to be ignored, or saving as a FITS file), the tile values
would not be displayed properly. Until now, we didn't use permutations, but
a very complicated and slow process of building a new array.
The main reason permutations were chosen was that after doing some
operations it was necessary to reverse the transformation above and put the
values back in tile-order. But in the old implementation, it was necessary
to write a completely new large and complicated function which was time
consuming and also buggy. So permutations were chosen to do the job,
because after defining one permutation, it is possible to apply itself or
its inverse with some very simple functions. The GNU Scientific Library
(GSL) had some great functions for the job. But unfortunately GSL works
based on subjective types like `int' and `long' not fixed width types. So
it was not easy to convert it to Gnuastro's types. Therefore, its functions
were re-implemented in a type-agnostic manner (using `memcpy') as part of
Gnuastro in a new `gnuastro/permutation.h' header.
As part of the process the following changes were also made:
- `gal_data_bit_print_stream' is a new function to print a bit stream of
any address in memory.
- `gal_dimension_increment' is a new function to calculate the distance
necessary to increment an index along all dimensions.
- `gal_dimension_num_neighbors' is a new function to return the maximum
number of neighbors to a data-element/pixel, given a specific dimension.
- Thanks to `gal_dimension_increment', `gal_dimension_index_to_coord' has
been re-implemented in a more robust/clean way.
- `gal_tile_full_interpolate' is a new function to do interpolation over
tiles. Since tile value interpolation is commonly necessary in a lot of
places in Gnuastro and it is not too trivial, this function will do
everything for any application.
- Only some ground-work has started for the `lib/interpolate.c'
library. But thanks to the major new infra-structure of this commit,
work will go on much more smoothly after this commit.
---
bin/statistics/args.h | 16 ++
bin/statistics/main.h | 1 +
bin/statistics/statistics.c | 33 ++-
bin/statistics/ui.h | 3 +-
lib/Makefile.am | 20 +-
lib/convolve.c | 10 +-
lib/data.c | 21 ++
lib/{multidim.c => dimension.c} | 93 +++++---
lib/gnuastro/data.h | 3 +
lib/gnuastro/dimension.h | 348 +++++++++++++++++++++++++++++
lib/gnuastro/{multidim.h => interpolate.h} | 29 +--
lib/gnuastro/{multidim.h => permutation.h} | 35 +--
lib/gnuastro/tile.h | 10 +-
lib/interpolate.c | 61 +++++
lib/permutation.c | 188 ++++++++++++++++
lib/tile.c | 303 +++++++++++++++----------
lib/wcs.c | 4 +-
17 files changed, 962 insertions(+), 216 deletions(-)
diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index 3d3f099..40eb667 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -100,6 +100,22 @@ struct argp_option program_options[] =
+ /* Tessellation */
+ {
+ "interpolate",
+ ARGS_OPTION_KEY_INTERPOLATE,
+ 0,
+ 0,
+ "Interpolate over blank tiles to fill them.",
+ GAL_OPTIONS_GROUP_TESSELLATION,
+ &p->interpolate,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+
+
{
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index a409f8c..3c449c1 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -65,6 +65,7 @@ struct statisticsparams
float quantmin; /* Quantile min or range: from Q to 1-Q. */
float quantmax; /* Quantile maximum. */
uint8_t ontile; /* Do single value calculations on tiles. */
+ uint8_t interpolate; /* Use interpolation to fill blank tiles. */
uint8_t asciihist; /* Print an ASCII histogram. */
uint8_t asciicfp; /* Print an ASCII cumulative frequency plot.*/
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index f1fa283..6a104f3 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -32,9 +32,12 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/fits.h>
#include <gnuastro/tile.h>
+#include <gnuastro/blank.h>
#include <gnuastro/multidim.h>
#include <gnuastro/arithmetic.h>
#include <gnuastro/statistics.h>
+#include <gnuastro/interpolate.h>
+#include <gnuastro/permutation.h>
#include <timing.h>
#include <checkset.h>
@@ -230,7 +233,29 @@ statistics_print_one_row(struct statisticsparams *p)
/*******************************************************************/
/************** Single value on tile ***************/
/*******************************************************************/
-void
+static void
+statistics_interpolate_and_write(struct statisticsparams *p,
+ gal_data_t *values)
+{
+ char *output=p->cp.output;
+ struct gal_tile_two_layer_params *tl=&p->cp.tl;
+
+ /* Do the interpolation (if necessary). */
+ if(p->interpolate && gal_blank_present(values))
+ gal_tile_full_interpolate(values, tl);
+
+ /* Write the values. */
+ gal_tile_full_values_write(values, tl, output, PROGRAM_STRING);
+
+ /* Clean up. */
+ gal_data_free(values);
+}
+
+
+
+
+
+static void
statistics_on_tile(struct statisticsparams *p)
{
double arg;
@@ -378,9 +403,9 @@ statistics_on_tile(struct statisticsparams *p)
gal_data_free(tmp);
}
- /* Save the values. */
- gal_tile_full_values_write(values, tl, cp->output, PROGRAM_STRING);
- gal_data_free(values);
+ /* Do the interpolation (if necessary) and write the array into the
+ output. */
+ statistics_interpolate_and_write(p, values);
/* Clean up. */
if(operation->v==ARGS_OPTION_KEY_QUANTFUNC) gal_data_free(tmpv);
diff --git a/bin/statistics/ui.h b/bin/statistics/ui.h
index aa55c42..b681f4f 100644
--- a/bin/statistics/ui.h
+++ b/bin/statistics/ui.h
@@ -29,7 +29,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Available letters for short options:
- a b e f i j k p v w x y z
+ a b e f j k p v w x y z
B G J L R W X Y
*/
enum option_keys_enum
@@ -51,6 +51,7 @@ enum option_keys_enum
ARGS_OPTION_KEY_SIGMACLIP = 's',
ARGS_OPTION_KEY_NORMALIZE = 'n',
ARGS_OPTION_KEY_ONTILE = 't',
+ ARGS_OPTION_KEY_INTERPOLATE = 'i',
/* Only with long version (start with a value 1000, the rest will be set
automatically). */
diff --git a/lib/Makefile.am b/lib/Makefile.am
index b3b4c73..180f879 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -47,10 +47,11 @@ 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 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
+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 interpolate.c linkedlist.c options.c permutation.c polygon.c \
+ qsort.c dimension.c statistics.c table.c threads.c tile.c timing.c txt.c \
+ wcs.c
@@ -62,12 +63,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)/convolve.h \
- $(headersdir)/data.h $(headersdir)/fits.h $(headersdir)/git.h
\
- $(headersdir)/linkedlist.h $(headersdir)/multidim.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)/interpolate.h $(headersdir)/linkedlist.h \
+ $(headersdir)/dimension.h $(headersdir)/permutation.h \
$(headersdir)/polygon.h $(headersdir)/qsort.h $(headersdir)/statistics.h \
- $(headersdir)/table.h $(headersdir)/threads.h $(headersdir)/tile.h \
+ $(headersdir)/table.h $(headersdir)/threads.h $(headersdir)/tile.h \
$(headersdir)/wcs.h $(headersdir)/txt.h
diff --git a/lib/convolve.c b/lib/convolve.c
index d741a40..e3b1e3d 100644
--- a/lib/convolve.c
+++ b/lib/convolve.c
@@ -30,8 +30,8 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/tile.h>
#include <gnuastro/threads.h>
-#include <gnuastro/multidim.h>
#include <gnuastro/convolve.h>
+#include <gnuastro/dimension.h>
@@ -213,8 +213,8 @@ convolve_spatial_overlap(int on_edge, size_t *parse_coords,
size_t *hsize,
/* 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) );
+ : gal_dimension_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
@@ -222,11 +222,11 @@ convolve_spatial_overlap(int on_edge, size_t
*parse_coords, size_t *hsize,
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);
+ gal_dimension_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_inc=gal_dimension_coord_to_index(ndim, block->dsize, overlap_start);
overlap->array=gal_data_ptr_increment(block->array, overlap_inc,
block->type);
return full_overlap;
diff --git a/lib/data.c b/lib/data.c
index a683ca0..e23ed40 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -280,6 +280,27 @@ gal_data_out_type(gal_data_t *first, gal_data_t *second)
/*********************************************************************/
/************* Size and allocation *******************/
/*********************************************************************/
+/* Print the contents of any place in memory that `in' points to for `size'
+ bytes. This solution was inspired from:
+
+
http://stackoverflow.com/questions/111928/is-there-a-printf-converter-to-print-in-binary-format
*/
+void
+gal_data_bit_print_stream(void *in, size_t size)
+{
+ size_t i;
+ char *byte=in;
+ for(i=0;i<size;++i)
+ printf("%c%c%c%c%c%c%c%c ",
+ (byte[i] & 0x80 ? '1' : '0'), (byte[i] & 0x40 ? '1' : '0'),
+ (byte[i] & 0x20 ? '1' : '0'), (byte[i] & 0x10 ? '1' : '0'),
+ (byte[i] & 0x08 ? '1' : '0'), (byte[i] & 0x04 ? '1' : '0'),
+ (byte[i] & 0x02 ? '1' : '0'), (byte[i] & 0x01 ? '1' : '0') );
+}
+
+
+
+
+
int
gal_data_dsize_is_different(gal_data_t *first, gal_data_t *second)
{
diff --git a/lib/multidim.c b/lib/dimension.c
similarity index 65%
rename from lib/multidim.c
rename to lib/dimension.c
index d092770..edc6672 100644
--- a/lib/multidim.c
+++ b/lib/dimension.c
@@ -1,5 +1,5 @@
/*********************************************************************
-multidim -- Functions for multi-dimensional operations.
+dimension -- Functions for multi-dimensional operations.
This is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
@@ -22,12 +22,13 @@ 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 <stdlib.h>
-#include <gnuastro/multidim.h>
+#include <gnuastro/dimension.h>
@@ -37,7 +38,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/******************** Info **********************/
/************************************************************************/
size_t
-gal_multidim_total_size(size_t ndim, size_t *dsize)
+gal_dimension_total_size(size_t ndim, size_t *dsize)
{
size_t i, num=1;
for(i=0;i<ndim;++i) num *= dsize[i];
@@ -48,6 +49,43 @@ gal_multidim_total_size(size_t ndim, size_t *dsize)
+/* Calculate the values necessary to increment/decrement along each
+ dimension of a dataset with size `dsize'. */
+size_t *
+gal_dimension_increment(size_t ndim, size_t *dsize)
+{
+ int i;
+ size_t *out=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
+
+ /* Along the fastest dimension, it is 1. */
+ out[ndim-1]=1;
+
+ /* For the rest of the dimensions, it is */
+ if(ndim>1)
+ for(i=ndim-2;i>=0;--i)
+ out[i]=dsize[i+1]*out[i+1];
+
+ /* Return the allocated array. */
+ return out;
+}
+
+
+
+
+
+size_t
+gal_dimension_num_neighbors(size_t ndim)
+{
+ if(ndim)
+ return pow(3, ndim)-1;
+ else
+ error(EXIT_FAILURE, 0, "ndim cannot be zero in "
+ "`gal_dimension_num_neighbors'");
+ return 0;
+}
+
+
+
@@ -67,7 +105,7 @@ gal_multidim_total_size(size_t ndim, size_t *dsize)
/******************** Coordinates **********************/
/************************************************************************/
void
-gal_multidim_add_coords(size_t *c1, size_t *c2, size_t *out, size_t ndim)
+gal_dimension_add_coords(size_t *c1, size_t *c2, size_t *out, size_t ndim)
{
size_t *end=c1+ndim;
do *out++ = *c1++ + *c2++; while(c1<end);
@@ -80,14 +118,14 @@ gal_multidim_add_coords(size_t *c1, size_t *c2, size_t
*out, size_t ndim)
/* 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)
+gal_dimension_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 "
+ error(EXIT_FAILURE, 0, "`gal_dimension_coord_to_index' doesn't accept "
"0 dimensional arrays");
case 1:
@@ -129,16 +167,16 @@ gal_multidim_coord_to_index(size_t ndim, size_t *dsize,
size_t *coord)
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)
+gal_dimension_index_to_coord(size_t ind, size_t ndim, size_t *dsize,
+ size_t *coord)
{
- size_t d, i, in_all_faster_dim;
+ size_t d, *dinc;
switch(ndim)
{
case 0:
error(EXIT_FAILURE, 0, "a 0-dimensional dataset is not defined in "
- "`gal_multidim_ind_to_coord'");
+ "`gal_dimension_ind_to_coord'");
/* One dimensional dataset. */
case 1:
@@ -153,29 +191,24 @@ gal_multidim_index_to_coord(size_t ind, size_t ndim,
size_t *dsize,
/* Higher dimensional datasets. */
default:
- /* We start with the slowest dimension (first in the C standard). */
+ /* Set the incrementation values for each dimension. */
+ dinc=gal_dimension_increment(ndim, dsize);
+
+ /* We start with the slowest dimension (first in the C standard) and
+ continue until (but not including) the fastest dimension. This is
+ because except for the fastest (coniguous) dimension, the other
+ coordinates can be found by division. */
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];
+ /* Set the coordinate value for this dimension. */
+ coord[d] = ind / dinc[d];
- /* 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;
- }
+ /* Replace the index with its remainder with the number of
+ elements in all faster dimensions. */
+ ind %= dinc[d];
}
+
+ /* Clean up. */
+ free(dinc);
}
}
diff --git a/lib/gnuastro/data.h b/lib/gnuastro/data.h
index de9cc89..4663c49 100644
--- a/lib/gnuastro/data.h
+++ b/lib/gnuastro/data.h
@@ -233,6 +233,9 @@ typedef struct gal_data_t
/*************************************************************
************** Type information ***************
*************************************************************/
+void
+gal_data_bit_print_stream(void *in, size_t size);
+
char *
gal_data_type_as_string(uint8_t type, int long_name);
diff --git a/lib/gnuastro/dimension.h b/lib/gnuastro/dimension.h
new file mode 100644
index 0000000..b27ff70
--- /dev/null
+++ b/lib/gnuastro/dimension.h
@@ -0,0 +1,348 @@
+/*********************************************************************
+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_dimension_total_size(size_t ndim, size_t *dsize);
+
+size_t *
+gal_dimension_increment(size_t ndim, size_t *dsize);
+
+size_t
+gal_dimension_num_neighbors(size_t ndim);
+
+
+
+
+
+/************************************************************************/
+/******************** Coordinates **********************/
+/************************************************************************/
+void
+gal_dimension_add_coords(size_t *c1, size_t *c2, size_t *out, size_t ndim);
+
+size_t
+gal_dimension_coord_to_index(size_t ndim, size_t *dsize, size_t *coord);
+
+void
+gal_dimension_index_to_coord(size_t ind, size_t ndim, size_t *dsize,
+ size_t *coord);
+
+
+
+
+
+/************************************************************************/
+/******************** Neighbors **********************/
+/************************************************************************/
+/* Purpose
+ -------
+
+ This macro will allow you to do a fixed operation on the neighbors of an
+ eleemnt. The identifier for the element is its dimension-agnostic index
+ (distance from start of array). It is defined as a macro (and not a
+ function) it is often necessary to loop over a very large number of
+ pixels/indexs and the number of neighbors differs (in different
+ dimensions and on the edges of the image).
+
+ Usage
+ -----
+
+ The inputs are:
+
+ `index': The index of the desired point.
+
+ `ndim': The number of dimensions in the dataset.
+
+ `dsize': The size of the dataset along each dimension (C order).
+
+ `connectivity': The connectivity of the neighbors. This is a single
+ integer with a value between `1' and `ndim' (it is
+ irrelevant for a 1D dataset). Below, the case for 3D data
+ is shown. Each `dn' corresponds to `dinc[n]'. See `dinc'
+ below for the `dinc' values.
+
+ `1': At most ONE addition/subtraction. In 1D: only the
+ first line should be used. In 2D, this is 4
+ connectivity. In 3D all the elements that share an
+ 2D-face with the cube at index `i'.
+
+ i - d0 i + d0 (1D, 2D & 3D)
+ i - d1 i + d1 (2D & 3D)
+ i - d2 i + d2 (3D)
+
+ `2': At most TWO additions/subtractions. In 2D: 8
+ connectivity. In 3D: all elements that share a 1D
+ edge with the cube at index `i'.
+
+ i - d0 - d1 i - d0 + d1 (2D & 3D)
+ i + d0 - d1 i + d0 + d1 (2D & 3D)
+ i - d0 - d2 i - d0 + d2 (3D)
+ i + d0 - d2 i + d0 + d2 (3D)
+ i - d1 - d2 i - d1 + d2 (3D)
+ i + d1 - d2 i + d1 + d2 (3D)
+
+ `3': At most THREE additions/subtractions (only for 3D
+ and higher dimensions). All cubes that share a 0-D
+ vertex with the cube at index `i'.
+
+ i - d0 - d1 - d2 i - d0 - d1 + d2
+ i - d0 + d1 - d2 i - d0 + d1 + d2
+ i + d0 - d1 - d2 i + d0 - d1 + d2
+ i + d0 + d1 - d2 i + d0 + d1 + d2
+
+ `dinc': An array keeping the length necessary to increment along
+ each dimension. You can make this array with the following
+ function:
+
+ dinc=gal_dimension_increment(ndim, dsize);
+
+ `operation': Any C operation. `nind' is a `size_t' type variable
+ that is defined by this macro and will have the index of
+ each neighbor. You can use this `nind' for any processing
+ that you like on the neighbor. Note that `op' will be
+ repeated the number of times there is a neighbor.
+
+
+ Implementation
+ --------------
+
+ To be most efficient (avoid as many `if's as possible), we will start
+ parsing the neighbors from the fastest dimension. When-ever the element
+ is on the edge of the dataset in any dimension, we will store it in a
+ bit-wise array (one bit for each dimension, to mark if it is on the edge
+ (either side 1, or in the middle 0). Far many more elements will be in
+ the middle, so there is no problem to check the edge. The good thing
+ with a bit-array is that it can take one register on the CPU and stay
+ there until the end. But if we want to have an array of values, multiple
+ registers will be necessary.
+
+ The bit information is in two two-byte spaces, so in theory, this works
+ for 16 dimensions.
+*/
+#define gal_dimension_neighbor_op(index, ndim, dsize, connectivity, \
+ dinc, operation) { \
+ uint32_t bitstr=0; \
+ size_t nind, ind=index; \
+ uint8_t D, *is_start, *is_end, *is_edge, one=1; \
+ \
+ /* Initialize the start/end. */ \
+ is_start=(uint8_t *)(&bitstr); \
+ is_end=(uint8_t *)(&bitstr)+1; \
+ \
+ /* Start with the slowest dimension and see if it is on the edge */ \
+ /* or not, similar to `gal_dimension_index_to_coord'. In the */ \
+ /* process, also fill the `connectivity==1' neighbors. */ \
+ for(D=0;D<ndim;++D) \
+ { \
+ /* If this dimension is only one element wide, no neighbors. */ \
+ if( (dsize)[D] == 1 ) \
+ { \
+ *is_start |= 1<<D; \
+ *is_end |= 1<<D; \
+ } \
+ else \
+ { \
+ if( ind / (dinc)[D] ) \
+ { \
+ /* We are at the end of this dimension. */ \
+ if( ind / (dinc)[D] == (dsize)[D]-1 ) \
+ { \
+ *is_end |= one<<D; \
+ nind = (index) - (dinc)[D]; {operation;}; \
+ } \
+ \
+ /* Middle of the dimension: both +1 and -1 possible. */ \
+ else \
+ { \
+ nind = (index) - (dinc)[D]; {operation;}; \
+ nind = (index) + (dinc)[D]; {operation;}; \
+ } \
+ } \
+ else \
+ { \
+ *is_start |= one<<D; \
+ nind = (index) + (dinc)[D]; {operation;}; \
+ } \
+ } \
+ \
+ /* Change `ind' to the remainder of previous dimensions. */ \
+ ind %= dinc[D]; \
+ } \
+ \
+ /* We now know if the index is on the edge or not. During the */ \
+ /* process above, we also finished the `connectivity==1' case. */ \
+ /* So we'll just have to find the rest of the terms. */ \
+ if(connectivity>1 && ndim>1) \
+ { \
+ /* Finalize `is_edge' (bit value 1 for respective dim.). */ \
+ is_edge=(uint8_t *)(&bitstr)+2; \
+ *is_edge = *is_start | *is_end; \
+ \
+ /* Shared between 2D and 3D datasets. */ \
+ if(*is_edge) \
+ { /* NOTE: these are bitwise operators, not conditionals. */ \
+ if( !( *is_start & ( one | one<<1 ) ) ) \
+ { nind=(index) - (dinc)[0] - (dinc)[1]; {operation;}; } \
+ if( !( *is_start & one ) && !( *is_end & one<<1 ) ) \
+ { nind=(index) - (dinc)[0] + (dinc)[1]; {operation;}; } \
+ if( !( *is_end & one ) && !( *is_start & one<<1 ) ) \
+ { nind=(index) + (dinc)[0] - (dinc)[1]; {operation;}; } \
+ if( !( *is_end & ( one | one<<1 ) ) ) \
+ { nind=(index) + (dinc)[0] + (dinc)[1]; {operation;}; } \
+ } \
+ else \
+ { \
+ nind=(index) - (dinc)[0] - (dinc)[1]; {operation;}; \
+ nind=(index) - (dinc)[0] + (dinc)[1]; {operation;}; \
+ nind=(index) + (dinc)[0] - (dinc)[1]; {operation;}; \
+ nind=(index) + (dinc)[0] + (dinc)[1]; {operation;}; \
+ } \
+ \
+ /* Only for 3D datasets. */ \
+ if(ndim>2) \
+ { \
+ /* Connectivity == 2. */ \
+ if(*is_edge) \
+ for(D=0;D<2;++D) \
+ { \
+ if( !( *is_start & ( one<<D | one<<2 ) ) ) \
+ { nind=(index) - (dinc)[D] - (dinc)[2]; {operation;}; } \
+ if( !( *is_start & one<<D ) && !( *is_end & one<<2 ) ) \
+ { nind=(index) - (dinc)[D] + (dinc)[2]; {operation;}; } \
+ if( !( *is_end & one<<D ) && !( *is_start & one<<2 ) ) \
+ { nind=(index) + (dinc)[D] - (dinc)[2]; {operation;}; } \
+ if( !( *is_end & ( one<<D | one<<2 ) ) ) \
+ { nind=(index) + (dinc)[D] + (dinc)[2]; {operation;}; } \
+ } \
+ else \
+ for(D=0;D<2;++D) \
+ { \
+ nind=(index) - (dinc)[D] - (dinc)[2]; {operation;}; \
+ nind=(index) - (dinc)[D] + (dinc)[2]; {operation;}; \
+ nind=(index) + (dinc)[D] - (dinc)[2]; {operation;}; \
+ nind=(index) + (dinc)[D] + (dinc)[2]; {operation;}; \
+ } \
+ \
+ /* Connectivity == 3. */ \
+ if(connectivity>2) \
+ { \
+ if(*is_edge) \
+ { \
+ if( !*is_start ) \
+ { nind=(index) - (dinc)[0] - (dinc)[1] - (dinc)[2]; \
+ {operation;}; } \
+ \
+ if( !(*is_start & one) && !(*is_start & one<<1) \
+ && !(*is_end & one<<2)) \
+ { nind=(index) - (dinc)[0] - (dinc)[1] + (dinc)[2]; \
+ {operation;}; } \
+ \
+ if( !(*is_start & one) && !(*is_end & one<<1) \
+ && !(*is_start & one<<2)) \
+ { nind=(index) - (dinc)[0] + (dinc)[1] - (dinc)[2]; \
+ {operation;}; } \
+ \
+ if( !(*is_start & one) && !(*is_end & one<<1) \
+ && !(*is_end & one<<2)) \
+ { nind=(index) - (dinc)[0] + (dinc)[1] + (dinc)[2]; \
+ {operation;}; } \
+ \
+ if( !(*is_end & one) && !(*is_start & one<<1) \
+ && !(*is_start & one<<2)) \
+ { nind=(index) + (dinc)[0] - (dinc)[1] - (dinc)[2]; \
+ {operation;}; } \
+ \
+ if( !(*is_end & one) && !(*is_start & one<<1) \
+ && !(*is_end & one<<2)) \
+ { nind=(index) + (dinc)[0] - (dinc)[1] + (dinc)[2]; \
+ {operation;}; } \
+ \
+ if( !(*is_end & one) && !(*is_end & one<<1) \
+ && !(*is_start & one<<2)) \
+ { nind=(index) + (dinc)[0] + (dinc)[1] - (dinc)[2]; \
+ {operation;}; } \
+ \
+ if( !*is_end ) \
+ { nind=(index) + (dinc)[0] + (dinc)[1] + (dinc)[2]; \
+ {operation;}; } \
+ } \
+ else \
+ { \
+ nind=(index) - (dinc)[0] - (dinc)[1] - (dinc)[2]; \
+ {operation;}; \
+ nind=(index) - (dinc)[0] - (dinc)[1] + (dinc)[2]; \
+ {operation;}; \
+ nind=(index) - (dinc)[0] + (dinc)[1] - (dinc)[2]; \
+ {operation;}; \
+ nind=(index) - (dinc)[0] + (dinc)[1] + (dinc)[2]; \
+ {operation;}; \
+ nind=(index) + (dinc)[0] - (dinc)[1] - (dinc)[2]; \
+ {operation;}; \
+ nind=(index) + (dinc)[0] - (dinc)[1] + (dinc)[2]; \
+ {operation;}; \
+ nind=(index) + (dinc)[0] + (dinc)[1] - (dinc)[2]; \
+ {operation;}; \
+ nind=(index) + (dinc)[0] + (dinc)[1] + (dinc)[2]; \
+ {operation;}; \
+ } \
+ } \
+ } \
+ } \
+ \
+ /* For a check. */ \
+ /* printf("\nEdge bit flags: "); */ \
+ /* gal_data_bit_print_stream(&bitstr, 3); printf("\n"); */ \
+ }
+
+
+__END_C_DECLS /* From C++ preparations */
+
+#endif
diff --git a/lib/gnuastro/multidim.h b/lib/gnuastro/interpolate.h
similarity index 63%
copy from lib/gnuastro/multidim.h
copy to lib/gnuastro/interpolate.h
index 3f37d8f..3c2495a 100644
--- a/lib/gnuastro/multidim.h
+++ b/lib/gnuastro/interpolate.h
@@ -1,5 +1,5 @@
/*********************************************************************
-multidim -- Functions for multi-dimensional operations.
+tile -- work with tesselations over a host dataset.
This is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
@@ -20,8 +20,8 @@ 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__
+#ifndef __GAL_INTERPOLATE_H__
+#define __GAL_INTERPOLATE_H__
/* Include other headers if necessary here. Note that other header files
must be included before the C++ preparations below */
@@ -42,29 +42,8 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* 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 **********************/
-/************************************************************************/
void
-gal_multidim_add_coords(size_t *c1, size_t *c2, size_t *out, size_t ndim);
-
-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);
-
+gal_interpolate(gal_data_t *input);
__END_C_DECLS /* From C++ preparations */
diff --git a/lib/gnuastro/multidim.h b/lib/gnuastro/permutation.h
similarity index 74%
rename from lib/gnuastro/multidim.h
rename to lib/gnuastro/permutation.h
index 3f37d8f..905855d 100644
--- a/lib/gnuastro/multidim.h
+++ b/lib/gnuastro/permutation.h
@@ -1,5 +1,5 @@
/*********************************************************************
-multidim -- Functions for multi-dimensional operations.
+Permutation -- Work on permutations (arrays of indexs).
This is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
@@ -20,8 +20,8 @@ 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__
+#ifndef __GAL_PERMUTATION_H__
+#define __GAL_PERMUTATION_H__
/* Include other headers if necessary here. Note that other header files
must be included before the C++ preparations below */
@@ -43,27 +43,28 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
__BEGIN_C_DECLS /* From C++ preparations */
-/************************************************************************/
-/******************** Info **********************/
-/************************************************************************/
-size_t
-gal_multidim_total_size(size_t ndim, size_t *dsize);
+/*********************************************************************/
+/*************** Permutation info *******************/
+/*********************************************************************/
+void
+gal_permutation_check(size_t *permutation, size_t size);
-/************************************************************************/
-/******************** Coordinates **********************/
-/************************************************************************/
-void
-gal_multidim_add_coords(size_t *c1, size_t *c2, size_t *out, size_t ndim);
-size_t
-gal_multidim_coord_to_index(size_t ndim, size_t *dsize, size_t *coord);
+
+/*********************************************************************/
+/*************** Apply permutation *******************/
+/*********************************************************************/
void
-gal_multidim_index_to_coord(size_t ind, size_t ndim, size_t *dsize,
- size_t *coord);
+gal_permutation_apply(gal_data_t *input, size_t *permutation);
+
+void
+gal_permutation_apply_inverse(gal_data_t *input, size_t *permutation);
+
+
__END_C_DECLS /* From C++ preparations */
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index eb4782a..f4c3420 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -107,6 +107,7 @@ struct gal_tile_two_layer_params
size_t *numtiles; /* Tile no. in each dim. over-all. */
size_t *numtilesinch; /* Tile no. in each dim. on one channel. */
char *tilecheckname; /* Name of file to check tiles. */
+ size_t *permutation; /* Tile pos. in memory --> pos. overall. */
/* Actual tile and channel data structures. */
gal_data_t *tiles; /* Tiles array (also linked with `next'). */
@@ -128,15 +129,18 @@ gal_tile_full_two_layers(gal_data_t *input,
void
gal_tile_full_free_contents(struct gal_tile_two_layer_params *tl);
-gal_data_t *
-gal_tile_full_values_for_disp(gal_data_t *tilevalues,
- struct gal_tile_two_layer_params *tl);
+void
+gal_tile_full_permutation(struct gal_tile_two_layer_params *tl);
void
gal_tile_full_values_write(gal_data_t *tilevalues,
struct gal_tile_two_layer_params *tl,
char *filename, char *program_string);
+void
+gal_tile_full_interpolate(gal_data_t *tilevalues,
+ struct gal_tile_two_layer_params *tl);
+
diff --git a/lib/interpolate.c b/lib/interpolate.c
new file mode 100644
index 0000000..a61a0f2
--- /dev/null
+++ b/lib/interpolate.c
@@ -0,0 +1,61 @@
+/*********************************************************************
+Interpolate - Fill blank values in a 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 <stdlib.h>
+
+#include <gnuastro/data.h>
+#include <gnuastro/fits.h>
+#include <gnuastro/blank.h>
+#include <gnuastro/dimension.h>
+
+
+
+
+
+/* Interpolate blank values in an array. */
+void
+gal_interpolate(gal_data_t *input)
+{
+ size_t F;
+ uint8_t *flagarr;
+ gal_data_t *flag;
+ size_t *dinc=gal_dimension_increment(input->ndim, input->dsize);
+
+ /* Set a value of 1 for every element that must be interpolated. */
+ flag=gal_blank_flag(input);
+
+ /* Find the proper value for each neighbor. */
+ flagarr=flag->array;
+ for(F=0; F<input->size; ++F)
+ if(flagarr[F])
+ {
+ printf("to be filled: %zu\n", F);
+ gal_dimension_neighbor_op(F, input->ndim, input->dsize, 2,
+ dinc, {printf("\tneighbor: %zu\n", nind);});
+ }
+
+ /* Clean up. */
+ gal_data_free(flag);
+}
diff --git a/lib/permutation.c b/lib/permutation.c
new file mode 100644
index 0000000..2675cb3
--- /dev/null
+++ b/lib/permutation.c
@@ -0,0 +1,188 @@
+/*********************************************************************
+Permutation -- Work on permutations (arrays of indexs).
+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 <string.h>
+#include <stdlib.h>
+
+#include <gnuastro/permutation.h>
+
+
+
+/*********************************************************************/
+/*************** Permutation info *******************/
+/*********************************************************************/
+void
+gal_permutation_check(size_t *permutation, size_t size)
+{
+ size_t i;
+ for(i=0; i<size; ++i)
+ printf("out[ %-4zu ] = in [ %-4zu ]\n", i, permutation[i]);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*********************************************************************/
+/*************** Apply permutation *******************/
+/*********************************************************************/
+/* Re-order the input dataset based on the given permutation. If
+ `permutation' is NULL, then the input won't be touched (no re-ordering).
+
+ This is a re-implementation of GSL's `gsl_permute' function (from its
+ `permutation/permute_source.c'). The reason we didn't use that function
+ was that it uses system-specific types (like `long' and `int') which are
+ not easily convertable to Gnuastro's width-based types. There is also a
+ separate function for each type, heavily using macros to allow a "base"
+ function to work on all the types. Thus it is hard to
+ read/understand. Since we use fixed-width types, we can easily use
+ `memcpy' and have a type-agnostic implementation (only needing the width
+ of the type).
+
+ As described in GSL's source code and manual, this implementation comes
+ from Knuth's "Art of computer programmin" book, the "Sorting and
+ Searching" chapter of Volume 3 (3rd ed). Section 5.2 Exercise 10
+ (answers), p 617. See there fore more explanations. The algorithm is a
+ little too abstract, but memory and CPU efficient.
+
+ Definition of permutations:
+
+ permute: OUT[ i ] = IN[ perm[i] ] i = 0 .. N-1
+ inverse: OUT[ perm[i] ] = IN[ i ] i = 0 .. N-1
+*/
+void
+gal_permutation_apply(gal_data_t *input, size_t *permutation)
+{
+ void *tmp;
+ size_t i, k, pk, width;
+ uint8_t *array=input->array;
+
+ /* If permutation is NULL, then it is assumed that the data doesn't need
+ to be re-ordered. */
+ if(permutation)
+ {
+ /* Necessary initializations. */
+ width=gal_data_sizeof(input->type);
+ tmp=gal_data_malloc_array(input->type, 1);
+
+ /* Do the permutation. */
+ for(i=0;i<input->size;++i)
+ {
+ k=permutation[i];
+
+ while(k>i) k=permutation[k];
+
+ if(k>=i)
+ {
+ pk = permutation[k];
+ if( pk != i )
+ {
+ memcpy(tmp, &array[i*width], width);
+
+ while(pk!=i)
+ {
+ memcpy(&array[k*width], &array[pk*width], width);
+ k = pk;
+ pk = permutation[k];
+ }
+
+ memcpy(&array[k*width], tmp, width);
+ }
+ }
+ }
+
+ /* Clean up. */
+ free(tmp);
+ }
+}
+
+
+
+
+/* Apply the inverse of given permutation on the input dataset, see
+ `gal_permutation_apply_inverse'. */
+void
+gal_permutation_apply_inverse(gal_data_t *input, size_t *permutation)
+{
+ void *tmp, *ttmp;
+ size_t i, k, pk, width;
+ uint8_t *array=input->array;
+
+ if(permutation)
+ {
+ /* Initializations */
+ width=gal_data_sizeof(input->type);
+ tmp=gal_data_malloc_array(input->type, 1);
+ ttmp=gal_data_malloc_array(input->type, 1);
+
+ /* Re-order the values. */
+ for(i=0;i<input->size;++i)
+ {
+ k=permutation[i];
+
+ while(k>i) k=permutation[k];
+
+ if(k>=i)
+ {
+ pk = permutation[k];
+
+ if(pk!=i)
+ {
+ memcpy(tmp, &array[k*width], width);
+
+ while(pk!=i)
+ {
+ memcpy(ttmp, &array[pk*width], width);
+ memcpy(&array[pk*width], tmp, width);
+ memcpy(tmp, ttmp, width);
+ k = pk;
+ pk = permutation[k];
+ }
+
+ memcpy(&array[pk*width], tmp, width);
+ }
+ }
+ }
+
+ /* Clean up. */
+ free(tmp);
+ free(ttmp);
+ }
+}
diff --git a/lib/tile.c b/lib/tile.c
index 0b28f8a..9dc2531 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -33,7 +33,9 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/blank.h>
#include <gnuastro/threads.h>
#include <gnuastro/convolve.h>
-#include <gnuastro/multidim.h>
+#include <gnuastro/dimension.h>
+#include <gnuastro/interpolate.h>
+#include <gnuastro/permutation.h>
#include "checkset.h"
@@ -64,7 +66,7 @@ gal_tile_start_coord(gal_data_t *tile, size_t *start_coord)
{
/* Calculate the coordinates of the first pixel of the tile. */
ind = gal_data_ptr_dist(block->array, tile->array, block->type);
- gal_multidim_index_to_coord(ind, ndim, block->dsize, start_coord);
+ gal_dimension_index_to_coord(ind, ndim, block->dsize, start_coord);
}
}
@@ -99,7 +101,7 @@ gal_tile_start_end_coord(gal_data_t *tile, size_t
*start_end, int rel_block)
/* Get the coordinates of the starting point relative to the allocated
block. */
- gal_multidim_index_to_coord(start_ind, ndim, block->dsize, start_end);
+ gal_dimension_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. */
@@ -111,13 +113,13 @@ gal_tile_start_end_coord(gal_data_t *tile, size_t
*start_end, int rel_block)
/* 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);
+ gal_dimension_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, end, ndim);
+ gal_dimension_add_coords(start_end, tile->dsize, end, ndim);
}
@@ -146,7 +148,7 @@ gal_tile_start_end_ind_inclusive(gal_data_t *tile,
gal_data_t *work,
/* To find the end index, we need to know the coordinates of the starting
point in the allocated block. */
- gal_multidim_index_to_coord(start_end_inc[0], ndim, block->dsize,
+ gal_dimension_index_to_coord(start_end_inc[0], ndim, block->dsize,
start_coord);
@@ -161,7 +163,8 @@ gal_tile_start_end_ind_inclusive(gal_data_t *tile,
gal_data_t *work,
/* Convert the (inclusive) ending point's coordinates into an index. */
- start_end_inc[1]=gal_multidim_coord_to_index(ndim, block->dsize, end_coord);
+ start_end_inc[1]=gal_dimension_coord_to_index(ndim, block->dsize,
+ end_coord);
/* For a check:
@@ -571,7 +574,7 @@ gal_tile_full(gal_data_t *input, size_t *regular,
dimension, then allocate the array of tiles. */
gal_tile_full_regular_first(input, regular, remainderfrac,
first, last, tsize);
- numtiles=gal_multidim_total_size(input->ndim, tsize);
+ numtiles=gal_dimension_total_size(input->ndim, tsize);
/* Allocate the necessary space for all the tiles (if necessary). */
@@ -594,7 +597,7 @@ gal_tile_full(gal_data_t *input, size_t *regular,
{
/* 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);
+ gal_dimension_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
@@ -617,7 +620,7 @@ gal_tile_full(gal_data_t *input, size_t *regular,
/* 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);
+ tind=gal_dimension_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'
@@ -787,7 +790,7 @@ gal_tile_full_two_layers(gal_data_t *input,
/* Initialize necessary values and do the channels tessellation. */
tl->numchannels = gal_tile_full(input, tl->channelsize, tl->remainderfrac,
&tl->channels, 1);
- tl->totchannels = gal_multidim_total_size(ndim, tl->numchannels);
+ tl->totchannels = gal_dimension_total_size(ndim, tl->numchannels);
/* Tile each channel. While tiling the first channel, we are also going
@@ -796,7 +799,7 @@ gal_tile_full_two_layers(gal_data_t *input,
tl->numtilesinch = gal_tile_full(tl->channels, tl->tilesize,
tl->remainderfrac, &tl->tiles,
tl->totchannels);
- tl->tottilesinch = gal_multidim_total_size(ndim, tl->numtilesinch);
+ tl->tottilesinch = gal_dimension_total_size(ndim, tl->numtilesinch);
for(i=1; i<tl->totchannels; ++i)
{
/* Set the first tile in this channel. Then use it it fill the `next'
@@ -818,137 +821,124 @@ gal_tile_full_two_layers(gal_data_t *input,
tl->numtiles = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
for(i=0;i<ndim;++i)
tl->numtiles[i] = tl->numtilesinch[i] * tl->numchannels[i];
- tl->tottiles = gal_multidim_total_size(ndim, tl->numtiles);
+ tl->tottiles = gal_dimension_total_size(ndim, tl->numtiles);
}
-/* If necessary (see below), make a new Re-ordered array that has one
- pixel/element for each tile, such that it is ready for displaying. When
- it isn't necessary, this function will just return the input data
- structure. So you should free the output only if it is different from
- the input.
+/* Usage
+ -----
- When there is only one channel OR one dimension, you can just save the
- array and everything will be fine. However, when there is more than one
- channel AND more than one dimension this won't work and the values will
- not be in the places corresponding to the original dataset. This happens
- because the tiles (and thus the value you assign to each tile) in each
- channel are taken to be a contiguous patch of memory. But when you look
- at the whole tessellation, you want the tiles to be ordered based on
- their position in the final dataset/image.
+ Make a permutation to allow the conversion of tile location in memory to
+ its location in the full input dataset and put it in the input's
+ `permutation' element. If a permutation has already been defined for the
+ tessellation, this function will not do anythin. If permutation won't be
+ necessary, then this function will just return (the permutation must
+ have been initialized to NULL).
+
+
+ Theory
+ ------
+
+ When there is only one channel OR one dimension, the tiles are allocated
+ in memory in the same order that they represent the input data. However,
+ to make channel-independent processing possible in a generic way, the
+ tiles of each channel are allocated contiguously. So, when there is more
+ than one channel AND more than one dimension, the index of the tile does
+ not correspond to its position in the grid covering the input dataset.
The example below may help clarify: assume you have a 6x6 tessellation
with two channels in the horizontal and one in the vertical. On the left
- you can see how the tiles (and their values) are allocated in memory. On
- the right is how they will be displayed if you just save it as a FITS
- file with no modification (this will not correspond to your input
- dataset).
+ you can see how the tile IDs correspond to the input dataset. NOTE how
+ `03' is on the second row, not on the first after `02'. On the right,
+ you can see how the tiles are stored in memory (and shown if you simply
+ write the array into a FITS file for example).
- Allocation map When displayed
- -------------- --------------
+ Corresponding to input In memory
+ ---------------------- --------------
15 16 17 33 34 35 30 31 32 33 34 35
12 13 14 30 31 32 24 25 26 27 28 29
09 10 11 27 28 29 18 19 20 21 22 23
- 06 07 08 24 25 26 12 13 14 15 16 17
+ 06 07 08 24 25 26 <-- 12 13 14 15 16 17
03 04 05 21 22 23 06 07 08 09 10 11
- 00 01 02 18 19 20 00 01 02 03 04 05 */
-gal_data_t *
-gal_tile_full_values_for_disp(gal_data_t *tilevalues,
- struct gal_tile_two_layer_params *tl)
+ 00 01 02 18 19 20 00 01 02 03 04 05
+
+ As a result, if your values are stored in same order as the tiles, and
+ you want them in over-all memory (for example to save as a FITS file),
+ you need to permute the values:
+
+ gal_permutation_apply(values, tl->permutation);
+
+ If you have values over-all and you want them in tile-order, you can
+ apply the inverse permutation:
+
+ gal_permutation_apply_inverse(values, tl->permutation);
+
+ Recall that this is the definition of permutation in this context:
+
+ permute: IN_ALL[ i ] = IN_MEMORY[ perm[i] ]
+ inverse: IN_ALL[ perm[i] ] = IN_MEMORY[ i ] */
+void
+gal_tile_full_permutation(struct gal_tile_two_layer_params *tl)
{
- void *start;
- gal_data_t *out, *fake;
- size_t o_inc, num_increment, start_ind;
- size_t i, ch_ind, contig, i_inc=0, ndim=tl->ndim;
- size_t *chstart=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
- size_t *s_e_inc=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
-
- /* Make sure that the input and tessellation correspond to each
- other. This is important, because the user might mistakenly give the
- input dataset as input, not the dataset that has one value per tile.*/
- if(tilevalues->ndim!=tl->ndim)
- error(EXIT_FAILURE, 0, "the input (`tilevalues') and tessellation "
- "(`tl') in `gal_tile_full_values_for_disp' have %zu and %zu "
- "dimensions respectively! They must have the same number of "
- "dimensions", tilevalues->ndim, tl->ndim);
- for(i=0; i<ndim; ++i)
- if(tilevalues->dsize[i]!=tl->numtiles[i])
- error(EXIT_FAILURE, 0, "the total number of tiles (in all channels) "
- "of dimension %zu does not correspond between the input "
- "(`tilevalues': %zu) and tessellation (`tl': %zu)", ndim-i,
- tilevalues->dsize[i], tl->numtiles[i]);
-
- /* If there is only one dimension or one channel, then just return the
- input `tilevalues'. */
- if(tilevalues->ndim==1 || tl->totchannels==1)
- return tilevalues;
-
- /* Allocate the output. */
- out=gal_data_alloc(NULL, tilevalues->type, tilevalues->ndim,
- tilevalues->dsize, tilevalues->wcs, 0,
- tilevalues->minmapsize, tilevalues->name,
- tilevalues->unit, tilevalues->comment);
-
- /* Allocate a fake tile (with a size equal to the number of tiles in one
- channel) for easy parsing with standard techniques. */
- fake=gal_data_alloc(NULL, GAL_DATA_TYPE_UINT8, tilevalues->ndim,
- tl->numtilesinch, NULL, 0, tilevalues->minmapsize,
- NULL, NULL, NULL);
-
- /* Initialize the fake array: we will set it's `block' to the output.*/
- fake->block=out;
- free(fake->array);
- fake->type=GAL_DATA_TYPE_INVALID;
-
- /* Put the values into the proper place. */
- contig=tl->numtilesinch[ndim-1];
- for(ch_ind=0; ch_ind<tl->totchannels; ++ch_ind)
+ size_t *ch_coord, *tinch_coord;
+ size_t i, p=0, t, ch, ind_in_all, ndim=tl->ndim;
+
+ /* If the permutation has already been defined for this tessellation,
+ then there is no need to do it again here. */
+ if(tl->permutation) return;
+
+ /* If there is only one channel or one dimension, return NULL. The
+ permutation functions know that the input and output indexs are the
+ same when the permutation is NULL. */
+ if( ndim==1 || tl->totchannels==1) return;
+
+ /* Allocate the space for the permutation and coordinates. */
+ ch_coord=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
+ tinch_coord=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
+ tl->permutation=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, tl->tottiles);
+
+ /* Fill in the permutation, we use the fact that the tiles are filled
+ from the first channel to the last. */
+ for(ch=0;ch<tl->totchannels;++ch)
{
- /* Set the starting pointer of this channel for the `fake' tile that
- represents it. */
- gal_multidim_index_to_coord(ch_ind, ndim, tl->numchannels, chstart);
- for(i=0;i<ndim;++i) chstart[i]*=tl->numtilesinch[i];
- start_ind=gal_multidim_coord_to_index(ndim, tl->numtiles, chstart);
- fake->array=gal_data_ptr_increment(out->array, start_ind, out->type);
-
- /* Find the starting and ending tile indexs (inclusive in the output)
- of this channel. */
- start=gal_tile_start_end_ind_inclusive(fake, out, s_e_inc);
-
- /* Start parsing the channel and putting it in the output. */
- o_inc=0;
- num_increment=1;
- while(s_e_inc[0] + o_inc <= s_e_inc[1])
+ /* Get the coordinates of this channel's first tile. */
+ gal_dimension_index_to_coord(ch, ndim, tl->numchannels, ch_coord);
+ for(i=0;i<ndim;++i) ch_coord[i] *= tl->numtilesinch[i];
+
+ /* Go over all the tiles in this channel. */
+ for(t=0;t<tl->tottilesinch;++t)
{
- /* Copy the contiguous region from the `tilevalues' array to the
- output array. Note that since they have the same type, we can
- just use `memcpy' and don't actually have to parse it.*/
- memcpy(gal_data_ptr_increment(start, o_inc, out->type),
- gal_data_ptr_increment(tilevalues->array, i_inc, out->type),
- gal_data_sizeof(out->type)*contig);
-
- /* Set the incrementation for the next contiguous patch. */
- o_inc += gal_tile_block_increment(out, fake->dsize,
- num_increment++, NULL);
- i_inc += contig;
+ /* Convert its index to coordinates and add them to the channel's
+ starting coordinates. */
+ gal_dimension_index_to_coord(t, ndim, tl->numtilesinch,
+ tinch_coord);
+ for(i=0;i<ndim;++i) tinch_coord[i] += ch_coord[i];
+
+ /* Convert the coordinates into an index. */
+ ind_in_all = gal_dimension_coord_to_index(ndim, tl->numtiles,
+ tinch_coord);
+ tl->permutation[ind_in_all] = p++;
}
}
/* Clean up and return. */
- fake->array=NULL;
- gal_data_free(fake);
- free(chstart);
- free(s_e_inc);
- return out;
+ free(tinch_coord);
+ free(ch_coord);
}
+/* Write one value for each tile into a file.
+
+ IMPORTANT: it is assumed that the values are in the same order as the
+ tiles.
+
+ tile[i] --> tilevalues[i] */
void
gal_tile_full_values_write(gal_data_t *tilevalues,
struct gal_tile_two_layer_params *tl,
@@ -956,14 +946,86 @@ gal_tile_full_values_write(gal_data_t *tilevalues,
{
gal_data_t *disp;
+
/* Make the dataset to be displayed. */
- disp = ( tl->oneelempertile
- ? gal_tile_full_values_for_disp(tilevalues, tl)
- : gal_tile_block_write_const_value(tilevalues, tl->tiles, 0) );
+ if(tl->oneelempertile)
+ {
+ if(tl->ndim>1 && tl->totchannels>1)
+ {
+ disp = gal_data_copy(tilevalues);
+ gal_permutation_apply(disp, tl->permutation);
+ }
+ else disp = tilevalues;
+ }
+ else
+ disp=gal_tile_block_write_const_value(tilevalues, tl->tiles, 0);
- /* Write the array as a file and then clean up. */
+
+ /* Write the array as a file and then clean up (if necessary). */
gal_fits_img_write(disp, filename, NULL, program_string);
- gal_data_free(disp);
+ if(disp!=tilevalues) gal_data_free(disp);
+}
+
+
+
+
+
+/* Interpolate the blank tile values while respecting the channels (if
+ requested). The output maybe a newly allocated array (if there are
+ channels), or the same array as input with interpolated values.
+
+ IMPORTANT: it is assumed that the values are in the same order as the
+ tiles.
+
+ tile[i] --> tilevalues[i] */
+void
+gal_tile_full_interpolate(gal_data_t *tilevalues,
+ struct gal_tile_two_layer_params *tl)
+{
+ size_t i;
+ void *chvstart;
+ gal_data_t *tointerp;
+
+ /* Ignore channels (if requested). */
+ if(tl->workoverch)
+ {
+ /* Prepare the permutation (if necessary/not already defined). */
+ gal_tile_full_permutation(tl);
+
+ /* Re-order values to ignore channels (if necessary). */
+ gal_permutation_apply(tilevalues, tl->permutation);
+
+ /* Interpolate the values. */
+ gal_interpolate(tilevalues);
+
+ /* Re-order values back to their original order. */
+ gal_permutation_apply_inverse(tilevalues, tl->permutation);
+ }
+ else
+ {
+ /* Interpolate each channel individually. */
+ for(i=0;i<tl->totchannels;++i)
+ {
+ /* Set the starting pointer for this channel IN THE VALUES
+ array. */
+ chvstart=gal_data_ptr_increment(tilevalues->array,
+ i*tl->tottilesinch,
+ tilevalues->type);
+
+ /* Make a cover-data-strcuture for the values in this
+ channel. */
+ tointerp=gal_data_alloc(chvstart, tilevalues->type,
+ tilevalues->ndim, tl->numtilesinch,
+ NULL, 0, 0, NULL, NULL, NULL);
+
+ /* Interpolate over the tiles in this channel. */
+ gal_interpolate(tointerp);
+
+ /* Clean up. */
+ tointerp->array=NULL;
+ gal_data_free(tointerp);
+ }
+ }
}
@@ -981,6 +1043,7 @@ gal_tile_full_free_contents(struct
gal_tile_two_layer_params *tl)
if(tl->numtiles) free(tl->numtiles);
if(tl->numtilesinch) free(tl->numtilesinch);
if(tl->tilecheckname) free(tl->tilecheckname);
+ if(tl->permutation) free(tl->permutation);
/* Free the arrays of `gal_data_c' for each tile and channel. */
if(tl->tiles) gal_data_array_free(tl->tiles, tl->tottiles, 0);
diff --git a/lib/wcs.c b/lib/wcs.c
index 8957f2e..b4c7e64 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -36,7 +36,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/wcs.h>
#include <gnuastro/tile.h>
#include <gnuastro/fits.h>
-#include <gnuastro/multidim.h>
+#include <gnuastro/dimension.h>
@@ -245,7 +245,7 @@ gal_wcs_on_tile(gal_data_t *tile)
/* Find the coordinates of the tile's starting index. */
start_ind=gal_data_ptr_dist(block->array, tile->array, block->type);
- gal_multidim_index_to_coord(start_ind, ndim, block->dsize, coord);
+ gal_dimension_index_to_coord(start_ind, ndim, block->dsize, coord);
/* Correct the copied WCS structure. Note that crpix is indexed in
the FITS/Fortran order while coord is ordered in C, it also starts
- [gnuastro-commits] master d6aa57f 081/125: Text library can write images also, besides tables, (continued)
- [gnuastro-commits] master d6aa57f 081/125: Text library can write images also, besides tables, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 67cc39d 107/125: Structure to keep tessellation information, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master b03f7f9 113/125: New inplementation of interpolation complete, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 5fa37eb 096/125: HDU modification now fully implemented in Fits program, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 38f401b 121/125: Clump S/N threshold is now calculated on Sky clumps, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 1e961f5 067/125: Starting to use gal_data_t in MakeProfiles, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e58ed29 108/125: One value per tile in multi-channel tessellation, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0e344fe 086/125: Warp now uses gal_data_t and all its fruits, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 8348f1e 119/125: NoiseChisel pseudo-detections identified on threads, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master fca5e80 117/125: First re-implementation of NoiseChisel's complete detection, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master a1ff97b 111/125: Permutations and neighbor finding algorithms,
Mohammad Akhlaghi <=
- [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