[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 81c2b7f 2/2: Dimension collapsing in Arithmeti
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 81c2b7f 2/2: Dimension collapsing in Arithmetic program and library |
|
Date: |
Wed, 4 Jul 2018 11:46:38 -0400 (EDT) |
branch: master
commit 81c2b7f9310d6293d250da7f0dcb5068cc005e72
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Dimension collapsing in Arithmetic program and library
Three new dimension collapsing operators have been added to the Arithmetic
program and the `dimension.h' library.
---
NEWS | 13 +-
bin/arithmetic/arithmetic.c | 90 +++++++++-
bin/arithmetic/arithmetic.h | 3 +
doc/gnuastro.texi | 109 ++++++++++--
lib/arithmetic.c | 139 +++++++--------
lib/dimension.c | 425 ++++++++++++++++++++++++++++++++++++++++++++
lib/gnuastro/dimension.h | 17 ++
lib/gnuastro/wcs.h | 3 +
lib/wcs.c | 170 ++++++++++++++++++
9 files changed, 883 insertions(+), 86 deletions(-)
diff --git a/NEWS b/NEWS
index 5a359f9..b4c0010 100644
--- a/NEWS
+++ b/NEWS
@@ -5,8 +5,19 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
** New features
+ Arithmetic:
+ - `collapse-sum': collapse/remove a dimension by summing over it.
+ - `collapse-mean': collapse/remove a dimension by averaging over it.
+ - `collapse-number': Number of elements included in the collapse.
+
Table:
- - `--colinfoinstdout': column information when writing to standard output.
+ --colinfoinstdout: column information when writing to standard output.
+
+ Library:
+ - gal_dimension_collapse_sum: collapse/remove a dimension by summing.
+ - gal_dimension_collapse_mean: collapse/remove a dimension by averaging.
+ - gal_dimension_collapse_number: collapse/remove a dimension by number.
+ - gal_wcs_remove_dimension: Remove a dimension in the given WCS structure.
** Removed features
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 44ab2e5..f89dbed 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -685,6 +685,82 @@ arithmetic_interpolate(struct arithmeticparams *p, char
*token)
+static void
+arithmetic_collapse(struct arithmeticparams *p, char *token, int operator)
+{
+ long dim;
+ gal_data_t *collapsed=NULL;
+
+ /* First popped operand is the dimension. */
+ gal_data_t *dimension = operands_pop(p, token);
+
+ /* The second popped operand is the desired input dataset. */
+ gal_data_t *input = operands_pop(p, token);
+
+
+ /* Small sanity check. */
+ if( dimension->ndim!=1 || dimension->size!=1)
+ error(EXIT_FAILURE, 0, "First popped operand of `collapse-*' operators "
+ "(dimension to collapse) must be a single number (single-element, "
+ "one-dimensional dataset). But it has %zu dimension(s) and %zu "
+ "element(s).", dimension->ndim, dimension->size);
+ if(dimension->type==GAL_TYPE_FLOAT32 || dimension->type==GAL_TYPE_FLOAT64)
+ error(EXIT_FAILURE, 0, "First popped operand of `collapse-*' operators "
+ "(dimension to collapse) must have an integer type, but it has "
+ "a floating point type (`%s')", gal_type_name(dimension->type,1));
+ dimension=gal_data_copy_to_new_type_free(dimension, GAL_TYPE_LONG);
+ dim=((long *)(dimension->array))[0];
+ if(dim<0 || dim==0)
+ error(EXIT_FAILURE, 0, "First popped operand of `collapse-*' operators "
+ "(dimension to collapse) must be positive (larger than zero), it "
+ "is %ld", dim);
+
+
+ /* If a WCS structure has been read, we'll need to pass it to
+ `gal_dimension_collapse', so it modifies it respectively. */
+ input->wcs=p->refdata.wcs;
+
+
+ /* Run the relevant library function. */
+ switch(operator)
+ {
+ case ARITHMETIC_OP_COLLAPSE_SUM:
+ collapsed=gal_dimension_collapse_sum(input, input->ndim-dim, NULL);
+ break;
+
+ case ARITHMETIC_OP_COLLAPSE_MEAN:
+ collapsed=gal_dimension_collapse_mean(input, input->ndim-dim, NULL);
+ break;
+
+ case ARITHMETIC_OP_COLLAPSE_NUMBER:
+ collapsed=gal_dimension_collapse_number(input, input->ndim-dim);
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
+ "problem. The operator code %d is not recognized", __func__,
+ PACKAGE_BUGREPORT, operator);
+ }
+
+
+ /* If a WCS structure existed, a modified WCS is now present in
+ `collapsed->wcs'. So we'll let the freeing of `input' free the old
+ `p->refdata.wcs' structure and we'll put the new one there, then we'll
+ set `collapsed->wcs' to `NULL', so the new one isn't freed. */
+ p->refdata.wcs = collapsed->wcs;
+ collapsed->wcs = NULL;
+
+
+ /* Clean up and add the collapsed dataset to the top of the operands. */
+ gal_data_free(input);
+ gal_data_free(dimension);
+ operands_add(p, NULL, collapsed);
+}
+
+
+
+
+
@@ -854,7 +930,7 @@ reversepolish(struct arithmeticparams *p)
else if (!strcmp(token->v, "float64"))
{ op=GAL_ARITHMETIC_OP_TO_FLOAT64; nop=1; }
- /* Filters. */
+ /* Library wrappers. */
else if (!strcmp(token->v, "filter-mean"))
{ op=ARITHMETIC_OP_FILTER_MEAN; nop=0; }
else if (!strcmp(token->v, "filter-median"))
@@ -873,6 +949,12 @@ reversepolish(struct arithmeticparams *p)
{ op=ARITHMETIC_OP_INVERT; nop=0; }
else if (!strcmp(token->v, "interpolate-medianngb"))
{ op=ARITHMETIC_OP_INTERPOLATE_MEDIANNGB; nop=0; }
+ else if (!strcmp(token->v, "collapse-sum"))
+ { op=ARITHMETIC_OP_COLLAPSE_SUM; nop=0; }
+ else if (!strcmp(token->v, "collapse-mean"))
+ { op=ARITHMETIC_OP_COLLAPSE_MEAN; nop=0; }
+ else if (!strcmp(token->v, "collapse-number"))
+ { op=ARITHMETIC_OP_COLLAPSE_NUMBER; nop=0; }
/* Finished checks with known operators */
@@ -964,6 +1046,12 @@ reversepolish(struct arithmeticparams *p)
arithmetic_interpolate(p, token->v);
break;
+ case ARITHMETIC_OP_COLLAPSE_SUM:
+ case ARITHMETIC_OP_COLLAPSE_MEAN:
+ case ARITHMETIC_OP_COLLAPSE_NUMBER:
+ arithmetic_collapse(p, token->v, op);
+ break;
+
default:
error(EXIT_FAILURE, 0, "%s: a bug! please contact us at "
"%s to fix the problem. The code %d is not "
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index e5eca69..4144006 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -40,6 +40,9 @@ enum arithmetic_prog_operators
ARITHMETIC_OP_CONNECTED_COMPONENTS,
ARITHMETIC_OP_INVERT,
ARITHMETIC_OP_INTERPOLATE_MEDIANNGB,
+ ARITHMETIC_OP_COLLAPSE_SUM,
+ ARITHMETIC_OP_COLLAPSE_MEAN,
+ ARITHMETIC_OP_COLLAPSE_NUMBER,
};
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 503d6e8..fa3d408 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -5010,7 +5010,7 @@ line.
@cindex Ubuntu
@cindex @command{apt-get}
@cindex Advanced Packaging Tool (APT, Debian)
-Debian-based GNU/Linux
address@hidden://en.wikipedia.org/wiki/Debian,Debian}-based GNU/Linux
address@hidden@url{https://en.wikipedia.org/wiki/List_of_Linux_distributions#Debian-based}}
(for example Ubuntu or its derivatives) are arguably the largest, and most
used, class of GNU/Linux distributions. All such GNU/Linux distributions
@@ -5028,12 +5028,11 @@ $ sudo apt-get install ghostscript libtool-bin
libjpeg-dev \
@cindex Homebrew
@cindex MacPorts
@cindex @command{brew}
-macOS is the operating system used on Apple devices
-(@url{https://en.wikipedia.org/wiki/MacOS}). macOS does not come with a
-package manager pre-installed, but several widely used, third-party package
-managers exist, such as Homebrew or MacPorts. Both are free
-software. Currently we have only tested Gnuastro's installation with
-Homebrew as described below.
address@hidden://en.wikipedia.org/wiki/MacOS,macOS} is the operating system
+used on Apple devices. macOS does not come with a package manager
+pre-installed, but several widely used, third-party package managers exist,
+such as Homebrew or MacPorts. Both are free software. Currently we have
+only tested Gnuastro's installation with Homebrew as described below.
If not already installed, first obtain Homebrew by following the
instructions at @url{https://brew.sh}. Homebrew manages packages in
@@ -5052,13 +5051,12 @@ $ brew install wcslib
@item @command{pacman} (Arch Linux)
@cindex Arch Linux
@cindex @command{pacman}
-Arch Linux is a smaller GNU/Linux distribution. As discussed in
address@hidden://en.wikipedia.org/wiki/Arch_Linux, Wikipedia}, it follows ``the
-KISS principle ("keep it simple, stupid") as the general guideline, and
-focuses on elegance, code correctness, minimalism and simplicity, and
-expects the user to be willing to make some effort to understand the
-system's operation''. Arch Linux uses ``Package manager'' (Pacman) to
-manage its packages/components.
address@hidden://en.wikipedia.org/wiki/Arch_Linux,Arch Linux} is a smaller
+GNU/Linux distribution, which follows the KISS principle (``keep it simple,
+stupid'') as a general guideline. It ``focuses on elegance, code
+correctness, minimalism and simplicity, and expects the user to be willing
+to make some effort to understand the system's operation''. Arch Linux uses
+``Package manager'' (Pacman) to manage its packages/components.
@example
$ sudo pacman -S ghostscript libtool libjpeg libtiff libgit2 \
\
@@ -11230,6 +11228,55 @@ non-blank neighbors used to calculate the median is
given by the first
popped operand. Note that the distance of the nearest non-blank neighbors
is irrelevant in this interpolation.
address@hidden collapse-sum
+Collapse the given dataset (second popped operand), by summing all elements
+along the first popped operand (a dimension in FITS standard: counting from
+one, from fastest dimension). The returned dataset has one dimension less
+compared to the input.
+
+The output will have a double-precision floating point type irrespective of
+the input dataset's type. Doing the operation in double-precision (64-bit)
+floating point will help the collapse (summation) be affected less by
+floating point errors. But afterwards, single-precision floating points are
+usually enough in real (noisy) datasets. So depending on the type of the
+input and its nature, it is recommended to use one of the type conversion
+operators on the returned dataset.
+
address@hidden World Coordinate System (WCS)
+If any WCS is present, the returned dataset will also lack the respective
+dimension in its WCS matrix. Therefore, when the WCS is important for later
+processing, be sure that the input is aligned with the respective axises: all
non-diagonal elements in the WCS matrix are zero.
+
address@hidden IFU
address@hidden Data cubes
address@hidden 3D data-cubes
address@hidden Cubes (3D data)
address@hidden Narrow-band image
address@hidden Integral field unit (IFU)
+One common application of this operator is the creation of pseudo
+broad-band or narrow-band 2D images from 3D data cubes. For example
+integral field unit (IFU) data products that have two spatial dimensions
+(first two FITS dimensions) and one spectral dimension (third FITS
+dimension). The command below will collapse the whole third dimension into
+a 2D array the size of the first two dimensions, and then convert the
+output to single-precision floating point (as discussed above).
+
address@hidden
+$ astarithmetic cube.fits 3 collapse-sum float32
address@hidden example
+
address@hidden collapse-mean
+Similar to @option{collapse-sum}, but the returned dataset will be the mean
+value along the collapsed dimension, not the sum.
+
address@hidden collapse-number
+Similar to @option{collapse-sum}, but the returned dataset will be the
+number of non-blank values along the collapsed dimension. The output will
+have a 32-bit signed integer type. If the input dataset doesn't have blank
+values, all the elements in the returned dataset will have a single value
+(the length of the collapsed dimension). Therefore this is mostly relevant
+when there are blank values in the dataset.
+
@item erode
@cindex Erosion
Erode the foreground pixels (with value @code{1}) of the input dataset
@@ -22771,6 +22818,40 @@ the two coordinates @code{a} and @code{b} (each an
array of @code{ndim}
elements).
@end deftypefun
address@hidden {gal_data_t *} gal_dimension_collapse_sum (gal_data_t
@code{*in}, size_t @code{c_dim}, gal_data_t @code{*weight})
+Collapse the input dataset (@code{in}) along the given dimension
+(@code{c_dim}, in C definition: starting from zero, from the slowest
+dimension), by summing all elements in that direction. If
address@hidden, it must be a single-dimensional array, with the same
+size as the dimension to be collapsed. The respective weight will be
+multiplied to each element during the collapse.
+
+For generality, the returned dataset will have a @code{GAL_TYPE_FLOAT64}
+type. See @ref{Copying datasets} for converting the returned dataset to a
+desired type. Also, for more on the application of this function, see the
+Arithmetic program's @option{collapse-sum} operator (which uses this
+function) in @ref{Arithmetic operators}.
address@hidden deftypefun
+
address@hidden {gal_data_t *} gal_dimension_collapse_mean (gal_data_t
@code{*in}, size_t @code{c_dim}, gal_data_t @code{*weight})
+Similar to @code{gal_dimension_collapse_sum} (above), but the collapse will
+be done by calculating the mean along the requested dimension, not summing
+over it.
address@hidden deftypefun
+
address@hidden {gal_data_t *} gal_dimension_collapse_number (gal_data_t
@code{*in}, size_t @code{c_dim})
+Collapse the input dataset (@code{in}) along the given dimension
+(@code{c_dim}, in C definition: starting from zero, from the slowest
+dimension), by counting how many non-blank elements there are along that
+dimension.
+
+For generality, the returned dataset will have a @code{GAL_TYPE_INT32}
+type. See @ref{Copying datasets} for converting the returned dataset to a
+desired type. Also, for more on the application of this function, see the
+Arithmetic program's @option{collapse-number} operator (which uses this
+function) in @ref{Arithmetic operators}.
address@hidden deftypefun
+
@deffn {Function-like macro} GAL_DIMENSION_NEIGHBOR_OP (@code{index},
@code{ndim}, @code{dsize}, @code{connectivity}, @code{dinc}, @code{operation})
Parse the neighbors of the element located at @code{index} and do the
requested operation on them. This is defined as a macro to allow easy
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 6f31f51..3e84f95 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -183,16 +183,16 @@ arithmetic_not(gal_data_t *data, int flags)
/* Go over the pixels and set the output values. */
switch(data->type)
{
- case GAL_TYPE_UINT8: TYPE_CASE_FOR_NOT(uint8_t); break;
- case GAL_TYPE_INT8: TYPE_CASE_FOR_NOT(int8_t); break;
- case GAL_TYPE_UINT16: TYPE_CASE_FOR_NOT(uint16_t); break;
- case GAL_TYPE_INT16: TYPE_CASE_FOR_NOT(int16_t); break;
- case GAL_TYPE_UINT32: TYPE_CASE_FOR_NOT(uint32_t); break;
- case GAL_TYPE_INT32: TYPE_CASE_FOR_NOT(int32_t); break;
- case GAL_TYPE_UINT64: TYPE_CASE_FOR_NOT(uint64_t); break;
- case GAL_TYPE_INT64: TYPE_CASE_FOR_NOT(int64_t); break;
- case GAL_TYPE_FLOAT32: TYPE_CASE_FOR_NOT(float); break;
- case GAL_TYPE_FLOAT64: TYPE_CASE_FOR_NOT(double); break;
+ case GAL_TYPE_UINT8: TYPE_CASE_FOR_NOT( uint8_t ); break;
+ case GAL_TYPE_INT8: TYPE_CASE_FOR_NOT( int8_t ); break;
+ case GAL_TYPE_UINT16: TYPE_CASE_FOR_NOT( uint16_t ); break;
+ case GAL_TYPE_INT16: TYPE_CASE_FOR_NOT( int16_t ); break;
+ case GAL_TYPE_UINT32: TYPE_CASE_FOR_NOT( uint32_t ); break;
+ case GAL_TYPE_INT32: TYPE_CASE_FOR_NOT( int32_t ); break;
+ case GAL_TYPE_UINT64: TYPE_CASE_FOR_NOT( uint64_t ); break;
+ case GAL_TYPE_INT64: TYPE_CASE_FOR_NOT( int64_t ); break;
+ case GAL_TYPE_FLOAT32: TYPE_CASE_FOR_NOT( float ); break;
+ case GAL_TYPE_FLOAT64: TYPE_CASE_FOR_NOT( double ); break;
case GAL_TYPE_BIT:
error(EXIT_FAILURE, 0, "%s: bit datatypes are not yet supported, "
@@ -1330,63 +1330,63 @@ gal_arithmetic_operator_string(int operator)
{
switch(operator)
{
- case GAL_ARITHMETIC_OP_PLUS: return "+";
- case GAL_ARITHMETIC_OP_MINUS: return "-";
- case GAL_ARITHMETIC_OP_MULTIPLY: return "*";
- case GAL_ARITHMETIC_OP_DIVIDE: return "/";
- case GAL_ARITHMETIC_OP_MODULO: return "%";
-
- case GAL_ARITHMETIC_OP_LT: return "<";
- case GAL_ARITHMETIC_OP_LE: return "<=";
- case GAL_ARITHMETIC_OP_GT: return ">";
- case GAL_ARITHMETIC_OP_GE: return ">=";
- case GAL_ARITHMETIC_OP_EQ: return "==";
- case GAL_ARITHMETIC_OP_NE: return "!=";
- case GAL_ARITHMETIC_OP_AND: return "and";
- case GAL_ARITHMETIC_OP_OR: return "or";
- case GAL_ARITHMETIC_OP_NOT: return "not";
- case GAL_ARITHMETIC_OP_ISBLANK: return "isblank";
- case GAL_ARITHMETIC_OP_WHERE: return "where";
-
- case GAL_ARITHMETIC_OP_BITAND: return "bitand";
- case GAL_ARITHMETIC_OP_BITOR: return "bitor";
- case GAL_ARITHMETIC_OP_BITXOR: return "bitxor";
- case GAL_ARITHMETIC_OP_BITLSH: return "lshift";
- case GAL_ARITHMETIC_OP_BITRSH: return "rshift";
- case GAL_ARITHMETIC_OP_BITNOT: return "bitnot";
-
- case GAL_ARITHMETIC_OP_ABS: return "abs";
- case GAL_ARITHMETIC_OP_POW: return "pow";
- case GAL_ARITHMETIC_OP_SQRT: return "sqrt";
- case GAL_ARITHMETIC_OP_LOG: return "log";
- case GAL_ARITHMETIC_OP_LOG10: return "log10";
-
- case GAL_ARITHMETIC_OP_MINVAL: return "minvalue";
- case GAL_ARITHMETIC_OP_MAXVAL: return "maxvalue";
- case GAL_ARITHMETIC_OP_NUMVAL: return "numvalue";
- case GAL_ARITHMETIC_OP_SUMVAL: return "sumvalue";
- case GAL_ARITHMETIC_OP_MEANVAL: return "meanvalue";
- case GAL_ARITHMETIC_OP_STDVAL: return "stdvalue";
- case GAL_ARITHMETIC_OP_MEDIANVAL: return "medianvalue";
-
- case GAL_ARITHMETIC_OP_MIN: return "min";
- case GAL_ARITHMETIC_OP_MAX: return "max";
- case GAL_ARITHMETIC_OP_NUM: return "num";
- case GAL_ARITHMETIC_OP_SUM: return "sum";
- case GAL_ARITHMETIC_OP_MEAN: return "mean";
- case GAL_ARITHMETIC_OP_STD: return "std";
- case GAL_ARITHMETIC_OP_MEDIAN: return "median";
-
- case GAL_ARITHMETIC_OP_TO_UINT8: return "uchar";
- case GAL_ARITHMETIC_OP_TO_INT8: return "char";
- case GAL_ARITHMETIC_OP_TO_UINT16: return "ushort";
- case GAL_ARITHMETIC_OP_TO_INT16: return "short";
- case GAL_ARITHMETIC_OP_TO_UINT32: return "uint";
- case GAL_ARITHMETIC_OP_TO_INT32: return "int";
- case GAL_ARITHMETIC_OP_TO_UINT64: return "ulong";
- case GAL_ARITHMETIC_OP_TO_INT64: return "long";
- case GAL_ARITHMETIC_OP_TO_FLOAT32: return "float32";
- case GAL_ARITHMETIC_OP_TO_FLOAT64: return "float64";
+ case GAL_ARITHMETIC_OP_PLUS: return "+";
+ case GAL_ARITHMETIC_OP_MINUS: return "-";
+ case GAL_ARITHMETIC_OP_MULTIPLY: return "*";
+ case GAL_ARITHMETIC_OP_DIVIDE: return "/";
+ case GAL_ARITHMETIC_OP_MODULO: return "%";
+
+ case GAL_ARITHMETIC_OP_LT: return "<";
+ case GAL_ARITHMETIC_OP_LE: return "<=";
+ case GAL_ARITHMETIC_OP_GT: return ">";
+ case GAL_ARITHMETIC_OP_GE: return ">=";
+ case GAL_ARITHMETIC_OP_EQ: return "==";
+ case GAL_ARITHMETIC_OP_NE: return "!=";
+ case GAL_ARITHMETIC_OP_AND: return "and";
+ case GAL_ARITHMETIC_OP_OR: return "or";
+ case GAL_ARITHMETIC_OP_NOT: return "not";
+ case GAL_ARITHMETIC_OP_ISBLANK: return "isblank";
+ case GAL_ARITHMETIC_OP_WHERE: return "where";
+
+ case GAL_ARITHMETIC_OP_BITAND: return "bitand";
+ case GAL_ARITHMETIC_OP_BITOR: return "bitor";
+ case GAL_ARITHMETIC_OP_BITXOR: return "bitxor";
+ case GAL_ARITHMETIC_OP_BITLSH: return "lshift";
+ case GAL_ARITHMETIC_OP_BITRSH: return "rshift";
+ case GAL_ARITHMETIC_OP_BITNOT: return "bitnot";
+
+ case GAL_ARITHMETIC_OP_ABS: return "abs";
+ case GAL_ARITHMETIC_OP_POW: return "pow";
+ case GAL_ARITHMETIC_OP_SQRT: return "sqrt";
+ case GAL_ARITHMETIC_OP_LOG: return "log";
+ case GAL_ARITHMETIC_OP_LOG10: return "log10";
+
+ case GAL_ARITHMETIC_OP_MINVAL: return "minvalue";
+ case GAL_ARITHMETIC_OP_MAXVAL: return "maxvalue";
+ case GAL_ARITHMETIC_OP_NUMVAL: return "numvalue";
+ case GAL_ARITHMETIC_OP_SUMVAL: return "sumvalue";
+ case GAL_ARITHMETIC_OP_MEANVAL: return "meanvalue";
+ case GAL_ARITHMETIC_OP_STDVAL: return "stdvalue";
+ case GAL_ARITHMETIC_OP_MEDIANVAL: return "medianvalue";
+
+ case GAL_ARITHMETIC_OP_MIN: return "min";
+ case GAL_ARITHMETIC_OP_MAX: return "max";
+ case GAL_ARITHMETIC_OP_NUM: return "num";
+ case GAL_ARITHMETIC_OP_SUM: return "sum";
+ case GAL_ARITHMETIC_OP_MEAN: return "mean";
+ case GAL_ARITHMETIC_OP_STD: return "std";
+ case GAL_ARITHMETIC_OP_MEDIAN: return "median";
+
+ case GAL_ARITHMETIC_OP_TO_UINT8: return "uchar";
+ case GAL_ARITHMETIC_OP_TO_INT8: return "char";
+ case GAL_ARITHMETIC_OP_TO_UINT16: return "ushort";
+ case GAL_ARITHMETIC_OP_TO_INT16: return "short";
+ case GAL_ARITHMETIC_OP_TO_UINT32: return "uint";
+ case GAL_ARITHMETIC_OP_TO_INT32: return "int";
+ case GAL_ARITHMETIC_OP_TO_UINT64: return "ulong";
+ case GAL_ARITHMETIC_OP_TO_INT64: return "long";
+ case GAL_ARITHMETIC_OP_TO_FLOAT32: return "float32";
+ case GAL_ARITHMETIC_OP_TO_FLOAT64: return "float64";
default:
error(EXIT_FAILURE, 0, "%s: operator code %d not recognized",
@@ -1412,7 +1412,7 @@ gal_arithmetic(int operator, int flags, ...)
/* Prepare the variable arguments (starting after the flags argument). */
va_start(va, flags);
- /* Depending on the operator do the job: */
+ /* Depending on the operator, do the job: */
switch(operator)
{
@@ -1474,13 +1474,12 @@ gal_arithmetic(int operator, int flags, ...)
out=arithmetic_from_statistics(operator, flags, d1);
break;
- /* Absolute operator */
+ /* Absolute operator. */
case GAL_ARITHMETIC_OP_ABS:
d1 = va_arg(va, gal_data_t *);
out=arithmetic_abs(flags, d1);
break;
-
/* Multi-operand operators */
case GAL_ARITHMETIC_OP_MIN:
case GAL_ARITHMETIC_OP_MAX:
diff --git a/lib/dimension.c b/lib/dimension.c
index 3d8d8e2..ee76413 100644
--- a/lib/dimension.c
+++ b/lib/dimension.c
@@ -28,6 +28,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <error.h>
#include <stdlib.h>
+#include <gnuastro/wcs.h>
#include <gnuastro/pointer.h>
#include <gnuastro/dimension.h>
@@ -265,3 +266,427 @@ gal_dimension_dist_manhattan(size_t *a, size_t *b, size_t
ndim)
for(i=0;i<ndim;++i) out += (a[i] > b[i]) ? (a[i]-b[i]) : (b[i]-a[i]);
return out;
}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/************************************************************************/
+/******************** Collapsing a dimension **********************/
+/************************************************************************/
+enum dimension_collapse_operation
+{
+ DIMENSION_COLLAPSE_INVALID, /* ==0 by C standard. */
+
+ DIMENSION_COLLAPSE_SUM,
+ DIMENSION_COLLAPSE_MEAN,
+ DIMENSION_COLLAPSE_NUMBER,
+};
+
+
+
+
+
+static gal_data_t *
+dimension_collapse_sanity_check(gal_data_t *in, gal_data_t *weight,
+ size_t c_dim, int hasblank, size_t *cnum,
+ double **warr)
+{
+ gal_data_t *wht=NULL;
+
+ /* The requested dimension to collapse cannot be larger than the input's
+ number of dimensions. */
+ if( c_dim > (in->ndim-1) )
+ error(EXIT_FAILURE, 0, "%s: the input has %zu dimensions, but you have "
+ "asked to collapse dimension %zu", __func__, in->ndim, c_dim);
+
+ /* If there is no blank value, there is no point in calculating the
+ number of points in each collapsed dataset (when necessary). In that
+ case, `cnum!=0'. */
+ if(hasblank==0)
+ *cnum=in->dsize[c_dim];
+
+ /* Weight sanity checks */
+ if(weight)
+ {
+ if( weight->ndim!=1 )
+ error(EXIT_FAILURE, 0, "%s: the weight dataset has %zu dimensions, "
+ "it must be one-dimensional", __func__, weight->ndim);
+ if( in->dsize[c_dim]!=weight->size )
+ error(EXIT_FAILURE, 0, "%s: the weight dataset has %zu elements, "
+ "but the input dataset has %zu elements in dimension %zu",
+ __func__, weight->size, in->dsize[c_dim], c_dim);
+ wht = ( weight->type == GAL_TYPE_FLOAT64
+ ? weight
+ : gal_data_copy_to_new_type(weight, GAL_TYPE_FLOAT64) );
+ *warr = wht->array;
+ }
+
+ /* Return the weight data structure. */
+ return wht;
+}
+
+
+
+
+/* Set the collapsed output sizes. */
+static void
+dimension_collapse_sizes(gal_data_t *in, size_t c_dim, size_t *outndim,
+ size_t *outdsize)
+{
+ size_t i, a=0;
+
+ if(in->ndim==1)
+ *outndim=outdsize[0]=1;
+ else
+ {
+ *outndim=in->ndim-1;
+ for(i=0;i<in->ndim;++i)
+ if(i!=c_dim) outdsize[a++]=in->dsize[i];
+ }
+}
+
+
+
+
+
+/* Depending on the operator, write the result into the output. */
+#define COLLAPSE_WRITE(OIND,IIND) { \
+ /* We need the sum when number operator is requested. */ \
+ if(farr) farr[ OIND ] += (warr ? warr[w] : 1) * inarr[ IIND ]; \
+ \
+ /* We don't need the number when the sum operator is requested. */ \
+ if(iarr) \
+ { \
+ if(num->type==GAL_TYPE_UINT8) iarr[ OIND ] = 1; \
+ else ++iarr[ OIND ]; \
+ } \
+ \
+ /* If the sum of weights for is needed, add it. */ \
+ if(wsumarr) wsumarr[ OIND ] += warr[w]; \
+ }
+
+
+
+/* Deal properly with blanks. */
+#define COLLAPSE_CHECKBLANK(OIND,IIND) { \
+ if(hasblank) \
+ { \
+ if(B==B) /* An integer type: blank can be checked with `=='. */ \
+ { \
+ if( inarr[IIND] != B ) COLLAPSE_WRITE(OIND,IIND); \
+ } \
+ else /* A floating point type where NAN != NAN. */ \
+ { \
+ if( inarr[IIND] == inarr[IIND] ) COLLAPSE_WRITE(OIND,IIND); \
+ } \
+ } \
+ else COLLAPSE_WRITE(OIND,IIND); \
+ }
+
+
+
+#define COLLAPSE_DIM(IT) { \
+ IT B, *inarr=in->array; \
+ if(hasblank) gal_blank_write(&B, in->type); \
+ switch(in->ndim) \
+ { \
+ /* 1D input dataset. */ \
+ case 1: \
+ for(i=0;i<in->dsize[0];++i) \
+ { \
+ if(weight) w=i; \
+ COLLAPSE_CHECKBLANK(0,i); \
+ } \
+ break; \
+ \
+ /* 2D input dataset. */ \
+ case 2: \
+ for(i=0;i<in->dsize[0];++i) \
+ for(j=0;j<in->dsize[1];++j) \
+ { \
+ /* In a more easy to understand format: \
+ dim==0 --> a=j; \
+ dim==1 --> a=i; */ \
+ a = c_dim==0 ? j : i; \
+ if(weight) w = c_dim == 0 ? i : j; \
+ COLLAPSE_CHECKBLANK(a, i*in->dsize[1] + j); \
+ } \
+ break; \
+ \
+ /* 3D input dataset. */ \
+ case 3: \
+ slice=in->dsize[1]*in->dsize[2]; \
+ for(i=0;i<in->dsize[0];++i) \
+ for(j=0;j<in->dsize[1];++j) \
+ for(k=0;k<in->dsize[2];++k) \
+ { \
+ /* In a more easy to understand format: \
+ dim==0 --> a=j; b=k; \
+ dim==1 --> a=i; b=k; \
+ dim==2 --> a=i; b=j; */ \
+ a = c_dim==0 ? j : i; \
+ b = c_dim==2 ? j : k; \
+ if(weight) w = c_dim==0 ? i : (c_dim==1 ? j : k); \
+ COLLAPSE_CHECKBLANK(a*outdsize[1]+b, \
+ i*slice + j*in->dsize[2] + k); \
+ } \
+ break; \
+ \
+ /* Input dataset's dimensionality not yet supported. */ \
+ default: \
+ error(EXIT_FAILURE, 0, "%s: %zu-dimensional datasets not yet " \
+ "supported, please contact us at %s to add this feature", \
+ __func__, in->ndim, PACKAGE_BUGREPORT); \
+ } \
+ }
+
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_sum(gal_data_t *in, size_t c_dim, gal_data_t *weight)
+{
+ double *wsumarr=NULL;
+ int8_t *ii, *iarr=NULL;
+ size_t a, b, i, j, k, w, cnum=0;
+ size_t outdsize[10], slice, outndim;
+ int hasblank=gal_blank_present(in, 0);
+ double *dd, *df, *warr=NULL, *farr=NULL;
+ gal_data_t *wht=NULL, *sum=NULL, *num=NULL;
+
+ /* Basic sanity checks. */
+ wht=dimension_collapse_sanity_check(in, weight, c_dim, hasblank,
+ &cnum, &warr);
+
+ /* Set the size of the collapsed output. */
+ dimension_collapse_sizes(in, c_dim, &outndim, outdsize);
+
+ /* Allocate the sum (output) dataset. */
+ sum=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, outndim, outdsize, in->wcs,
+ 1, in->minmapsize, NULL, NULL, NULL);
+
+ /* The number dataset (when there are blank values).*/
+ if(hasblank)
+ num=gal_data_alloc(NULL, GAL_TYPE_INT8, outndim, outdsize, NULL,
+ 1, in->minmapsize, NULL, NULL, NULL);
+
+ /* Set the array pointers. */
+ if(sum) farr=sum->array;
+ if(num) iarr=num->array;
+
+ /* Parse the dataset. */
+ switch(in->type)
+ {
+ case GAL_TYPE_UINT8: COLLAPSE_DIM( uint8_t ); break;
+ case GAL_TYPE_INT8: COLLAPSE_DIM( int8_t ); break;
+ case GAL_TYPE_UINT16: COLLAPSE_DIM( uint16_t ); break;
+ case GAL_TYPE_INT16: COLLAPSE_DIM( int16_t ); break;
+ case GAL_TYPE_UINT32: COLLAPSE_DIM( uint32_t ); break;
+ case GAL_TYPE_INT32: COLLAPSE_DIM( int32_t ); break;
+ case GAL_TYPE_UINT64: COLLAPSE_DIM( uint64_t ); break;
+ case GAL_TYPE_INT64: COLLAPSE_DIM( int64_t ); break;
+ case GAL_TYPE_FLOAT32: COLLAPSE_DIM( float ); break;
+ case GAL_TYPE_FLOAT64: COLLAPSE_DIM( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: type value (%d) not recognized",
+ __func__, in->type);
+ }
+
+ /* If `num' is zero on any element, set its sum to NaN. */
+ if(num)
+ {
+ ii = num->array;
+ df = (dd=sum->array) + sum->size;
+ do if(*ii++==0) *dd=NAN; while(++dd<df);
+ }
+
+ /* Remove the respective dimension in the WCS structure also (if any
+ exists). Note that `sum->ndim' has already been changed. So we'll use
+ `in->wcs'. */
+ gal_wcs_remove_dimension(sum->wcs, in->ndim-c_dim);
+
+ /* Clean up and return. */
+ if(wht!=weight) gal_data_free(wht);
+ if(num) gal_data_free(num);
+ return sum;
+}
+
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_mean(gal_data_t *in, size_t c_dim,
+ gal_data_t *weight)
+{
+ double wsum=NAN;
+ double *wsumarr=NULL;
+ int8_t *ii, *iarr=NULL;
+ size_t outdsize[10], slice, outndim;
+ int hasblank=gal_blank_present(in, 0);
+ size_t a, b, i, j, k, w, cnum=0;
+ gal_data_t *wht=NULL, *sum=NULL, *num=NULL;
+ double *dd, *dw, *df, *warr=NULL, *farr=NULL;
+
+ /* Basic sanity checks. */
+ wht=dimension_collapse_sanity_check(in, weight, c_dim, hasblank,
+ &cnum, &warr);
+
+ /* Set the size of the collapsed output. */
+ dimension_collapse_sizes(in, c_dim, &outndim, outdsize);
+
+ /* The sum array. */
+ sum=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, outndim, outdsize, in->wcs,
+ 1, in->minmapsize, NULL, NULL, NULL);
+
+ /* If a weighted mean is requested. */
+ if( weight )
+ {
+ /* There are blank values, so we'll need to keep the sums of the
+ weights for each collapsed dimension */
+ if( hasblank )
+ wsumarr=gal_pointer_allocate(GAL_TYPE_FLOAT64, sum->size, 1,
+ __func__, "wsumarr");
+
+ /* There aren't any blank values, so one summation over the
+ weights is enough to calculate the weighted mean. */
+ else
+ {
+ wsum=0.0f;
+ df=(dd=weight->array)+weight->size;
+ do wsum += *dd++; while(dd<df);
+ }
+ }
+ /* No weight is given, so we'll need the number of elements. */
+ else if( hasblank )
+ num=gal_data_alloc(NULL, GAL_TYPE_INT32, outndim, outdsize, NULL,
+ 1, in->minmapsize, NULL, NULL, NULL);
+
+ /* Set the array pointers. */
+ if(sum) farr=sum->array;
+ if(num) iarr=num->array;
+
+ /* Parse the dataset. */
+ switch(in->type)
+ {
+ case GAL_TYPE_UINT8: COLLAPSE_DIM( uint8_t ); break;
+ case GAL_TYPE_INT8: COLLAPSE_DIM( int8_t ); break;
+ case GAL_TYPE_UINT16: COLLAPSE_DIM( uint16_t ); break;
+ case GAL_TYPE_INT16: COLLAPSE_DIM( int16_t ); break;
+ case GAL_TYPE_UINT32: COLLAPSE_DIM( uint32_t ); break;
+ case GAL_TYPE_INT32: COLLAPSE_DIM( int32_t ); break;
+ case GAL_TYPE_UINT64: COLLAPSE_DIM( uint64_t ); break;
+ case GAL_TYPE_INT64: COLLAPSE_DIM( int64_t ); break;
+ case GAL_TYPE_FLOAT32: COLLAPSE_DIM( float ); break;
+ case GAL_TYPE_FLOAT64: COLLAPSE_DIM( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: type value (%d) not recognized",
+ __func__, in->type);
+ }
+
+ /* If `num' is zero on any element, set its sum to NaN. */
+ if(num)
+ {
+ ii = num->array;
+ df = (dd=sum->array) + sum->size;
+ do if(*ii++==0) *dd=NAN; while(++dd<df);
+ }
+
+ /* Divide the sum by the number. */
+ df = (dd=sum->array) + sum->size;
+ if(weight)
+ {
+ if(hasblank) { dw=wsumarr; do *dd /= *dw++; while(++dd<df); }
+ else do *dd /= wsum; while(++dd<df);
+ }
+ else
+ if(num) { ii = num->array; do *dd /= *ii++; while(++dd<df); }
+ else do *dd /= cnum; while(++dd<df);
+
+ /* Correct the WCS, clean up and return. */
+ gal_wcs_remove_dimension(sum->wcs, in->ndim-c_dim);
+ if(wht!=weight) gal_data_free(wht);
+ if(wsumarr) free(wsumarr);
+ gal_data_free(num);
+ return sum;
+}
+
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_number(gal_data_t *in, size_t c_dim)
+{
+ double *wsumarr=NULL;
+ double *warr=NULL, *farr=NULL;
+ int32_t *ii, *iif, *iarr=NULL;
+ size_t a, b, i, j, k, w, cnum=0;
+ size_t outdsize[10], slice, outndim;
+ int hasblank=gal_blank_present(in, 0);
+ gal_data_t *weight=NULL, *wht=NULL, *num=NULL;
+
+ /* Basic sanity checks. */
+ wht=dimension_collapse_sanity_check(in, weight, c_dim, hasblank,
+ &cnum, &warr);
+
+ /* Set the size of the collapsed output. */
+ dimension_collapse_sizes(in, c_dim, &outndim, outdsize);
+
+ /* The number dataset (when there are blank values).*/
+ num=gal_data_alloc(NULL, GAL_TYPE_INT32, outndim, outdsize, in->wcs,
+ 1, in->minmapsize, NULL, NULL, NULL);
+
+ /* Set the array pointers. */
+ iarr=num->array;
+
+ /* Parse the input dataset (if necessary). */
+ if(hasblank)
+ switch(in->type)
+ {
+ case GAL_TYPE_UINT8: COLLAPSE_DIM( uint8_t ); break;
+ case GAL_TYPE_INT8: COLLAPSE_DIM( int8_t ); break;
+ case GAL_TYPE_UINT16: COLLAPSE_DIM( uint16_t ); break;
+ case GAL_TYPE_INT16: COLLAPSE_DIM( int16_t ); break;
+ case GAL_TYPE_UINT32: COLLAPSE_DIM( uint32_t ); break;
+ case GAL_TYPE_INT32: COLLAPSE_DIM( int32_t ); break;
+ case GAL_TYPE_UINT64: COLLAPSE_DIM( uint64_t ); break;
+ case GAL_TYPE_INT64: COLLAPSE_DIM( int64_t ); break;
+ case GAL_TYPE_FLOAT32: COLLAPSE_DIM( float ); break;
+ case GAL_TYPE_FLOAT64: COLLAPSE_DIM( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: type value (%d) not recognized",
+ __func__, in->type);
+ }
+ else
+ {
+ iif=(ii=num->array)+num->size;
+ do *ii++ = cnum; while(ii<iif);
+ }
+
+ /* Remove the respective dimension in the WCS structure also (if any
+ exists). Note that `sum->ndim' has already been changed. So we'll use
+ `in->wcs'. */
+ gal_wcs_remove_dimension(num->wcs, in->ndim-c_dim);
+
+ /* Return. */
+ if(wht!=weight) gal_data_free(wht);
+ return num;
+}
diff --git a/lib/gnuastro/dimension.h b/lib/gnuastro/dimension.h
index 5c83a0d..a0623a2 100644
--- a/lib/gnuastro/dimension.h
+++ b/lib/gnuastro/dimension.h
@@ -97,6 +97,23 @@ gal_dimension_dist_manhattan(size_t *a, size_t *b, size_t
ndim);
/************************************************************************/
+/******************** Collapsing a dimension **********************/
+/************************************************************************/
+gal_data_t *
+gal_dimension_collapse_sum(gal_data_t *in, size_t c_dim, gal_data_t *weight);
+
+gal_data_t *
+gal_dimension_collapse_mean(gal_data_t *in, size_t c_dim,
+ gal_data_t *weight);
+
+gal_data_t *
+gal_dimension_collapse_number(gal_data_t *in, size_t c_dim);
+
+
+
+
+
+/************************************************************************/
/******************** Neighbors **********************/
/************************************************************************/
/* Purpose
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 72be5a0..13c2bb5 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -74,6 +74,9 @@ struct wcsprm *
gal_wcs_copy(struct wcsprm *wcs);
void
+gal_wcs_remove_dimension(struct wcsprm *wcs, size_t fitsdim);
+
+void
gal_wcs_on_tile(gal_data_t *tile);
double *
diff --git a/lib/wcs.c b/lib/wcs.c
index c7d7376..3fe19cf 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -258,6 +258,176 @@ gal_wcs_copy(struct wcsprm *wcs)
+/* Remove the algorithm part of CTYPE (anything after, and including, a
+ `-') if necessary. */
+static void
+wcs_ctype_noalgorithm(char *str)
+{
+ size_t i, len=strlen(str);
+ for(i=0;i<len;++i) if(str[i]=='-') { str[i]='\0'; break; }
+}
+
+
+
+
+/* See if the CTYPE string ends with TAN. */
+static int
+wcs_ctype_has_tan(char *str)
+{
+ size_t len=strlen(str);
+
+ return !strcmp(&str[len-3], "TAN");
+}
+
+
+
+
+
+/* Remove dimension. */
+#define WCS_REMOVE_DIM_CHECK 0
+void
+gal_wcs_remove_dimension(struct wcsprm *wcs, size_t fitsdim)
+{
+ size_t c, i, j, naxis=wcs->naxis;
+
+ /* Sanity check. */
+ if(fitsdim==0 || fitsdim>wcs->naxis)
+ error(EXIT_FAILURE, 0, "%s: requested dimension (fitsdim=%zu) must be "
+ "larger than zero and smaller than the number of dimensions in "
+ "the given WCS structure (%zu)", __func__, fitsdim, naxis);
+
+ /* If the WCS structure is NULL, just return. */
+ if(wcs==NULL) return;
+
+ /**************************************************/
+#if WCS_REMOVE_DIM_CHECK
+ printf("\n\nfitsdim: %zu\n", fitsdim);
+ printf("\n##################\n");
+ /*
+ wcs->pc[0]=0; wcs->pc[1]=1; wcs->pc[2]=2;
+ wcs->pc[3]=3; wcs->pc[4]=4; wcs->pc[5]=5;
+ wcs->pc[6]=6; wcs->pc[7]=7; wcs->pc[8]=8;
+ */
+ for(i=0;i<wcs->naxis;++i)
+ {
+ for(j=0;j<wcs->naxis;++j)
+ printf("%-5g", wcs->pc[i*wcs->naxis+j]);
+ printf("\n");
+ }
+#endif
+ /**************************************************/
+
+
+ /* First loop over the arrays. */
+ for(i=0;i<naxis;++i)
+ {
+ /* The dimensions are in FITS order, but counting starts from 0, so
+ we'll have to subtract 1 from `fitsdim'. */
+ if(i>fitsdim-1)
+ {
+ /* 1-D arrays. */
+ if(wcs->crpix) wcs->crpix[i-1] = wcs->crpix[i];
+ if(wcs->cdelt) wcs->cdelt[i-1] = wcs->cdelt[i];
+ if(wcs->crval) wcs->crval[i-1] = wcs->crval[i];
+ if(wcs->crota) wcs->crota[i-1] = wcs->crota[i];
+ if(wcs->crder) wcs->crder[i-1] = wcs->crder[i];
+ if(wcs->csyer) wcs->csyer[i-1] = wcs->csyer[i];
+
+ /* The strings are all statically allocated, so we don't need to
+ check. */
+ memcpy(wcs->cunit[i-1], wcs->cunit[i], 72);
+ memcpy(wcs->ctype[i-1], wcs->ctype[i], 72);
+ memcpy(wcs->cname[i-1], wcs->cname[i], 72);
+
+ /* For 2-D arrays, just bring up all the rows. We'll fix the
+ columns in a second loop. */
+ for(j=0;j<naxis;++j)
+ {
+ if(wcs->pc) wcs->pc[ (i-1)*naxis+j ] = wcs->pc[ i*naxis+j ];
+ if(wcs->cd) wcs->cd[ (i-1)*naxis+j ] = wcs->cd[ i*naxis+j ];
+ }
+ }
+ }
+
+
+ /**************************************************/
+#if WCS_REMOVE_DIM_CHECK
+ printf("\n###### Respective row removed (replaced).\n");
+ for(i=0;i<wcs->naxis;++i)
+ {
+ for(j=0;j<wcs->naxis;++j)
+ printf("%-5g", wcs->pc[i*wcs->naxis+j]);
+ printf("\n");
+ }
+#endif
+ /**************************************************/
+
+
+ /* Second loop for 2D arrays. */
+ c=0;
+ for(i=0;i<naxis;++i)
+ for(j=0;j<naxis;++j)
+ if(j!=fitsdim-1)
+ {
+ if(wcs->pc) wcs->pc[ c ] = wcs->pc[ i*naxis+j ];
+ if(wcs->cd) wcs->cd[ c ] = wcs->cd[ i*naxis+j ];
+ ++c;
+ }
+
+
+ /* Correct the total number of dimensions in the WCS structure. */
+ naxis = wcs->naxis -= 1;
+
+
+ /* The `TAN' algorithm needs two dimensions. So we need to remove it when
+ it can cause confusion. */
+ switch(naxis)
+ {
+ /* The `TAN' algorithm cannot be used for any single-dimensional
+ dataset. So we'll have to remove it if it exists. */
+ case 1:
+ wcs_ctype_noalgorithm(wcs->ctype[0]);
+ break;
+
+ /* For any other dimensionality, `TAN' should be kept only when exactly
+ two dimensions have it. */
+ default:
+
+ c=0;
+ for(i=0;i<naxis;++i)
+ if( wcs_ctype_has_tan(wcs->ctype[i]) )
+ ++c;
+
+ if(c!=2)
+ for(i=0;i<naxis;++i)
+ if( wcs_ctype_has_tan(wcs->ctype[i]) )
+ wcs_ctype_noalgorithm(wcs->ctype[i]);
+ break;
+ }
+
+
+
+ /**************************************************/
+#if WCS_REMOVE_DIM_CHECK
+ printf("\n###### Respective column removed.\n");
+ for(i=0;i<naxis;++i)
+ {
+ for(j=0;j<naxis;++j)
+ printf("%-5g", wcs->pc[i*naxis+j]);
+ printf("\n");
+ }
+ printf("\n###### One final string\n");
+ for(i=0;i<naxis;++i)
+ printf("%s\n", wcs->ctype[i]);
+ exit(0);
+#endif
+ /**************************************************/
+}
+
+
+
+
+
/* Using the block data structure of the tile, add a WCS structure for
it. In many cases, tiles are created for internal processing, so there
is no need to keep their WCS. Hence for preformance reasons, when